diff --git a/surface_integral/interface/SurfaceIntegrator.hpp b/surface_integral/interface/SurfaceIntegrator.hpp
index b9580c4..3f2fc14 100644
--- a/surface_integral/interface/SurfaceIntegrator.hpp
+++ b/surface_integral/interface/SurfaceIntegrator.hpp
@@ -3,66 +3,72 @@
#include
#include
#include
-
#include
namespace internal
{
-// 每个活动子面应具有以下结构
-struct parametric_plane {
- stl_vector_mp chain_vertices{}; // 链顶点
- flat_index_group chains{}; // 链组
- stl_vector_mp singularity_vertices{}; // 奇异顶点，即链的交点
- stl_vector_mp polar_vertices{}; // 极性顶点，即两个连接顶点周围的最小/最大 x/y
- stl_vector_mp parallel_start_vertices{}; // 平行起始顶点，即边 {v, v + 1} 平行于 x/y 轴
-};
-
+/**
+ * @brief Numerical integrator for parametric surfaces with trimming curves.
+ *
+ * This class computes integrals (e.g., area, mass, etc.) over trimmed parametric surfaces
+ * by subdividing the parameter domain and applying Gaussian quadrature.
+ *
+ * The integrator does not own the input data; it holds const references to ensure zero-copy semantics.
+ * Users must ensure that the lifetime of input data exceeds that of the integrator.
+ */
class SI_API integrator_t
{
public:
- // 构造函数，接受对 subface 的引用
- explicit integrator_t(const subface& surface);
- [[deprecated("Use calculate_new() instead")]] integrator_t(const subface& surface,
- const stl_vector_mp& u_breaks,
- double umin,
- double umax,
- double vmin,
- double vmax);
- integrator_t(const subface& surface, const parametric_plane& uv_plane);
+ /**
+ * @note This constructor does not copy the data; it stores const references.
+ * The caller is responsible for ensuring the validity of the referenced data
+ * throughout the lifetime of this integrator.
+ */
+ integrator_t(const stl_vector_mp>& surfaces,
+ const flat_hash_map_mp& uv_planes)
+ : m_surfaces(surfaces), m_uv_planes(uv_planes)
+ {
+ }
+ /// Default destructor
+ ~integrator_t() = default;
- // 设置 u_breaks（裁剪/分割线）
- void set_ubreaks(const stl_vector_mp& u_breaks);
- // 计算面积主函数
template
- double calculate(Func&& func, int gauss_order = 3) const;
- double compute_volume(int gauss_order = 4) const;
+ double calculate(int gauss_order, Func&& func) const;
+
+
+ template
+ double calculate_one_subface(const subface& subface,
+ const parametric_plane_t& param_plane,
+ int gauss_order,
+ Func&& func) const;
private:
- const subface& m_surface; // 引用原始曲面
- stl_vector_mp m_u_breaks; // 分割线信息（可选）
- parametric_plane m_uv_plane;
- double Umin = 0.0; // 参数域范围
- double Umax = 1.0;
- double Vmin = 0.0;
- double Vmax = 1.0;
-
- // 私有辅助函数
- // 直线u=u_val与边界的交点
- std::vector find_vertical_intersections(double u_val) const;
-
- 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) const;
-};
+ /// Non-owning reference to the list of subfaces
+ const stl_vector_mp>& m_surfaces;
+
+ /// Non-owning reference to the map of parametric planes (ID -> parametric_plane_t)
+ const flat_hash_map_mp& m_uv_planes;
+ stl_vector_mp compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max);
+
+ stl_vector_mp find_v_intersections_at_u(const subface& subface,
+ const parametric_plane_t& param_plane,
+ double u_val) const;
+
+ bool is_point_inside_domain(const subface& subface, const parametric_plane_t& param_plane, double u, double v) const;
+
+ bool is_u_near_singularity(double u, double tol = 1e-6) const;
+
+
+ void sort_and_unique_with_tol(stl_vector_mp& vec, double epsilon = 1e-8) const;
+};
-SI_API double newton_method(const std::function& F,
- double v_initial,
- double tolerance = 1e-8,
- int max_iterations = 100);
+double newton_method(const std::function& F,
+ double v_initial,
+ double tolerance = 1e-8,
+ int max_iterations = 100);
} // namespace internal
\ No newline at end of file
diff --git a/surface_integral/src/SurfaceIntegrator.cpp b/surface_integral/src/SurfaceIntegrator.cpp
index 3ee133e..a044294 100644
--- a/surface_integral/src/SurfaceIntegrator.cpp
+++ b/surface_integral/src/SurfaceIntegrator.cpp
@@ -8,111 +8,53 @@
namespace internal
{
-// Constructor 1: Initialize only with a reference to the surface
-integrator_t::integrator_t(const subface& surface) : m_surface(surface) {}
-
-// Constructor 2: Initialize with surface and u-breaks (e.g., trimming curves)
-integrator_t::integrator_t(const subface& surface,
- const stl_vector_mp& u_breaks,
- double umin,
- double umax,
- double vmin,
- double vmax)
- : m_surface(surface), m_u_breaks(u_breaks), Umin(umin), Umax(umax), Vmin(vmin), Vmax(vmax)
-{
-}
-
-integrator_t::integrator_t(const subface& surface, const parametric_plane& uv_plane)
- : m_surface(surface), m_uv_plane(uv_plane), Umin(0.0), Umax(0.0), Vmin(0.0), Vmax(0.0)
+// Main entry point to compute surface area
+// Main entry point to compute surface area
+template
+double integrator_t::calculate(int gauss_order, Func&& func) const
{
- if (!uv_plane.chain_vertices.empty()) {
- // 初始化为第一个点的坐标
- double min_u = uv_plane.chain_vertices[0].x();
- double max_u = uv_plane.chain_vertices[0].x();
- double min_v = uv_plane.chain_vertices[0].y();
- double max_v = uv_plane.chain_vertices[0].y();
-
- // 遍历所有链顶点
- for (const auto& pt : uv_plane.chain_vertices) {
- double u = pt.x();
- double v = pt.y();
-
- if (u < min_u) min_u = u;
- if (u > max_u) max_u = u;
- if (v < min_v) min_v = v;
- if (v > max_v) max_v = v;
- }
-
- Umin = min_u;
- Umax = max_u;
- Vmin = min_v;
- Vmax = max_v;
- } else {
- // 没有顶点时使用默认范围 [0, 1] × [0, 1]
- Umin = 0.0;
- Umax = 1.0;
- Vmin = 0.0;
- Vmax = 1.0;
+ total_integral = 0.0;
+ for (const auto& [subface_index, param_plane] : m_uv_planes) {
+ const auto& subface = m_subfaces[subface_index].object_ptr.get();
+ total_integral += calculate_one_subface(subface, param_plane, gauss_order, std::forward(func));
}
-
- std::set unique_vertex_indices;
-
- // 插入所有类型的顶点索引
- unique_vertex_indices.insert(uv_plane.singularity_vertices.begin(), uv_plane.singularity_vertices.end());
- unique_vertex_indices.insert(uv_plane.polar_vertices.begin(), uv_plane.polar_vertices.end());
- unique_vertex_indices.insert(uv_plane.parallel_start_vertices.begin(), uv_plane.parallel_start_vertices.end());
-
- std::set unique_u_values;
-
- for (uint32_t idx : unique_vertex_indices) {
- if (idx < uv_plane.chain_vertices.size()) {
- double u = uv_plane.chain_vertices[idx].x();
- unique_u_values.insert(u);
- }
- }
-
- // 转换为 vector 并排序（set 默认是有序的）
- m_u_breaks = stl_vector_mp(unique_u_values.begin(), unique_u_values.end());
+ return total_integral;
}
-// Set u-breaks (optional trimming or partitioning lines)
-void integrator_t::set_ubreaks(const stl_vector_mp& u_breaks) { m_u_breaks = u_breaks; }
-
-// Main entry point to compute surface area
template
-double integrator_t::calculate(Func&& func, int gauss_order) const
+double integrator_t::calculate_one_subface(const subface& subface, const parametric_plane_t& param_plane, int gauss_order, Func&& func) const
{
- auto solver = m_surface.fetch_solver_evaluator();
- // 在u方向进行高斯积分
+ auto solver = subface.fetch_solver_evaluator();
+ // Gaussian integration in u direction
auto u_integrand = [&](double u) {
- // 对每个u，找到v方向的精确交点
+ // Find exact v intersections for each u
std::vector v_breaks = find_vertical_intersections(u);
- // 在v方向进行高斯积分
+ // Gaussian integration in v direction
auto v_integrand = [&](double v) {
- // 判断点是否在有效域内
- if (is_point_inside_domain(u, v)) {
+ // Check if point is inside domain
+ if (is_point_inside_domain(subface, param_plane, u, v)) {
try {
- // 获取两个方向的 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
+ // Get evaluators for both directions
+ auto eval_du = subface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // Fix v, vary u → ∂/∂u
+ auto eval_dv = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // Fix u, vary v → ∂/∂v
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 p = res_u.f.template head<3>(); // Point (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||
+ // Area element: ||dU × dV||
Eigen::Vector3d cross = dU.cross(dV);
- double jacobian = cross.norm(); // 雅可比行列式（面积缩放因子）
+ double jacobian = cross.norm(); // Jacobian (area scaling factor)
return func(u, v, p, dU, dV) * jacobian;
} catch (...) {
- return 0.0; // 跳过奇异点
+ return 0.0; // Skip singular points
}
}
- return 0.0; // 点不在有效域内
+ return 0.0; // Point not in domain
};
double v_integral = 0.0;
@@ -120,9 +62,9 @@ double integrator_t::calculate(Func&& func, int gauss_order) const
double a = v_breaks[i];
double b = v_breaks[i + 1];
- // 检查区间中点是否有效
+ // Check midpoint validity
double mid_v = (a + b) / 2.0;
- if (is_point_inside_domain(u, mid_v)) {
+ if (is_point_inside_domain(subface, param_plane, 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;
@@ -132,18 +74,17 @@ double integrator_t::calculate(Func&& func, int gauss_order) const
return v_integral;
};
- // 在u方向积分
+ // Integrate in u direction
+ const auto& u_breaks = compute_u_breaks(param_plane);
double integral = 0.0;
for (size_t i = 0; i < u_breaks.size() - 1; ++i) {
double a = u_breaks[i];
double b = u_breaks[i + 1];
- // 检查区间中点是否有效
double mid_u = (a + b) / 2.0;
auto v_intersections = find_vertical_intersections(mid_u, self.outerEdges, self.Vmin, self.Vmax);
- if (!v_intersections.empty()) { // 确保该u区间有有效区域
-
+ if (!v_intersections.empty()) {
integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u));
}
}
@@ -152,6 +93,7 @@ double integrator_t::calculate(Func&& func, int gauss_order) const
}
// 在 integrator_t 类中添加：
+/*
double integrator_t::compute_volume(int gauss_order) const
{
double total_volume = 0.0;
@@ -163,7 +105,7 @@ double integrator_t::compute_volume(int gauss_order) const
// 内层：对 v 积分
auto v_integrand = [&](double v) -> double {
- if (!is_point_inside_domain(u, v)) {
+ if (!is_point_inside_domain(subface, u, v)) {
return 0.0;
}
@@ -193,7 +135,7 @@ double integrator_t::compute_volume(int gauss_order) const
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)) {
+ if (is_point_inside_domain(subface, u, mid_v)) {
v_integral += integrate_1D(a, b, v_integrand, gauss_order);
}
}
@@ -214,56 +156,126 @@ double integrator_t::compute_volume(int gauss_order) const
// 乘以 1/3
return std::abs(total_volume) / 3.0;
+}*/
+
+/**
+ * @brief Compute the u-parameter breakpoints for integration.
+ * @note The function currently uses std::set for uniqueness and sorting,
+ * but floating-point precision may cause near-duplicate values to be
+ * treated as distinct. A tolerance-based comparison is recommended.
+ * TODO: Use a tolerance-based approach to avoid floating-point precision issues
+ * when inserting u-values (e.g., merge values within 1e-12).
+ */
+stl_vector_mp compute_u_breaks(
+ const parametric_plane_t& param_plane,
+ double u_min,
+ double u_max)
+{
+ std::set break_set;
+
+ // Insert domain boundaries
+ break_set.insert(u_min);
+ break_set.insert(u_max);
+
+ // Insert u-values from special vertices (e.g., trimming curve vertices)
+ for (size_t i = 0; i < param_plane.chain_vertices.size(); ++i) {
+ const auto& vertices = param_plane.chain_vertices[i];
+ auto& vertex_flags = param_plane.vertex_special_flags[i];
+ for (size_t j = 0; j < vertices.size(); ++j) {
+ if (vertex_flags[j]) { break_set.insert( vertices[j].x());} // Special vertex u
+ }
+ }
+ // Return as vector (sorted and unique due to set)
+ return stl_vector_mp(break_set.begin(), break_set.end());
}
-// 直线u=u_val与边界的交点
-std::vector integrator_t::find_vertical_intersections(double u_val) const
+
+stl_vector_mp integrator_t::find_v_intersections_at_u(
+ const subface& subface,
+ const parametric_plane_t& param_plane,
+ double u_val) const
{
- std::vector intersections;
- std::vector uPositionFlags;
- uPositionFlags.reserve(m_uv_plane.chain_vertices.size());
-
- std::transform(m_uv_plane.chain_vertices.begin(),
- m_uv_plane.chain_vertices.end(),
- std::back_inserter(uPositionFlags),
- [&](const auto& currentVertex) {
- double uDifference = currentVertex.x() - u_val;
- if (uDifference < 0) return -1; // 在参考值左侧
- if (uDifference > 0) return 1; // 在参考值右侧
- return 0; // 等于参考值
- });
-
- uint32_t group_idx = 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++;
- if (element_idx + 1 < m_uv_plane.chains.start_indices[group_idx + 1]) {
- uint32_t vertex_idx1 = m_uv_plane.chains.index_group[element_idx];
- uint32_t vertex_idx2 = m_uv_plane.chains.index_group[element_idx + 1];
- if (uPositionFlags[vertex_idx1] * uPositionFlags[vertex_idx2] <= 0) {
- auto v1 = m_uv_plane.chain_vertices[vertex_idx1], v2 = m_uv_plane.chain_vertices[vertex_idx2];
- 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(parameter_u_t{}, u_val);
- auto solver_evaluator = m_surface.fetch_solver_evaluator();
- auto target_function = [&](double v) -> internal::implicit_equation_intermediate {
- 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);
- };
+ stl_vector_mp intersections;
+
+ // Iterate over each boundary chain
+ for (const auto& chain : param_plane.chain_vertices) {
+ const size_t n_vertices = chain.size();
+ if (n_vertices < 2) continue; // Skip degenerate chains
+
+ // Iterate over each edge in the chain (including closing edge if desired)
+ for (size_t i = 0; i < n_vertices; ++i) {
+ size_t j = (i + 1) % n_vertices; // Next vertex index (wraps around for closed chain)
+
+ const Eigen::Vector2d& v1 = chain[i]; // Current vertex: (u1, v1)
+ const Eigen::Vector2d& v2 = chain[j]; // Next vertex: (u2, v2)
+
+ double u1 = v1.x(), v1_val = v1.y();
+ double u2 = v2.x(), v2_val = v2.y();
+
+ // Classify position relative to u = u_val: -1 (left), 0 (on), +1 (right)
+ const double eps = 1e-12;
+ auto sign_cmp = [u_val, eps](double u) -> int {
+ if (u < u_val - eps) return -1;
+ if (u > u_val + eps) return 1;
+ return 0;
+ };
+
+ // Then use it
+ int pos1 = sign_cmp(u1);
+ int pos2 = sign_cmp(u2);
+
+ // Case 1: Both endpoints on the same side (and not on the line) → no intersection
+ if (pos1 * pos2 > 0) {
+ continue;
+ }
+
+ // Case 2: Both endpoints lie exactly on u = u_val → add both v-values
+ if (pos1 == 0 && pos2 == 0) {
+ intersections.push_back(v1_val);
+ intersections.push_back(v2_val);
+ }
+ // Case 3: One endpoint on u = u_val or segment crosses u = u_val
+ else if (std::abs(u1 - u2) < eps) {
+ // Vertical segment: if aligned with u_val, treat as overlapping
+ if (pos1 == 0 || pos2 == 0) {
+ intersections.push_back(v1_val);
+ intersections.push_back(v2_val);
+ }
+ }
+ else {
+ // General case: non-vertical segment crossing u = u_val
+ // Compute linear interpolation parameter t
+ double t = (u_val - u1) / (u2 - u1);
+ double v_initial = v1_val + t * (v2_val - v1_val); // Initial guess for v
+
+ // Fetch evaluators from subface for constraint solving
+ auto curve_evaluator = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u_val);
+ auto solver_evaluator = subface.fetch_solver_evaluator();
+
+ // Define target function: v ↦ residual of implicit equation
+ auto target_function = [&](double v) -> internal::implicit_equation_intermediate {
+ constraint_curve_intermediate temp_res = curve_evaluator(v);
+ auto full_res = solver_evaluator(std::move(temp_res));
+ return std::get(full_res);
+ };
+
+ // Refine solution using Newton-Raphson method
+ try {
double v_solution = newton_method(target_function, v_initial);
intersections.push_back(v_solution);
- } else {
- intersections.push_back(v1.y());
- intersections.push_back(v2.y());
+ } catch (...) {
+ // If Newton's method fails (e.g., divergence), fall back to linear approximation
+ intersections.push_back(v_initial);
}
}
}
}
- // 去重排序
- sort_and_unique_with_tol(intersections);
+ // Final step: sort and remove duplicates within tolerance
+ if (!intersections.empty()) {
+ sort_and_unique_with_tol(intersections, 1e-8);
+ }
+
return intersections;
}
@@ -274,65 +286,35 @@ std::vector integrator_t::find_vertical_intersections(double u_val) cons
NOTE: when v_intersect - v < threshold, further checks are required to accurately determine if the intersection point lies
precisely on the boundary segment.
*/
-bool integrator_t::is_point_inside_domain(double u, double v) const
+bool integrator_t::is_point_inside_domain(
+ const subface& subface,
+ const parametric_plane_t& param_plane,
+ 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++;
- if (element_idx + 1 < m_uv_plane.chains.start_indices[group_idx + 1]) {
- uint32_t vertex_idx1 = m_uv_plane.chains.index_group[element_idx];
- 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_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_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) -> internal::implicit_equation_intermediate {
- 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;
- }
-
- }
- }
- /*
- case v1.x() == v2.x() == u, do not count as intersection.but will cout in next iteration
- */
+ auto intersections = find_v_intersections_at_u(subface, param_plane, u);
+ const double tol_near = 1e-8;
+ const double tol_above = 1e-12;
+
+ uint32_t count = 0;
+ for (double v_int : intersections) {
+ double diff = v_int - v;
+ if (diff > tol_above) {
+ count++;
+ }
+ else if (std::abs(diff) < tol_near) {
+ return true; // on boundary → inside
}
}
- return intersection_count % 2 == 1; // in domain
+ return (count % 2) == 1;
}
bool integrator_t::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();
- if (std::abs(u - singular_u) < tol) { return true; }
- }
- // 可扩展：判断是否靠近极点、极性顶点等
- for (auto idx : m_uv_plane.polar_vertices) {
- double polar_u = m_uv_plane.chain_vertices[idx].x();
- if (std::abs(u - polar_u) < tol) { return true; }
- }
return false;
}
-void integrator_t::sort_and_unique_with_tol(std::vector& vec, double epsilon) const
+void integrator_t::sort_and_unique_with_tol(stl_vector_mp& vec, double epsilon) const
{
if (vec.empty()) return;