style: style: format subface_integrator header and source files

surface_integral/interface/mesh_algorithm.hpp

@@ -1,104 +1,96 @@
 #include <cmath>
 #include <cstdio>
-#include <solve.h>  // polymesh_t
-
+#include <solve.h> // 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
 }

surface_integral/interface/subface_integrator.hpp

@@ -62,13 +62,12 @@ private:
     stl_vector_mp<double> 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<double>&     intersections,
-                                                    stl_vector_mp<uint16_t>&   intersected_chains
-                                                ) const;
+                                   const parametric_plane_t& param_plane,
+                                   double                    u_val,
+                                   double                    v_min,
+                                   double                    v_max,
+                                   stl_vector_mp<double>&    intersections,
+                                   stl_vector_mp<uint16_t>&  intersected_chains) const;
 
     bool is_point_inside_domain(const subface& subface, const parametric_plane_t& param_plane, double u, double v) const;
 

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<PolygonRing>& holes)
+bool point_in_polygon_with_holes(const Point2d& P, const PolygonRing& outer_ring, const std::vector<PolygonRing>& 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<PolygonRing>& holes,
-    const std::string& context) 
+bool validate_polygon_input(const Point2d&                  query_point,
+                            const PolygonRing&              outer_ring,
+                            const std::vector<PolygonRing>& 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;
 }

surface_integral/src/subface_integrator.cpp

@@ -9,7 +9,7 @@ namespace internal
 {
 
 integrator_t(const stl_vector_mp<object_with_index_mapping<subface>>& surfaces,
-                const flat_hash_map_mp<uint32_t, parametric_plane_t>&    uv_planes)
+             const flat_hash_map_mp<uint32_t, parametric_plane_t>&    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<object_with_index_mapping<subface>>& 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<double> v_intersections;
+        stl_vector_mp<double>   v_intersections;
         stl_vector_mp<uint16_t> 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<double> v_intersections;
+        double                  mid_u = (a + b) / 2.0;
+        stl_vector_mp<double>   v_intersections;
         stl_vector_mp<uint16_t> 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<double> integrator_t::compute_u_breaks(const parametric_plane_t& p
     return stl_vector_mp<double>(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<Eigen::AlignedBox2d>& chain_bboxs,
-                                                              double                    u_val,
-                                                              double                    v_min,
-                                                              double                    v_max,
-                                                              stl_vector_mp<double>&     intersections,
-                                                              stl_vector_mp<uint16_t>&   intersected_chains
-                                                            ) const
+void integrator_t::find_v_intersections_at_u(const subface&                            subface,
+                                             const parametric_plane_t&                 param_plane,
+                                             const stl_vector_mp<Eigen::AlignedBox2d>& chain_bboxs,
+                                             double                                    u_val,
+                                             double                                    v_min,
+                                             double                                    v_max,
+                                             stl_vector_mp<double>&                    intersections,
+                                             stl_vector_mp<uint16_t>&                  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 || 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<stl_vector_mp<uint16_t>> & chain_group_indices,
-                           uint16_t chain_idx1,
-                           uint16_t chain_idx2) const
+bool is_edge_inside_domain(const stl_vector_mp<stl_vector_mp<uint16_t>>& 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<double>& vec, double epsilon) const

surface_integral/test/SurfaceIntegrator_test.cpp

@@ -5,13 +5,14 @@
 #include <SurfaceIntegrator.hpp>
 
 // ✅ 封装成函数：通过引用填充一个 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<object_with_index_mapping<subface>> subfaces;
-    flat_hash_map_mp<uint32_t, parametric_plane_t> uv_planes;
+    flat_hash_map_mp<uint32_t, parametric_plane_t>    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.}

surface_integral/test/test_mesh.cpp

@@ -1,37 +1,44 @@
 #include <iostream>
 #include <mesh_algorithm.hpp>
 
-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;
 }

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<PolygonRing> 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<TestCase> create_test_cases() {
+std::vector<TestCase> create_test_cases()
+{
     std::vector<TestCase> 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<TestCase> 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<TestCase> 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<TestCase> 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<TestCase> 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<TestCase> 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<TestCase> 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";