feat(surface-integral): add volume computation and improve domain checks in SurfaceAreaCalculator

primitive_process/interface/base/subface.hpp

@@ -67,6 +67,7 @@ EXTERN_C struct PE_API subface {
         double                  v) const = 0;
     virtual std::function<internal::equation_intermediate_t(internal::constraint_curve_intermediate &&)>
         fetch_solver_evaluator() const = 0;
+    virtual bool is_implicit_equation_intermediate() const {return true;};
 
     primitive_data_center_ptr_t         data_center{};
     internal::paired_model_matrix_ptr_t model_matrices{};

surface_integral/interface/SurfaceIntegrator.hpp

@@ -3,7 +3,7 @@
 #include <base/subface.hpp>
 #include <container/stl_alias.hpp>
 
-#include "flat_index_group.hpp"
+#include <container/wrapper/flat_index_group.hpp>
 
 namespace internal
 {
@@ -37,6 +37,7 @@ public:
     // 计算面积主函数
     template <typename Func>
     double calculate(Func&& func, int gauss_order = 3) const;
+    double compute_volume(int gauss_order = 4) const;
 
 private:
     const subface&        m_surface;  // 引用原始曲面
@@ -49,14 +50,15 @@ private:
 
     // 私有辅助函数
     // 直线u=u_val与边界的交点
-    std::vector<double> find_vertical_intersections(double u_val);
+    std::vector<double> find_vertical_intersections(double u_val) const;
 
-    bool   is_point_inside_domain(double u, double v);
+    bool   is_point_inside_domain(double u, double v) const;
-    bool   is_u_near_singularity(double u, double tol = 1e-6);
+    bool   is_u_near_singularity(double u, double tol = 1e-6) const;
 
-    void   sort_and_unique_with_tol(std::vector<double>& vec, double epsilon = 1e-8);
+    void   sort_and_unique_with_tol(std::vector<double>& vec, double epsilon = 1e-8) const;
 };
 
+
 double newton_method(const std::function<equation_intermediate_t(double)>& F,
                         double                                                v_initial,
                         double                                                tolerance      = 1e-8,

surface_integral/interface/flat_index_group.hpp

@@ -1,52 +0,0 @@
-#pragma once
-
-#include <container/stl_alias.hpp>
-#include <iterator/counting_iterator.hpp>
-
-struct flat_index_group {
-    stl_vector_mp<uint32_t> index_group{};   /// a list of indices in the group
-    stl_vector_mp<uint32_t> start_indices{}; /// a list of start indices for each group in the flat index group
-
-    inline size_t size() const { return start_indices.size() - 1; }
-
-    struct group_index_looper {
-        auto begin() const { return counting_iterator<uint32_t>{0}; }
-
-        auto end() const { return counting_iterator<uint32_t>{static_cast<uint32_t>(parent_group.size())}; }
-
-        const flat_index_group& parent_group{};
-    };
-
-    auto group_indices() const { return group_index_looper{*this}; }
-
-    struct group_looper {
-        auto begin() const
-        {
-            return std::next(parent_group.index_group.begin(), *(parent_group.start_indices.begin() + group_idx));
-        }
-
-        auto end() const
-        {
-            return std::next(parent_group.index_group.begin(), *(parent_group.start_indices.begin() + group_idx + 1));
-        }
-
-        const flat_index_group& parent_group{};
-        const uint32_t          group_idx{};
-    };
-
-    auto group(uint32_t group_idx) const { return group_looper{*this, group_idx}; }
-
-    inline auto group_index_iter_begin() const { return counting_iterator<uint32_t>{0}; }
-
-    inline auto group_index_iter_end() const { return counting_iterator<uint32_t>{static_cast<uint32_t>(size())}; }
-
-    inline auto group_begin(uint32_t group_idx) const
-    {
-        return std::next(index_group.begin(), *(start_indices.begin() + group_idx));
-    }
-
-    inline auto group_end(uint32_t group_idx) const
-    {
-        return std::next(index_group.begin(), *(start_indices.begin() + group_idx + 1));
-    }
-};

surface_integral/src/SurfaceIntegrator.cpp

@@ -3,6 +3,7 @@
 #include <Eigen/Geometry> // For vector and cross product operations
 #include <cmath>          // For math functions like sqrt
 #include <set>
+#include <iostream>
 
 namespace internal
 {
@@ -81,6 +82,7 @@ void SurfaceAreaCalculator::set_ubreaks(const stl_vector_mp<double>& u_breaks) {
 template <typename Func>
 double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
 {
+    auto solver = m_surface.fetch_solver_evaluator();
     // 在u方向进行高斯积分
     auto u_integrand = [&](double u) {
         // 对每个u，找到v方向的精确交点
@@ -89,13 +91,22 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
         // 在v方向进行高斯积分
         auto v_integrand = [&](double v) {
             // 判断点是否在有效域内
-            if (IsPointInsideself(u, v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) {
+            if (is_point_inside_domain(u, v)) {
                 try {
-                    gp_Pnt p;
+                    // 获取两个方向的 evaluator
-                    gp_Vec dU, dV;
+                    auto eval_du =  m_surface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // 固定 v，变 u → 得到 ∂/∂u
-                    surface->D1(u, v, p, dU, dV);
+                    auto eval_dv =  m_surface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // 固定 u，变 v → 得到 ∂/∂v
 
-                    const double jacobian = dU.Crossed(dV).Magnitude();
+                    auto res_u = eval_du(u); // f(u,v), grad_f = ∂r/∂u
+                    auto res_v = eval_dv(v); // f(u,v), grad_f = ∂r/∂v
+
+                    Eigen::Vector3d p  = res_u.f.template head<3>();      // 点坐标 (x,y,z)
+                    Eigen::Vector3d dU = res_u.grad_f.template head<3>(); // ∂r/∂u
+                    Eigen::Vector3d dV = res_v.grad_f.template head<3>(); // ∂r/∂v
+
+                    // ✅ 计算面积元：||dU × dV||
+                    Eigen::Vector3d cross = dU.cross(dV);
+                    double jacobian = cross.norm(); // 雅可比行列式（面积缩放因子）
                     return func(u, v, p, dU, dV) * jacobian;
                 } catch (...) {
                     return 0.0; // 跳过奇异点
@@ -111,7 +122,7 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
 
             // 检查区间中点是否有效
             double mid_v = (a + b) / 2.0;
-            if (IsPointInsideself(u, mid_v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) {
+            if (is_point_inside_domain(u, mid_v)) {
                 v_integral += integrate_1D(a, b, v_integrand, gauss_order);
             } else {
                 std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl;
@@ -140,8 +151,73 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const
     return integral;
 }
 
+// 在 SurfaceAreaCalculator 类中添加：
+double SurfaceAreaCalculator::compute_volume(int gauss_order) const
+{
+    double total_volume = 0.0;
+
+    // 外层：对 u 分段积分
+    auto u_integrand = [&](double u) {
+        std::vector<double> v_breaks = find_vertical_intersections(u);
+        double v_integral = 0.0;
+
+        // 内层：对 v 积分
+        auto v_integrand = [&](double v) -> double {
+            if (!is_point_inside_domain(u, v)) {
+                return 0.0;
+            }
+
+            try {
+                // 获取偏导数
+                auto eval_du = m_surface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // ∂/∂u
+                auto eval_dv = m_surface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // ∂/∂v
+
+                auto res_u = eval_du(u);
+                auto res_v = eval_dv(v);
+
+                Eigen::Vector3d p  = res_u.f.template head<3>();      // r(u,v)
+                Eigen::Vector3d dU = res_u.grad_f.template head<3>(); // r_u
+                Eigen::Vector3d dV = res_v.grad_f.template head<3>(); // r_v
+
+                // 计算 r · (r_u × r_v)
+                Eigen::Vector3d cross = dU.cross(dV);
+                double mixed_product = p.dot(cross);
+
+                return mixed_product; // 注意：不是 norm，是点积
+            } catch (...) {
+                return 0.0;
+            }
+        };
+
+        for (size_t i = 0; i < v_breaks.size() - 1; ++i) {
+            double a = v_breaks[i];
+            double b = v_breaks[i + 1];
+            double mid_v = (a + b) / 2.0;
+            if (is_point_inside_domain(u, mid_v)) {
+                v_integral += integrate_1D(a, b, v_integrand, gauss_order);
+            }
+        }
+
+        return v_integral;
+    };
+
+    // 在 u 方向积分
+    for (size_t i = 0; i < m_u_breaks.size() - 1; ++i) {
+        double a = m_u_breaks[i];
+        double b = m_u_breaks[i + 1];
+        double mid_u = (a + b) / 2.0;
+        auto v_intersections = find_vertical_intersections(mid_u);
+        if (!v_intersections.empty()) {
+            total_volume += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u));
+        }
+    }
+
+    // 乘以 1/3
+    return std::abs(total_volume) / 3.0;
+}
+
 // 直线u=u_val与边界的交点
-std::vector<double> SurfaceAreaCalculator::find_vertical_intersections(double u_val)
+std::vector<double> SurfaceAreaCalculator::find_vertical_intersections(double u_val) const
 {
     std::vector<double>  intersections;
     std::vector<int32_t> uPositionFlags;
@@ -168,7 +244,7 @@ std::vector<double> SurfaceAreaCalculator::find_vertical_intersections(double u_
                 if (v1.x() != v2.x()) {
                     // "The line segment is vertical (u₁ == u₂), so there is no unique v value corresponding to the given u."
                     double v_initial        = v1.y() + (v2.y() - v1.y()) * (u_val - v1.x()) / (v2.x() - v1.x());
-                    auto   curve_evaluator  = m_surface.fetch_curve_constraint_evaluator(u_val);
+                    auto   curve_evaluator  = m_surface.fetch_curve_constraint_evaluator(parameter_u_t{}, u_val); 
                     auto   solver_evaluator = m_surface.fetch_solver_evaluator();
                     auto   target_function  = [&](double v) -> equation_intermediate_t {
                         constraint_curve_intermediate temp_res = curve_evaluator(v);
@@ -196,8 +272,9 @@ std::vector<double> SurfaceAreaCalculator::find_vertical_intersections(double u_
  NOTE:  when v_intersect - v < threshold, further checks are required to accurately determine if the intersection point lies
  precisely on the boundary segment.
 */
-bool is_point_inside_domain(double u, double v)
+bool SurfaceAreaCalculator::is_point_inside_domain(double u, double v) const
 {
+    bool is_implicit_equation_intermediate = m_surface.is_implicit_equation_intermediate();
     uint32_t group_idx = 0, intersection_count = 0;
     for (uint32_t element_idx = 0; element_idx < m_uv_plane.chains.index_group.size() - 1; element_idx++) {
         if (element_idx > m_uv_plane.chains.start_indices[group_idx + 1]) group_idx++;
@@ -206,19 +283,30 @@ bool is_point_inside_domain(double u, double v)
             uint32_t vertex_idx2 = m_uv_plane.chains.index_group[element_idx + 1];
             auto     v1 = m_uv_plane.chain_vertices[vertex_idx1], v2 = m_uv_plane.chain_vertices[vertex_idx2];
             if ((v1.x() <= u && v2.x() > u) || (v2.x() < u && v1.x() >= u)) {
-                double v_intersected = v1.y() + (v2.y() - v1.y()) * (u_val - v1.x()) / (v2.x() - v1.x());
+                double v_initial = v1.y() + (v2.y() - v1.y()) * (u - v1.x()) / (v2.x() - v1.x());
-                if (v_interdected - v >= 1e-6) {
+                if (v_initial - v >= 1e-6) {
                     intersection_count++;
-                } else if (std::abs(v_intersected - v) < 1e-6) {
+                } else if (std::abs(v_initial - v) < 1e-6) {
-                    auto curve_evaluator  = m_surface.fetch_curve_constraint_evaluator(u);
+                    // Only use Newton's method for implicit surfaces (scalar equation)
+                    //    Newton requires f(v) and df/dv as scalars — only implicit provides this.
+                    //    Skip parametric surfaces (vector residual) — treat initial guess as final
+                    if (is_implicit_equation_intermediate) {
+                        auto curve_evaluator  = m_surface.fetch_curve_constraint_evaluator(parameter_u_t{}, u);
                         auto solver_evaluator = m_surface.fetch_solver_evaluator();
                         auto target_function  = [&](double v) -> equation_intermediate_t {
                             constraint_curve_intermediate temp_res = curve_evaluator(v);
-                        return solver_evaluator(std::move(temp_res));
+                            auto full_res = solver_evaluator(std::move(temp_res));
+                            // ensure solver_eval returns implicit_equation_intermediate)
+                            return std::get<internal::implicit_equation_intermediate>(full_res);
                         };
                         double v_solution = newton_method(target_function, v_initial);
                         if (std::abs(v_solution - v) > 0) { intersection_count++; }
                     }
+                    else{
+                        continue;
+                    }
+                    
+                }
             }
             /*
             case v1.x() == v2.x() == u, do not count as intersection.but will cout in next iteration
@@ -228,7 +316,7 @@ bool is_point_inside_domain(double u, double v)
     return intersection_count % 2 == 1; // in domain
 }
 
-bool is_u_near_singularity(double u, double tol = 1e-6)
+bool SurfaceAreaCalculator::is_u_near_singularity(double u, double tol) const
 {
     for (auto idx : m_uv_plane.singularity_vertices) {
         double singular_u = m_uv_plane.chain_vertices[idx].x();
@@ -242,7 +330,7 @@ bool is_u_near_singularity(double u, double tol = 1e-6)
     return false;
 }
 
-void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector<double>& vec, double epsilon)
+void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector<double>& vec, double epsilon) const
 {
     if (vec.empty()) return;
 
@@ -259,35 +347,34 @@ void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector<double>& vec, d
     vec.resize(write_index + 1);
 }
 
-// 牛顿法求解器
+// Only accepts functions that return implicit_equation_intermediate
-double newton_method(const std::function<equation_intermediate_t(double)>& F,
+double newton_method(
+    const std::function<internal::implicit_equation_intermediate(double)>& F,
     double v_initial,
-                     double                                                tolerance,
+    double tolerance = 1e-6,
-                     int                                                   max_iterations)
+    int max_iterations = 20)
 {
     double v = v_initial;
 
     for (int i = 0; i < max_iterations; ++i) {
-        equation_intermediate_t res = F(v);
+        auto res = F(v);           //  Known type: implicit_equation_intermediate
-        double                  f   = res.f;
+        double f_val = res.f;
-        double                  df  = res.df;
+        double df_val = res.df;
 
-        if (std::abs(f) < tolerance) {
+        if (std::abs(f_val) < tolerance) {
-            std::cout << "✅ Converged in " << i + 1 << " iterations. v = " << v << std::endl;
+            std::cout << "Converged at v = " << v << std::endl;
             return v;
         }
 
-        if (std::abs(df) < 1e-10) {
+        if (std::abs(df_val) < 1e-10) {
-            std::cerr << "⚠️ Derivative near zero. No convergence." << std::endl;
+            std::cerr << "Derivative near zero." << std::endl;
             return v;
         }
 
-        v -= f / df;
+        v = v - f_val / df_val;
-
-        std::cout << "Iteration " << i + 1 << ": v = " << v << ", f = " << f << std::endl;
     }
 
-    std::cerr << "❌ Did not converge within " << max_iterations << " iterations." << std::endl;
+    std::cerr << "Newton failed to converge." << std::endl;
     return v;
 }