diff --git a/primitive_process/interface/base/subface.hpp b/primitive_process/interface/base/subface.hpp index 5b4228f..e08c731 100644 --- a/primitive_process/interface/base/subface.hpp +++ b/primitive_process/interface/base/subface.hpp @@ -67,6 +67,7 @@ EXTERN_C struct PE_API subface { double v) const = 0; virtual std::function 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{}; diff --git a/surface_integral/interface/SurfaceIntegrator.hpp b/surface_integral/interface/SurfaceIntegrator.hpp index 8a59bed..61de135 100644 --- a/surface_integral/interface/SurfaceIntegrator.hpp +++ b/surface_integral/interface/SurfaceIntegrator.hpp @@ -3,7 +3,7 @@ #include #include -#include "flat_index_group.hpp" +#include namespace internal { @@ -37,6 +37,7 @@ public: // 计算面积主函数 template 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 find_vertical_intersections(double u_val); + std::vector find_vertical_intersections(double u_val) const; - bool is_point_inside_domain(double u, double v); - bool is_u_near_singularity(double u, double tol = 1e-6); + bool is_point_inside_domain(double u, double v) const; + bool is_u_near_singularity(double u, double tol = 1e-6) const; - void sort_and_unique_with_tol(std::vector& vec, double epsilon = 1e-8); + void sort_and_unique_with_tol(std::vector& vec, double epsilon = 1e-8) const; }; + double newton_method(const std::function& F, double v_initial, double tolerance = 1e-8, diff --git a/surface_integral/interface/flat_index_group.hpp b/surface_integral/interface/flat_index_group.hpp deleted file mode 100644 index ca1bb91..0000000 --- a/surface_integral/interface/flat_index_group.hpp +++ /dev/null @@ -1,52 +0,0 @@ -#pragma once - -#include -#include - -struct flat_index_group { - stl_vector_mp index_group{}; /// a list of indices in the group - stl_vector_mp 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{0}; } - - auto end() const { return counting_iterator{static_cast(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{0}; } - - inline auto group_index_iter_end() const { return counting_iterator{static_cast(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)); - } -}; \ No newline at end of file diff --git a/surface_integral/src/SurfaceIntegrator.cpp b/surface_integral/src/SurfaceIntegrator.cpp index 9d9b002..3cb85f0 100644 --- a/surface_integral/src/SurfaceIntegrator.cpp +++ b/surface_integral/src/SurfaceIntegrator.cpp @@ -3,6 +3,7 @@ #include // For vector and cross product operations #include // For math functions like sqrt #include +#include namespace internal { @@ -81,6 +82,7 @@ void SurfaceAreaCalculator::set_ubreaks(const stl_vector_mp& u_breaks) { template 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; - gp_Vec dU, dV; - surface->D1(u, v, p, dU, dV); + // 获取两个方向的 evaluator + auto eval_du = m_surface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // 固定 v，变 u → 得到 ∂/∂u + 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 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 SurfaceAreaCalculator::find_vertical_intersections(double u_val) +std::vector SurfaceAreaCalculator::find_vertical_intersections(double u_val) const { std::vector intersections; std::vector uPositionFlags; @@ -168,7 +244,7 @@ std::vector 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 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,18 +283,29 @@ 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()); - if (v_interdected - v >= 1e-6) { + double v_initial = v1.y() + (v2.y() - v1.y()) * (u - v1.x()) / (v2.x() - v1.x()); + if (v_initial - v >= 1e-6) { intersection_count++; - } else if (std::abs(v_intersected - v) < 1e-6) { - auto curve_evaluator = m_surface.fetch_curve_constraint_evaluator(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)); - }; - double v_solution = newton_method(target_function, v_initial); - if (std::abs(v_solution - v) > 0) { intersection_count++; } + } else if (std::abs(v_initial - v) < 1e-6) { + // 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); + auto full_res = solver_evaluator(std::move(temp_res)); + // ensure solver_eval returns implicit_equation_intermediate) + return std::get(full_res); + }; + double v_solution = newton_method(target_function, v_initial); + if (std::abs(v_solution - v) > 0) { intersection_count++; } + } + else{ + continue; + } + } } /* @@ -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& vec, double epsilon) +void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector& vec, double epsilon) const { if (vec.empty()) return; @@ -259,35 +347,34 @@ void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector& vec, d vec.resize(write_index + 1); } -// 牛顿法求解器 -double newton_method(const std::function& F, - double v_initial, - double tolerance, - int max_iterations) +// Only accepts functions that return implicit_equation_intermediate +double newton_method( + const std::function& F, + double v_initial, + double tolerance = 1e-6, + int max_iterations = 20) { double v = v_initial; for (int i = 0; i < max_iterations; ++i) { - equation_intermediate_t res = F(v); - double f = res.f; - double df = res.df; + auto res = F(v); // Known type: implicit_equation_intermediate + double f_val = res.f; + double df_val = res.df; - if (std::abs(f) < tolerance) { - std::cout << "✅ Converged in " << i + 1 << " iterations. v = " << v << std::endl; + if (std::abs(f_val) < tolerance) { + std::cout << "Converged at v = " << v << std::endl; return v; } - if (std::abs(df) < 1e-10) { - std::cerr << "⚠️ Derivative near zero. No convergence." << std::endl; + if (std::abs(df_val) < 1e-10) { + std::cerr << "Derivative near zero." << std::endl; return v; } - v -= f / df; - - std::cout << "Iteration " << i + 1 << ": v = " << v << ", f = " << f << std::endl; + v = v - f_val / df_val; } - std::cerr << "❌ Did not converge within " << max_iterations << " iterations." << std::endl; + std::cerr << "Newton failed to converge." << std::endl; return v; }