diff --git a/surface_integral/interface/mesh_algorithm.hpp b/surface_integral/interface/mesh_algorithm.hpp index aacf571..a2e52e3 100644 --- a/surface_integral/interface/mesh_algorithm.hpp +++ b/surface_integral/interface/mesh_algorithm.hpp @@ -1,104 +1,96 @@ #include #include -#include // polymesh_t - +#include // polymesh_t // Vector subtraction: a - b -inline vector3d vec3_sub(const vector3d& a, const vector3d& b) { - return { a.x - b.x, a.y - b.y, a.z - b.z }; -} +inline vector3d vec3_sub(const vector3d& a, const vector3d& b) { return {a.x - b.x, a.y - b.y, a.z - b.z}; } // Vector cross product: a × b -inline vector3d vec3_cross(const vector3d& a, const vector3d& b) { - return { - a.y * b.z - a.z * b.y, - a.z * b.x - a.x * b.z, - a.x * b.y - a.y * b.x - }; +inline vector3d vec3_cross(const vector3d& a, const vector3d& b) +{ + return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; } // Vector dot product: a · b -inline double vec3_dot(const vector3d& a, const vector3d& b) { - return a.x * b.x + a.y * b.y + a.z * b.z; -} +inline double vec3_dot(const vector3d& a, const vector3d& b) { return a.x * b.x + a.y * b.y + a.z * b.z; } // Vector length: ||a|| -inline double vec3_norm(const vector3d& a) { - return std::sqrt(a.x*a.x + a.y*a.y + a.z*a.z); -} +inline double vec3_norm(const vector3d& a) { return std::sqrt(a.x * a.x + a.y * a.y + a.z * a.z); } // Compute triangle area: 0.5 * || (b-a) × (c-a) || -inline double triangle_area(const vector3d& a, const vector3d& b, const vector3d& c) { - vector3d ab = vec3_sub(b, a); - vector3d ac = vec3_sub(c, a); +inline double triangle_area(const vector3d& a, const vector3d& b, const vector3d& c) +{ + vector3d ab = vec3_sub(b, a); + vector3d ac = vec3_sub(c, a); vector3d cross_product = vec3_cross(ab, ac); return 0.5 * vec3_norm(cross_product); } // Compute triangle signed volume contribution (w.r.t origin): (a · (b × c)) / 6 -inline double triangle_signed_volume(const vector3d& a, const vector3d& b, const vector3d& c) { +inline double triangle_signed_volume(const vector3d& a, const vector3d& b, const vector3d& c) +{ vector3d cross_product = vec3_cross(b, c); return vec3_dot(a, cross_product) / 6.0; } // Compute polygon face area (triangulate as fan) -double polygon_area(const vector3d* vertices, const uint32_t* face_indices, uint32_t n) { +double polygon_area(const vector3d* vertices, const uint32_t* face_indices, uint32_t n) +{ if (n < 3) return 0.0; - double area = 0.0; - const vector3d& v0 = vertices[face_indices[0]]; + double area = 0.0; + const vector3d& v0 = vertices[face_indices[0]]; for (uint32_t i = 1; i < n - 1; ++i) { - const vector3d& v1 = vertices[face_indices[i]]; - const vector3d& v2 = vertices[face_indices[i + 1]]; - area += triangle_area(v0, v1, v2); + const vector3d& v1 = vertices[face_indices[i]]; + const vector3d& v2 = vertices[face_indices[i + 1]]; + area += triangle_area(v0, v1, v2); } return area; } // Compute polygon signed volume contribution (sum of triangle contributions) -double polygon_signed_volume(const vector3d* vertices, const uint32_t* face_indices, uint32_t n) { +double polygon_signed_volume(const vector3d* vertices, const uint32_t* face_indices, uint32_t n) +{ if (n < 3) return 0.0; - double volume = 0.0; - const vector3d& v0 = vertices[face_indices[0]]; + double volume = 0.0; + const vector3d& v0 = vertices[face_indices[0]]; for (uint32_t i = 1; i < n - 1; ++i) { - const vector3d& v1 = vertices[face_indices[i]]; - const vector3d& v2 = vertices[face_indices[i + 1]]; - volume += triangle_signed_volume(v0, v1, v2); + const vector3d& v1 = vertices[face_indices[i]]; + const vector3d& v2 = vertices[face_indices[i + 1]]; + volume += triangle_signed_volume(v0, v1, v2); } return volume; } // 🟩 Main: compute total mesh surface area -double compute_surface_area(const polymesh_t* mesh) { - if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) { - return 0.0; - } +double compute_surface_area(const polymesh_t* mesh) +{ + if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) { return 0.0; } - double total_area = 0.0; - const uint32_t* face_ptr = mesh->faces; // current index into faces array + double total_area = 0.0; + const uint32_t* face_ptr = mesh->faces; // current index into faces array for (uint32_t i = 0; i < mesh->num_faces; ++i) { - uint32_t n = mesh->vertex_counts[i]; + uint32_t n = mesh->vertex_counts[i]; total_area += polygon_area(mesh->vertices, face_ptr, n); - face_ptr += n; // move to the next face's start index + face_ptr += n; // move to the next face's start index } return total_area; } // 🟦 Main: compute total mesh volume (requires closed mesh with outward normals) -double compute_volume(const polymesh_t* mesh) { - if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) { - return 0.0; - } +double compute_volume(const polymesh_t* mesh) +{ + if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) { return 0.0; } - double total_volume = 0.0; - const uint32_t* face_ptr = mesh->faces; + double total_volume = 0.0; + const uint32_t* face_ptr = mesh->faces; for (uint32_t i = 0; i < mesh->num_faces; ++i) { - uint32_t n = mesh->vertex_counts[i]; + uint32_t n = mesh->vertex_counts[i]; total_volume += polygon_signed_volume(mesh->vertices, face_ptr, n); - face_ptr += n; + face_ptr += n; } - return std::abs(total_volume); // ensure positive volume + return std::abs(total_volume); // ensure positive volume } \ No newline at end of file diff --git a/surface_integral/interface/subface_integrator.hpp b/surface_integral/interface/subface_integrator.hpp index 4ff62a3..8bb3b9c 100644 --- a/surface_integral/interface/subface_integrator.hpp +++ b/surface_integral/interface/subface_integrator.hpp @@ -62,13 +62,12 @@ private: stl_vector_mp compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max) const; void find_v_intersections_at_u(const subface& subface, - const parametric_plane_t& param_plane, - double u_val, - double v_min, - double v_max, - stl_vector_mp& intersections, - stl_vector_mp& intersected_chains - ) const; + const parametric_plane_t& param_plane, + double u_val, + double v_min, + double v_max, + stl_vector_mp& intersections, + stl_vector_mp& intersected_chains) const; bool is_point_inside_domain(const subface& subface, const parametric_plane_t& param_plane, double u, double v) const; diff --git a/surface_integral/src/polygon_winding_number.cpp b/surface_integral/src/polygon_winding_number.cpp index 7e41073..2e07715 100644 --- a/surface_integral/src/polygon_winding_number.cpp +++ b/surface_integral/src/polygon_winding_number.cpp @@ -16,29 +16,24 @@ const double eps = 1e-8; * @param B End point of the edge. * @return +1 if upward crossing, -1 if downward, 0 otherwise. */ -int winding_contribution(const Point2d& P, const Point2d& A, const Point2d& B) { +int winding_contribution(const Point2d& P, const Point2d& A, const Point2d& B) +{ // 1. If the edge is horizontal, no contribution - if (std::abs(A.y() - B.y()) < eps) { - return 0; - } + if (std::abs(A.y() - B.y()) < eps) { return 0; } // 2. Determine edge direction: does it cross P.y from below to above or vice versa? bool a_below = (A.y() < P.y() - eps); bool b_below = (B.y() < P.y() - eps); // 3. If both endpoints are on the same side of the ray, no crossing - if (a_below == b_below) { - return 0; - } + if (a_below == b_below) { return 0; } // 4. Compute x-coordinate of intersection with y = P.y - double t = (P.y() - A.y()) / (B.y() - A.y()); + double t = (P.y() - A.y()) / (B.y() - A.y()); double x_intersect = A.x() + t * (B.x() - A.x()); // 5. Only consider intersections to the right of P (or directly on P) - if (x_intersect < P.x() - eps) { - return 0; - } + if (x_intersect < P.x() - eps) { return 0; } // 6. Upward crossing: from below to above → +1 if (a_below && !b_below) { @@ -55,31 +50,31 @@ int winding_contribution(const Point2d& P, const Point2d& A, const Point2d& B) { /** * @brief Computes the total winding number of point P with respect to a closed ring. */ -int compute_winding_number_ring(const Point2d& P, const PolygonRing& ring) { +int compute_winding_number_ring(const Point2d& P, const PolygonRing& ring) +{ if (ring.size() < 3) return 0; int wn = 0; - int n = ring.size(); + int n = ring.size(); for (int i = 0; i < n; ++i) { - const Point2d& A = ring[i]; - const Point2d& B = ring[(i + 1) % n]; - wn += winding_contribution(P, A, B); + const Point2d& A = ring[i]; + const Point2d& B = ring[(i + 1) % n]; + wn += winding_contribution(P, A, B); } return wn; } // Helper: Check if P is on segment A-B -bool is_point_on_segment(const Point2d& P, const Point2d& A, const Point2d& B) { - if (P.x() < std::min(A.x(), B.x()) - eps || - P.x() > std::max(A.x(), B.x()) + eps || - P.y() < std::min(A.y(), B.y()) - eps || - P.y() > std::max(A.y(), B.y()) + eps) { +bool is_point_on_segment(const Point2d& P, const Point2d& A, const Point2d& B) +{ + if (P.x() < std::min(A.x(), B.x()) - eps || P.x() > std::max(A.x(), B.x()) + eps || P.y() < std::min(A.y(), B.y()) - eps + || P.y() > std::max(A.y(), B.y()) + eps) { return false; } double cross = (P.y() - A.y()) * (B.x() - A.x()) - (P.x() - A.x()) * (B.y() - A.y()); if (std::abs(cross) > eps) return false; - double dot = (P.x() - A.x()) * (B.x() - A.x()) + (P.y() - A.y()) * (B.y() - A.y()); + double dot = (P.x() - A.x()) * (B.x() - A.x()) + (P.y() - A.y()) * (B.y() - A.y()); double len_sq = (B.x() - A.x()) * (B.x() - A.x()) + (B.y() - A.y()) * (B.y() - A.y()); return (dot >= -eps && dot <= len_sq + eps); } @@ -91,18 +86,13 @@ bool is_point_on_segment(const Point2d& P, const Point2d& A, const Point2d& B) { * - Points on hole boundary → outside * - Otherwise: wn(outer) + sum(wn(hole)) != 0 → inside */ -bool point_in_polygon_with_holes( - const Point2d& P, - const PolygonRing& outer_ring, - const std::vector& holes) +bool point_in_polygon_with_holes(const Point2d& P, const PolygonRing& outer_ring, const std::vector& holes) { // Check boundary first auto on_ring = [&](const PolygonRing& ring) { int n = ring.size(); for (int i = 0; i < n; ++i) { - if (is_point_on_segment(P, ring[i], ring[(i+1)%n])) { - return true; - } + if (is_point_on_segment(P, ring[i], ring[(i + 1) % n])) { return true; } } return false; }; @@ -118,11 +108,9 @@ bool point_in_polygon_with_holes( total_wn += compute_winding_number_ring(P, hole); // hole 是 CW，返回负值 } - return total_wn != 0; // 非零即内部 + return total_wn != 0; // 非零即内部 } - - /** * @brief Validates the polygon input for correct orientation and sufficient vertices. * @@ -132,71 +120,69 @@ bool point_in_polygon_with_holes( * @param context Optional context string for logging. * @return true if valid, false otherwise. */ -bool validate_polygon_input( - const Point2d& query_point, - const PolygonRing& outer_ring, - const std::vector& holes, - const std::string& context) +bool validate_polygon_input(const Point2d& query_point, + const PolygonRing& outer_ring, + const std::vector& holes, + const std::string& context) { - std::string label = context.empty() ? "Polygon" : context; - bool valid = true; + std::string label = context.empty() ? "Polygon" : context; + bool valid = true; std::string issues = ""; // Check outer ring has enough points if (outer_ring.size() < 3) { issues += "Outer ring has fewer than 3 vertices."; - valid = false; + valid = false; } else { double outer_area = compute_ring_area(outer_ring); if (outer_area <= 0) { issues += "Outer ring is not CCW (area = " + std::to_string(outer_area) + ")"; - valid = false; + valid = false; } } // Check each hole for (size_t i = 0; i < holes.size(); ++i) { - const auto& hole = holes[i]; + const auto& hole = holes[i]; std::string hole_label = "Hole " + std::to_string(i); if (hole.size() < 3) { if (!issues.empty()) issues += ", "; issues += hole_label + " has fewer than 3 vertices"; - valid = false; + valid = false; } else { double hole_area = compute_ring_area(hole); if (hole_area >= 0) { if (!issues.empty()) issues += ", "; issues += hole_label + " is not CW (area = " + std::to_string(hole_area) + ")"; - valid = false; + valid = false; } } } // Output issues if invalid if (!valid) { - std::cerr << "[\033[33mVALIDATION ERROR\033[0m] " << label - << " | Query: " << query_point - << " | Issues: " << issues << "\n"; + std::cerr << "[\033[33mVALIDATION ERROR\033[0m] " << label << " | Query: " << query_point << " | Issues: " << issues + << "\n"; } return valid; } - /** * @brief Computes the signed area of a polygon ring. * * @param ring The polygon ring (closed). * @return The signed area (positive if CCW, negative if CW). */ -double compute_ring_area(const PolygonRing& ring) { +double compute_ring_area(const PolygonRing& ring) +{ if (ring.size() < 3) return 0.0; double area = 0.0; - int n = ring.size(); + int n = ring.size(); for (int i = 0; i < n; ++i) { - int j = (i + 1) % n; - area += ring[i].x() * ring[j].y() - ring[j].x() * ring[i].y(); + int j = (i + 1) % n; + area += ring[i].x() * ring[j].y() - ring[j].x() * ring[i].y(); } return area * 0.5; } \ No newline at end of file diff --git a/surface_integral/src/subface_integrator.cpp b/surface_integral/src/subface_integrator.cpp index 85b9d66..fc3d60c 100644 --- a/surface_integral/src/subface_integrator.cpp +++ b/surface_integral/src/subface_integrator.cpp @@ -9,7 +9,7 @@ namespace internal { integrator_t(const stl_vector_mp>& surfaces, - const flat_hash_map_mp& uv_planes) + const flat_hash_map_mp& uv_planes) : m_subfaces(surfaces), m_uv_planes(uv_planes) { std::cout << "Integrator initialized with " << surfaces.size() << " subfaces." << std::endl; @@ -19,17 +19,15 @@ integrator_t(const stl_vector_mp>& surfaces, for (int subface_index = 0; subface_index < m_subfaces.size(); ++subface_index) { if (m_uv_planes.find(subface_index) != m_uv_planes.end()) { const auto& param_plane = m_uv_planes.at(subface_index); - auto& chain_bboxs = m_chain_bboxes_hash[subface_index]; + auto& chain_bboxs = m_chain_bboxes_hash[subface_index]; chain_bboxes.clear(); chain_bboxs.reverse(chain_vertices.size()); - for (const auto& chain: param_plane.chain_vertices){ + for (const auto& chain : param_plane.chain_vertices) { Eigen::AlignedBox2d bbox; - bbox.setEmpty(); - + bbox.setEmpty(); + if (!chain.empty()) { - for (const auto& pt : chain) { - bbox.extend(pt); - } + for (const auto& pt : chain) { bbox.extend(pt); } } chain_bboxes.push_back(bbox); } @@ -46,9 +44,9 @@ double integrator_t::calculate( double total_integral = 0.0; for (int subface_index = 0; subface_index < m_subfaces.size(); ++subface_index) { if (m_uv_planes.find(subface_index) != m_uv_planes.end()) { - const auto& param_plane = m_uv_planes.at(subface_index); - const auto& subface = m_subfaces[subface_index].object_ptr.get(); - total_integral += calculate_one_subface(subface, param_plane, gauss_order, func); + const auto& param_plane = m_uv_planes.at(subface_index); + const auto& subface = m_subfaces[subface_index].object_ptr.get(); + total_integral += calculate_one_subface(subface, param_plane, gauss_order, func); } else { // the subface has no associated parametric plane, meaning it is compete face without trimming } @@ -62,16 +60,16 @@ double integrator_t::calculate_one_subface( int gauss_order, double (*func)(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)) const { - auto solver = subface.fetch_solver_evaluator(); - const auto& uv_bounds = param_plane.uv_bounds; - auto u_min = uv_bounds.min().x(); - auto u_max = uv_bounds.max().x(); - auto v_min = uv_bounds.min().y(); - auto v_max = uv_bounds.max().y(); + auto solver = subface.fetch_solver_evaluator(); + const auto& uv_bounds = param_plane.uv_bounds; + auto u_min = uv_bounds.min().x(); + auto u_max = uv_bounds.max().x(); + auto v_min = uv_bounds.min().y(); + auto v_max = uv_bounds.max().y(); // Gaussian integration in u direction - auto u_integrand = [&](double u) { + auto u_integrand = [&](double u) { // Find exact v intersections for each u - stl_vector_mp v_intersections; + stl_vector_mp v_intersections; stl_vector_mp intersected_chains; find_v_intersections_at_u(subface, param_plane, u, v_min, v_max, v_intersections, intersected_chains); @@ -123,8 +121,8 @@ double integrator_t::calculate_one_subface( double a = u_breaks[i]; double b = u_breaks[i + 1]; - double mid_u = (a + b) / 2.0; - stl_vector_mp v_intersections; + double mid_u = (a + b) / 2.0; + stl_vector_mp v_intersections; stl_vector_mp intersected_chains; find_v_intersections_at_u(subface, param_plane, mid_u, v_min, v_max, v_intersections, intersected_chains); @@ -165,40 +163,39 @@ stl_vector_mp integrator_t::compute_u_breaks(const parametric_plane_t& p return stl_vector_mp(break_set.begin(), break_set.end()); } -void integrator_t::find_v_intersections_at_u(const subface& subface, - const parametric_plane_t& param_plane, - const stl_vector_mp& chain_bboxs, - double u_val, - double v_min, - double v_max, - stl_vector_mp& intersections, - stl_vector_mp& intersected_chains - ) const +void integrator_t::find_v_intersections_at_u(const subface& subface, + const parametric_plane_t& param_plane, + const stl_vector_mp& chain_bboxs, + double u_val, + double v_min, + double v_max, + stl_vector_mp& intersections, + stl_vector_mp& intersected_chains) const { intersections.clear(); intersected_chains.clear(); - uint16_t chain_idx =0; + uint16_t chain_idx = 0; // Iterate over each boundary chain for (const auto& chain : param_plane.chain_vertices) { - const size_t n_vertices = chain.size()-1; - if (n_vertices < 2){ + const size_t n_vertices = chain.size() - 1; + if (n_vertices < 2) { std::cout << "Warning: Degenerate chain with less than 2 vertices." << std::endl; chain_idx++; continue; // Skip degenerate chains } // filter bbox - if (chain_bbox[chain_idx].isEmpty() - || (chain_bbox[chain_idx].x()u_val+1e-8)) { + if (chain_bbox[chain_idx].isEmpty() + || (chain_bbox[chain_idx].x() < u_val - 1e-8 || chain_bbox[chain_idx].x() > u_val + 1e-8)) { chain_idx++; continue; } // Iterate over each edge in the chain (including closing edge if desired) for (size_t i = 0; i < n_vertices; ++i) { - const Eigen::Vector2d& v1 = chain[i]; // Current vertex: (u1, v1) + const Eigen::Vector2d& v1 = chain[i]; // Current vertex: (u1, v1) const Eigen::Vector2d& v2 = chain[i + 1]; // Next vertex: (u2, v2) double u1 = v1.x(), v1_val = v1.y(); @@ -264,29 +261,23 @@ void integrator_t::find_v_intersections_at_u(const subface& subface, } } } - chain_idx ++; + chain_idx++; } // Final step: sort and remove duplicates within tolerance - //if (!intersections.empty()) { sort_and_unique_with_tol(intersections, 1e-8); } - + // if (!intersections.empty()) { sort_and_unique_with_tol(intersections, 1e-8); } } - - -bool is_edge_inside_domain(const stl_vector_mp> & chain_group_indices, - uint16_t chain_idx1, - uint16_t chain_idx2) const +bool is_edge_inside_domain(const stl_vector_mp>& chain_group_indices, + uint16_t chain_idx1, + uint16_t chain_idx2) const { const auto& group1 = chain_group_indices[chain_idx1]; const auto& group2 = chain_group_indices[chain_idx2]; - if (group1.find(chain_idx2) != group1.end() && group2.find(chain_idx1) != group2.end()) { - return true; - } + if (group1.find(chain_idx2) != group1.end() && group2.find(chain_idx1) != group2.end()) { return true; } return false; } - bool integrator_t::is_u_near_singularity(double u, double tol) const { return false; } void integrator_t::sort_and_unique_with_tol(stl_vector_mp& vec, double epsilon) const diff --git a/surface_integral/test/SurfaceIntegrator_test.cpp b/surface_integral/test/SurfaceIntegrator_test.cpp index 181536e..327679b 100644 --- a/surface_integral/test/SurfaceIntegrator_test.cpp +++ b/surface_integral/test/SurfaceIntegrator_test.cpp @@ -5,13 +5,14 @@ #include // ✅ 封装成函数：通过引用填充一个 parametric_plane_t 实例 -void populate_parametric_plane(parametric_plane_t& plane) { +void populate_parametric_plane(parametric_plane_t& plane) +{ // 确保是空的（可选） // plane = parametric_plane_t(); - auto& chain_vertices = plane.chain_vertices; + auto& chain_vertices = plane.chain_vertices; auto& chain_group_indices = plane.chain_group_indices; - auto& v_flags = plane.vertex_special_flags; - auto& e_flags = plane.edge_near_parallel_flags; + auto& v_flags = plane.vertex_special_flags; + auto& e_flags = plane.edge_near_parallel_flags; int vertex_count = 4; @@ -20,10 +21,18 @@ void populate_parametric_plane(parametric_plane_t& plane) { e_flags.reserve(vertex_count - 1); // 1. 添加一条链的顶点 (chain 0) - chain_vertices.push_back({{0.0, 0.0}}); - chain_vertices.push_back({{1.0, 0.0}}); - chain_vertices.push_back({{1.0, 1.0}}); - chain0_vertices.push_back({{0.0, 1.0}}); + chain_vertices.push_back({ + {0.0, 0.0} + }); + chain_vertices.push_back({ + {1.0, 0.0} + }); + chain_vertices.push_back({ + {1.0, 1.0} + }); + chain0_vertices.push_back({ + {0.0, 1.0} + }); chain_vertices.reserve(unique_chain_indices.size()); @@ -45,7 +54,7 @@ void populate_parametric_plane(parametric_plane_t& plane) { } v_flags.reverse(); - + dynamic_bitset_mp<> v_flags; v_flags.set(0); v_flags.set(3); @@ -64,7 +73,8 @@ void populate_parametric_plane(parametric_plane_t& plane) { } // --- 使用示例 --- -int main() { +int main() +{ std::cout << "=== 使用函数填充 parametric_plane_t 实例 ===" << std::endl; // 先创建实例 @@ -81,9 +91,7 @@ int main() { std::cout << "第一条链的顶点特殊标志: " << plane.vertex_special_flags[0] << std::endl; std::cout << "第一条链的边平行标志: " << plane.edge_near_parallel_flags[0] << std::endl; std::cout << "第一条链所属的组: "; - for (auto idx : plane.chain_group_indices[0]) { - std::cout << idx << " "; - } + for (auto idx : plane.chain_group_indices[0]) { std::cout << idx << " "; } std::cout << std::endl; } @@ -98,16 +106,15 @@ int main() internal::sphere_face_t sphere_face; stl_vector_mp> subfaces; - flat_hash_map_mp uv_planes; + flat_hash_map_mp uv_planes; subfaces.reserve(1); subfaces.emplace_back(&sphere_face, 0); - uv_planes[0] = parametric_plane_t{}; + uv_planes[0] = parametric_plane_t{}; uv_planes[0].u_min = 0.0; - // box_descriptor_t box{ // {0., 0., 0.}, // {1., 1., 1.} diff --git a/surface_integral/test/test_mesh.cpp b/surface_integral/test/test_mesh.cpp index 204eee7..f64ec4e 100644 --- a/surface_integral/test/test_mesh.cpp +++ b/surface_integral/test/test_mesh.cpp @@ -1,37 +1,44 @@ #include #include -int main() { +int main() +{ // 8 vertices: unit cube vector3d verts[8] = { - {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0}, - {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1} + {0, 0, 0}, + {1, 0, 0}, + {1, 1, 0}, + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, + {1, 1, 1}, + {0, 1, 1} }; // 6 faces, each with 4 vertices (counter-clockwise order when viewed from outside) uint32_t face_data[] = { - 0,1,2,3, // bottom face - 7,6,5,4, // top face (note CCW when seen from outside) - 0,4,5,1, // front face - 1,5,6,2, // right face - 3,2,6,7, // back face - 0,3,7,4 // left face + 0, 1, 2, 3, // bottom face + 7, 6, 5, 4, // top face (note CCW when seen from outside) + 0, 4, 5, 1, // front face + 1, 5, 6, 2, // right face + 3, 2, 6, 7, // back face + 0, 3, 7, 4 // left face }; uint32_t vertex_counts[6] = {4, 4, 4, 4, 4, 4}; polymesh_t mesh; - mesh.vertices = verts; - mesh.faces = face_data; + mesh.vertices = verts; + mesh.faces = face_data; mesh.vertex_counts = vertex_counts; - mesh.num_vertices = 8; - mesh.num_faces = 6; + mesh.num_vertices = 8; + mesh.num_faces = 6; - double area = compute_surface_area(&mesh); + double area = compute_surface_area(&mesh); double volume = compute_volume(&mesh); - std::cout << "area: " << area << std::endl; // expected: 6 - std::cout << "volume: " << volume << std::endl; // expected: 1 + std::cout << "area: " << area << std::endl; // expected: 6 + std::cout << "volume: " << volume << std::endl; // expected: 1 return 0; } \ No newline at end of file diff --git a/surface_integral/test/test_winding_number.cpp b/surface_integral/test/test_winding_number.cpp index 05be6c0..c53cf3e 100644 --- a/surface_integral/test/test_winding_number.cpp +++ b/surface_integral/test/test_winding_number.cpp @@ -10,35 +10,42 @@ // ----------------------------- struct TestCase { - std::string name; - PolygonRing outer_ring; + std::string name; + PolygonRing outer_ring; std::vector holes; - Point2d query_point; - bool expected; // true = inside, false = outside + Point2d query_point; + bool expected; // true = inside, false = outside }; // ----------------------------- // Helper: Print Point // ----------------------------- -std::ostream& operator<<(std::ostream& os, const Point2d& p) { - return os << "(" << p.x() << ", " << p.y() << ")"; -} - - +std::ostream& operator<<(std::ostream& os, const Point2d& p) { return os << "(" << p.x() << ", " << p.y() << ")"; } // ----------------------------- // Helper: Validate Ring Orientation // ----------------------------- -bool validate_orientation() { +bool validate_orientation() +{ // Test CCW (positive area) - PolygonRing ccw = {{0,0}, {1,0}, {1,1}, {0,1}}; + PolygonRing ccw = { + {0, 0}, + {1, 0}, + {1, 1}, + {0, 1} + }; if (compute_ring_area(ccw) <= 0) { std::cerr << "Error: CCW ring has non-positive area!\n"; return false; } // Test CW (negative area) - PolygonRing cw = {{0,0}, {0,1}, {1,1}, {1,0}}; + PolygonRing cw = { + {0, 0}, + {0, 1}, + {1, 1}, + {1, 0} + }; if (compute_ring_area(cw) >= 0) { std::cerr << "Error: CW ring has non-negative area!\n"; return false; @@ -49,15 +56,16 @@ bool validate_orientation() { // ----------------------------- // Validate Test Case: Orientation of outer ring (CCW) and holes (CW) // ----------------------------- -bool validate_test_case(const TestCase& tc) { - bool valid = true; +bool validate_test_case(const TestCase& tc) +{ + bool valid = true; std::string status = ""; // Validate outer ring orientation (should be CCW → positive area) double outer_area = compute_ring_area(tc.outer_ring); if (outer_area <= 0) { status += "Outer ring not CCW (area = " + std::to_string(outer_area) + ")"; - valid = false; + valid = false; } // Validate holes orientation (should be CW → negative area) @@ -66,18 +74,16 @@ bool validate_test_case(const TestCase& tc) { if (hole_area >= 0) { if (!status.empty()) status += ", "; status += "Hole " + std::to_string(i) + " not CW (area = " + std::to_string(hole_area) + ")"; - valid = false; + valid = false; } } if (valid) return true; std::cout << "[" - << "\033[33mINVALID\033[0m" // make "INVALID" yellow - << "] " - << tc.name << "\n" - << " Query: " << tc.query_point - << " | Status: \033[33mInvalid Test Case\033[0m" + << "\033[33mINVALID\033[0m" // make "INVALID" yellow + << "] " << tc.name << "\n" + << " Query: " << tc.query_point << " | Status: \033[33mInvalid Test Case\033[0m" << " | Reason: " << (valid ? "OK" : status) << "\n"; return valid; @@ -86,16 +92,13 @@ bool validate_test_case(const TestCase& tc) { // ----------------------------- // Run Single Test // ----------------------------- -bool run_test(const TestCase& tc) { +bool run_test(const TestCase& tc) +{ bool result = point_in_polygon_with_holes(tc.query_point, tc.outer_ring, tc.holes); bool passed = (result == tc.expected); - std::cout << "[" - << (passed ? "\033[32mPASS\033[0m" : "\033[31mFAIL\033[0m") - << "] " - << tc.name << "\n" - << " Query: " << tc.query_point - << " | Expected: " << (tc.expected ? "Inside" : "Outside") + std::cout << "[" << (passed ? "\033[32mPASS\033[0m" : "\033[31mFAIL\033[0m") << "] " << tc.name << "\n" + << " Query: " << tc.query_point << " | Expected: " << (tc.expected ? "Inside" : "Outside") << " | Got: " << (result ? "Inside" : "Outside") << "\n"; return passed; @@ -105,19 +108,32 @@ bool run_test(const TestCase& tc) { // Test Case Factory Functions // ----------------------------- -std::vector create_test_cases() { +std::vector create_test_cases() +{ std::vector tests; - // -------------------------------------------------- // Test 1: Simple hole - point inside hole → Outside // -------------------------------------------------- - PolygonRing outer1 = {{0,0}, {2,0}, {2,2}, {0,2}}; - PolygonRing hole1 = {{0.5,0.5}, {0.5,1.5}, {1.5,1.5}, {1.5,0.5}}; + PolygonRing outer1 = { + {0, 0}, + {2, 0}, + {2, 2}, + {0, 2} + }; + PolygonRing hole1 = { + {0.5, 0.5}, + {0.5, 1.5}, + {1.5, 1.5}, + {1.5, 0.5} + }; tests.push_back({ "Inside Hole", - outer1, {hole1}, {1.0, 1.0}, false + outer1, + {hole1}, + {1.0, 1.0}, + false }); // -------------------------------------------------- @@ -125,7 +141,10 @@ std::vector create_test_cases() { // -------------------------------------------------- tests.push_back({ "Inside Outer, Outside Hole", - outer1, {hole1}, {0.1, 0.1}, true + outer1, + {hole1}, + {0.1, 0.1}, + true }); // -------------------------------------------------- @@ -133,7 +152,10 @@ std::vector create_test_cases() { // -------------------------------------------------- tests.push_back({ "Outside Outer", - outer1, {hole1}, {3.0, 3.0}, false + outer1, + {hole1}, + {3.0, 3.0}, + false }); // -------------------------------------------------- @@ -142,7 +164,10 @@ std::vector create_test_cases() { // -------------------------------------------------- tests.push_back({ "On Outer Edge (Bottom)", - outer1, {}, {1.0, 0.0}, true + outer1, + {}, + {1.0, 0.0}, + true }); // -------------------------------------------------- @@ -150,70 +175,138 @@ std::vector create_test_cases() { // -------------------------------------------------- tests.push_back({ "On Outer Vertex Offset Slightly Outside", - outer1, {}, {0.0, 0.0-1e-8}, false + outer1, + {}, + {0.0, 0.0 - 1e-8}, + false }); // -------------------------------------------------- // Test 6: Point on hole edge // Should be outside (on hole boundary still "outside" the solid) // -------------------------------------------------- - PolygonRing hole6 = {{1.0,1.0}, {1.0,2.0}, {2.0,2.0}, {2.0,1.0}}; // CW + PolygonRing hole6 = { + {1.0, 1.0}, + {1.0, 2.0}, + {2.0, 2.0}, + {2.0, 1.0} + }; // CW tests.push_back({ "On Hole Edge", - {{0,0},{3,0},{3,3},{0,3}}, {hole6}, {1.5, 1.0}, false + {{0, 0}, {3, 0}, {3, 3}, {0, 3}}, + {hole6}, + {1.5, 1.0}, + false }); // -------------------------------------------------- // Test 7: Multiple holes // -------------------------------------------------- - PolygonRing hole7a = {{0.5,0.5}, {0.5,1.0}, {1.0,1.0}, {1.0,0.5}}; // CW - PolygonRing hole7b = {{1.5,1.5}, {1.5,2.0}, {2.0,2.0}, {2.0,1.5}}; // CW + PolygonRing hole7a = { + {0.5, 0.5}, + {0.5, 1.0}, + {1.0, 1.0}, + {1.0, 0.5} + }; // CW + PolygonRing hole7b = { + {1.5, 1.5}, + {1.5, 2.0}, + {2.0, 2.0}, + {2.0, 1.5} + }; // CW tests.push_back({ "Multiple Holes - Inside", - outer1, {hole7a, hole7b}, {0.1, 0.1}, true + outer1, + {hole7a, hole7b}, + {0.1, 0.1 }, + true }); tests.push_back({ "Multiple Holes - In First Hole", - outer1, {hole7a, hole7b}, {0.7, 0.7}, false + outer1, + {hole7a, hole7b}, + {0.7, 0.7 }, + false }); tests.push_back({ "Multiple Holes - In Second Hole", - outer1, {hole7a, hole7b}, {1.7, 1.7}, false + outer1, + {hole7a, hole7b}, + {1.7, 1.7 }, + false }); // -------------------------------------------------- // Test 8: Nested holes (hole inside hole) - NOT standard, but test behavior // Note: This is not valid in simple polygons, but test robustness // -------------------------------------------------- - PolygonRing outer8 = {{0,0}, {3,0}, {3,3}, {0,3}}; - PolygonRing hole8a = {{0.5,0.5}, {0.5,2.5}, {2.5,2.5}, {2.5,0.5}}; // large hole - PolygonRing hole8b = {{1.0,1.0}, {1.0,2.0}, {2.0,2.0}, {2.0,1.0}}; // small hole inside hole8a + PolygonRing outer8 = { + {0, 0}, + {3, 0}, + {3, 3}, + {0, 3} + }; + PolygonRing hole8a = { + {0.5, 0.5}, + {0.5, 2.5}, + {2.5, 2.5}, + {2.5, 0.5} + }; // large hole + PolygonRing hole8b = { + {1.0, 1.0}, + {1.0, 2.0}, + {2.0, 2.0}, + {2.0, 1.0} + }; // small hole inside hole8a // This creates a "island" in the hole tests.push_back({ "Nested Holes - On Island", - outer8, {hole8a, hole8b}, {1.5, 1.5}, true // inside outer, inside hole8a, but outside hole8b → solid + outer8, + {hole8a, hole8b}, + {1.5, 1.5 }, + true // inside outer, inside hole8a, but outside hole8b → solid }); // -------------------------------------------------- // Test 9: Complex L-shaped polygon with hole // -------------------------------------------------- PolygonRing outer9 = { - {0,0}, {3,0}, {3,1}, {2,1}, {2,3}, {0,3}, {0,0} // L-shape, explicitly closed and no self-intersection + {0, 0}, + {3, 0}, + {3, 1}, + {2, 1}, + {2, 3}, + {0, 3}, + {0, 0} // L-shape, explicitly closed and no self-intersection }; // CCW - PolygonRing hole9 = {{0.5,0.5}, {0.5,1.5}, {1.5,1.5}, {1.5,0.5}}; // CW + PolygonRing hole9 = { + {0.5, 0.5}, + {0.5, 1.5}, + {1.5, 1.5}, + {1.5, 0.5} + }; // CW tests.push_back({ "L-Shape with Hole - Inside Leg", - outer9, {hole9}, {0.1, 0.1}, true + outer9, + {hole9}, + {0.1, 0.1}, + true }); tests.push_back({ "L-Shape with Hole - In Hole", - outer9, {hole9}, {1.0, 1.0}, false + outer9, + {hole9}, + {1.0, 1.0}, + false }); tests.push_back({ "L-Shape with Hole - In Arm", - outer9, {hole9}, {2.0, 0.5}, true + outer9, + {hole9}, + {2.0, 0.5}, + true }); // -------------------------------------------------- @@ -221,7 +314,10 @@ std::vector create_test_cases() { // -------------------------------------------------- tests.push_back({ "Near Degenerate - Close to Corner", - outer1, {hole1}, {0.5000001, 0.5000001}, false // just outside hole + outer1, + {hole1}, + {0.5000001, 0.5000001}, + false // just outside hole }); return tests; @@ -230,7 +326,8 @@ std::vector create_test_cases() { // ----------------------------- // Main Function // ----------------------------- -int main() { +int main() +{ std::cout << "Running Polygon Winding Number Tests...\n\n"; // Optional: Validate orientation logic @@ -242,15 +339,11 @@ int main() { auto test_cases = create_test_cases(); int passed = 0; - int total = test_cases.size(); + int total = test_cases.size(); for (const auto& tc : test_cases) { - if (!validate_test_case(tc)) { - continue; - } - if (run_test(tc)) { - passed++; - } + if (!validate_test_case(tc)) { continue; } + if (run_test(tc)) { passed++; } } std::cout << "\nSummary: " << passed << " / " << total << " tests passed.\n";