From d1b82fdcd7615ac8fcaafe1f7db5aa1fb96f6701 Mon Sep 17 00:00:00 2001 From: mckay Date: Fri, 18 Jul 2025 18:03:32 +0800 Subject: [PATCH] Refactor SurfaceAreaCalculator: streamline constructors, improve function naming, and enhance integration methods --- .../interface/SurfaceIntegrator.hpp | 52 ++-- surface_integral/src/SurfaceIntegrator.cpp | 253 +++++++++--------- 2 files changed, 147 insertions(+), 158 deletions(-) diff --git a/surface_integral/interface/SurfaceIntegrator.hpp b/surface_integral/interface/SurfaceIntegrator.hpp index 40f036c..b30f1c9 100644 --- a/surface_integral/interface/SurfaceIntegrator.hpp +++ b/surface_integral/interface/SurfaceIntegrator.hpp @@ -1,7 +1,7 @@ #pragma once #include - +#include "surface_integral/include/quadrature.hpp" namespace internal { @@ -19,13 +19,12 @@ class SurfaceAreaCalculator public: // 构造函数，接受对 subface 的引用 explicit SurfaceAreaCalculator(const subface& surface); - [[deprecated("Use calculate_new() instead")]] - SurfaceAreaCalculator(const subface& surface, - const stl_vector_map& u_breaks, - double umin, - double umax, - double vmin, - double vmax); + [[deprecated("Use calculate_new() instead")]] SurfaceAreaCalculator(const subface& surface, + const stl_vector_map& u_breaks, + double umin, + double umax, + double vmin, + double vmax); SurfaceAreaCalculator(const subface& surface, const parametric_plane& uv_plane); @@ -39,31 +38,32 @@ public: private: const subface& m_surface; // 引用原始曲面 stl_vector_map m_u_breaks; // 分割线信息（可选） - parametric_plane m_uv_plane; + parametric_plane m_uv_plane; double Umin = 0.0; // 参数域范围 double Umax = 1.0; double Vmin = 0.0; double Vmax = 1.0; - - // 私有辅助函数（可扩展） - template - double SurfaceAreaCalculator::GaussIntegrate1D(double a, double b, Func&& func, int q) const; + // 私有辅助函数 // 直线u=u_val与边界的交点 - std::vector FindVerticalIntersectionsOCCT(double u_val, - const std::vector& edges, - double v_min, - double v_max); + std::vector find_vertical_intersections(double u_val, + const std::vector& edges, + double v_min, + double v_max); - bool IsPointInsideDomain(double u, - double v, - const std::vector& outerEdges, - const std::vector& innerEdges, - double u_min, - double u_max, - double v_min, - double v_max); - double integrate_over_uv(int num_samples) const; + bool is_point_inside_domain(double u, + double v, + const std::vector& outerEdges, + const std::vector& innerEdges, + double u_min, + double u_max, + double v_min, + double v_max); + double newton_method(const std::function& F, + double v_initial, + double tolerance = 1e-8, + int max_iterations = 100); + void sort_and_unique_with_tol(std::vector& vec, double epsilon = 1e-8); }; } // namespace internal \ No newline at end of file diff --git a/surface_integral/src/SurfaceIntegrator.cpp b/surface_integral/src/SurfaceIntegrator.cpp index f96eb62..4917f25 100644 --- a/surface_integral/src/SurfaceIntegrator.cpp +++ b/surface_integral/src/SurfaceIntegrator.cpp @@ -20,9 +20,8 @@ SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surfa } SurfaceAreaCalculator::SurfaceAreaCalculator(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) { - + : m_surface(surface), m_uv_plane(uv_plane), Umin(0.0), Umax(0.0), Vmin(0.0), Vmax(0.0) +{ if (!uv_plane.chain_vertices.empty()) { // 初始化为第一个点的坐标 double min_u = uv_plane.chain_vertices[0].x(); @@ -83,28 +82,12 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const // 在u方向进行高斯积分 auto u_integrand = [&](double u) { // 对每个u，找到v方向的精确交点 - std::vector v_breaks = {self.Vmin, self.Vmax}; - - auto outer_intersections = FindVerticalIntersectionsOCCT(u, self.outerEdges, self.Vmin, self.Vmax); - v_breaks.insert(v_breaks.end(), outer_intersections.begin(), outer_intersections.end()); - - auto inner_intersections = FindVerticalIntersectionsOCCT(u, self.innerEdges, self.Vmin, self.Vmax); - v_breaks.insert(v_breaks.end(), inner_intersections.begin(), inner_intersections.end()); - - // 排序并去重 - sort_and_unique_with_tol(v_breaks, 1e-6); + std::vector v_breaks = find_vertical_intersections(u); // 在v方向进行高斯积分 auto v_integrand = [&](double v) { // 判断点是否在有效域内 - if (IsPointInsideself(u, - v, - self.outerEdges, - self.innerEdges, - self.Umin, - self.Umax, - self.Vmin, - self.Vmax)) { + if (IsPointInsideself(u, v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) { try { gp_Pnt p; gp_Vec dU, dV; @@ -126,15 +109,8 @@ 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)) { - v_integral += GaussIntegrate1D(a, b, v_integrand, gauss_order); + if (IsPointInsideself(u, mid_v, self.outerEdges, self.innerEdges, self.Umin, self.Umax, self.Vmin, self.Vmax)) { + v_integral += gauss_integrate_1D(a, b, v_integrand, gauss_order); } else { std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl; } @@ -151,11 +127,11 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const // 检查区间中点是否有效 double mid_u = (a + b) / 2.0; - auto v_intersections = FindVerticalIntersectionsOCCT(mid_u, self.outerEdges, self.Vmin, self.Vmax); + auto v_intersections = find_vertical_intersections(mid_u, self.outerEdges, self.Vmin, self.Vmax); if (!v_intersections.empty()) { // 确保该u区间有有效区域 - double integralp = GaussIntegrate1D(a, b, u_integrand, gauss_order); + double integralp = gauss_integrate_1D(a, b, u_integrand, gauss_order); integral += integralp; std::cout << "integral " << i << ": " << integralp << std::endl; } @@ -164,131 +140,144 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const return integral; } -// 在区间[a,b]上进行高斯积分采样 -template -double SurfaceAreaCalculator::GaussIntegrate1D(double a, double b, Func&& func, int q) const { - - double sum = 0.0; - for (int i = 0; i < q; ++i) { - double x=a+(b-a)*GaussQuad::x(q, i); - double w=GaussQuad::w(q, i); - sum += w * func(x); - } - return sum * (b-a); -} - // 直线u=u_val与边界的交点 -std::vector SurfaceAreaCalculator::FindVerticalIntersectionsOCCT( - double u_val, - const std::vector& edges, - double v_min, - double v_max) +std::vector SurfaceAreaCalculator::find_vertical_intersections(double u_val) { - std::vector intersections; - - // 创建垂直线(u=u_val) - gp_Lin2d vertical_line(gp_Pnt2d(u_val, v_min), gp_Dir2d(0, 1)); - Handle(Geom2d_Line) hVerticalLine = new Geom2d_Line(vertical_line); - - for (const auto& edge : edges) { - if (edge.IsNull()) { - continue; - } - - // 使用Geom2dAPI_InterCurveCurve求交 - Geom2dAPI_InterCurveCurve intersector(edge, hVerticalLine); - - if (intersector.NbPoints() > 0) { - for (int i = 1; i <= intersector.NbPoints(); ++i) { - gp_Pnt2d pt = intersector.Point(i); - double v = pt.Y(); - if (v >= v_min && v <= v_max) { - intersections.push_back(v); + 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(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); + intersections.push_back(v_solution); + } else { + intersections.push_back(v1.y()); + intersections.push_back(v2.y()); } } } - - // 补充采样法确保可靠性 - Geom2dAdaptor_Curve adaptor(edge); - double first = adaptor.FirstParameter(); - double last = adaptor.LastParameter(); - - const int samples = 20; - gp_Pnt2d prev_p; - adaptor.D0(first, prev_p); - double prev_v = prev_p.Y(); - - for (int i = 1; i <= samples; ++i) { - double param = first + i * (last - first) / samples; - gp_Pnt2d curr_p; - adaptor.D0(param, curr_p); - double curr_v = curr_p.Y(); - - if ((prev_v - u_val) * (curr_v - u_val) <= 0) { - // 牛顿迭代法 - double t = param - (last - first)/samples; - for (int iter = 0; iter < 5; ++iter) { - gp_Pnt2d p; - gp_Vec2d deriv; - adaptor.D1(t, p, deriv); - if (std::abs(p.X() - u_val) < Precision::Confusion()) { - if (p.Y() >= v_min && p.Y() <= v_max) { - intersections.push_back(p.Y()); - } - break; - } - t -= (p.X() - u_val) / deriv.X(); - } - } - prev_v = curr_v; - } } // 去重排序 - std::sort(intersections.begin(), intersections.end()); - intersections.erase(std::unique(intersections.begin(), intersections.end()), intersections.end()); + sort_and_unique_with_tol(intersections); return intersections; } +/* + point (u, v) is inside the domain by ray-casting algorithm + To determine whether a point (u, v) is inside or outside a domain by counting the intersections of a vertical ray starting + from the point and extending upwards, + 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) +{ + 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_intersected = v1.y() + (v2.y() - v1.y()) * (u_val - v1.x()) / (v2.x() - v1.x()); + if (v_interdected - 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++; + } + } + } + /* + case v1.x() == v2.x() == u, do not count as intersection.but will cout in next iteration + */ + } + } + return intersection_count % 2 == 1; // in domain +} -// Private method: Numerical integration over the UV parameter domain -double SurfaceAreaCalculator::integrate_over_uv(int num_samples) const +// 牛顿法求解器 +double SurfaceAreaCalculator::newton_method(const std::function& F, + double v_initial, + double tolerance, + int max_iterations) { - double area = 0.0; + double v = v_initial; - double du = 1.0 / num_samples; - double dv = 1.0 / num_samples; + for (int i = 0; i < max_iterations; ++i) { + equation_intermediate_t res = F(v); + double f = res.f; + double df = res.df; - // Fetch evaluators for point and derivative calculation - auto point_func = m_surface.fetch_point_by_param_evaluator(); - auto solver_func = m_surface.fetch_solver_evaluator(); + if (std::abs(f) < tolerance) { + std::cout << "✅ Converged in " << i + 1 << " iterations. v = " << v << std::endl; + return v; + } - for (int i = 0; i < num_samples; ++i) { - double u = (i + 0.5) * du; // Midpoint sampling in u-direction + if (std::abs(df) < 1e-10) { + std::cerr << "⚠️ Derivative near zero. No convergence." << std::endl; + return v; + } - for (int j = 0; j < num_samples; ++j) { - double v = (j + 0.5) * dv; // Midpoint sampling in v-direction + v -= f / df; - // Evaluate the surface point at (u, v) - Eigen::Vector4d pt = point_func(Eigen::Vector2d(u, v)); + std::cout << "Iteration " << i + 1 << ": v = " << v << ", f = " << f << std::endl; + } - // Get constraint and solve intermediate result - auto constraint = m_surface.fetch_curve_constraint_evaluator(internal::parameter_v_t{}, v)(u); - auto result = std::get(solver_func(std::move(constraint))); + std::cerr << "❌ Did not converge within " << max_iterations << " iterations." << std::endl; + return v; +} - // Extract partial derivatives ∂S/∂u and ∂S/∂v - Eigen::Vector3d dSdu = result.grad_f.col(0); - Eigen::Vector3d dSdv = result.grad_f.col(1); +void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector& vec, double epsilon) +{ + if (vec.empty()) return; - // Compute the differential area element (cross product norm) - double dA = dSdu.cross(dSdv).norm(); + std::sort(vec.begin(), vec.end()); - // Accumulate area using numerical integration - area += dA * du * dv; + size_t write_index = 0; + for (size_t read_index = 1; read_index < vec.size(); ++read_index) { + if (std::fabs(vec[read_index] - vec[write_index]) > epsilon) { + ++write_index; + vec[write_index] = vec[read_index]; } } - return area; + vec.resize(write_index + 1); } + } // namespace internal \ No newline at end of file