diff --git a/surface_integral/include/quadrature.hpp b/surface_integral/include/quadrature.hpp index 2a4a813..243f815 100644 --- a/surface_integral/include/quadrature.hpp +++ b/surface_integral/include/quadrature.hpp @@ -25,4 +25,13 @@ double tanh_sinh_integrate_1D(double a, double b, Func&& func, int q) { result += w * func(x); } return result; +} + +template +double integrate_1D(double a, double b, Func&& func, int q, bool use_tanh_sinh = false) { + if (use_tanh_sinh) { + return tanh_sinh_integrate_1D(a, b, std::forward(func), q); + } else { + return gauss_integrate_1D(a, b, std::forward(func), q); + } } \ No newline at end of file diff --git a/surface_integral/interface/SurfaceIntegrator.hpp b/surface_integral/interface/SurfaceIntegrator.hpp index 5c5587c..8a59bed 100644 --- a/surface_integral/interface/SurfaceIntegrator.hpp +++ b/surface_integral/interface/SurfaceIntegrator.hpp @@ -52,6 +52,7 @@ private: std::vector find_vertical_intersections(double u_val); bool is_point_inside_domain(double u, double v); + bool is_u_near_singularity(double u, double tol = 1e-6); void sort_and_unique_with_tol(std::vector& vec, double epsilon = 1e-8); }; diff --git a/surface_integral/src/SurfaceIntegrator.cpp b/surface_integral/src/SurfaceIntegrator.cpp index ac35e8f..9d9b002 100644 --- a/surface_integral/src/SurfaceIntegrator.cpp +++ b/surface_integral/src/SurfaceIntegrator.cpp @@ -1,7 +1,7 @@ #include "SurfaceIntegrator.hpp" #include "quadrature.hpp" -#include // For vector and cross product operations -#include // For math functions like sqrt +#include // For vector and cross product operations +#include // For math functions like sqrt #include namespace internal @@ -11,12 +11,12 @@ namespace internal SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface) : m_surface(surface) {} // Constructor 2: Initialize with surface and u-breaks (e.g., trimming curves) -SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface, +SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface, const stl_vector_mp& u_breaks, - double umin, - double umax, - double vmin, - double vmax) + double umin, + double umax, + double vmin, + double vmax) : m_surface(surface), m_u_breaks(u_breaks), Umin(umin), Umax(umax), Vmin(vmin), Vmax(vmax) { } @@ -112,7 +112,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)) { - v_integral += gauss_integrate_1D(a, b, v_integrand, gauss_order); + v_integral += integrate_1D(a, b, v_integrand, gauss_order); } else { std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl; } @@ -133,9 +133,7 @@ double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const if (!v_intersections.empty()) { // 确保该u区间有有效区域 - double integralp = gauss_integrate_1D(a, b, u_integrand, gauss_order); - integral += integralp; - std::cout << "integral " << i << ": " << integralp << std::endl; + integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); } } @@ -208,30 +206,41 @@ 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) { 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 { + 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++; - } + 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 } +bool is_u_near_singularity(double u, double tol = 1e-6) +{ + 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 SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector& vec, double epsilon) { @@ -252,9 +261,9 @@ void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector& vec, d // 牛顿法求解器 double newton_method(const std::function& F, - double v_initial, - double tolerance, - int max_iterations) + double v_initial, + double tolerance, + int max_iterations) { double v = v_initial;