diff --git a/surface_integral/interface/subface_integrator.hpp b/surface_integral/interface/subface_integrator.hpp index b8688fe..96fadd7 100644 --- a/surface_integral/interface/subface_integrator.hpp +++ b/surface_integral/interface/subface_integrator.hpp @@ -5,8 +5,6 @@ #include #include - - namespace internal { @@ -39,21 +37,22 @@ public: ~integrator_t() = default; - - double calculate(int gauss_order, double (*func)(double u, double v, - const Eigen::Vector3d& p, - const Eigen::Vector3d& dU, - const Eigen::Vector3d& dV)) const; - + double calculate(int gauss_order, + double (*func)(double u, + double v, + const Eigen::Vector3d& p, + const Eigen::Vector3d& dU, + const Eigen::Vector3d& dV)) const; double calculate_one_subface(const subface& subface, const parametric_plane_t& param_plane, int gauss_order, - double (*func)(double u, double v, - const Eigen::Vector3d& p, - const Eigen::Vector3d& dU, - const Eigen::Vector3d& dV)) const; + double (*func)(double u, + double v, + const Eigen::Vector3d& p, + const Eigen::Vector3d& dU, + const Eigen::Vector3d& dV)) const; private: /// Non-owning reference to the list of subfaces @@ -81,15 +80,10 @@ double newton_method(const std::function double tolerance = 1e-8, int max_iterations = 100); -// 例如：计算面积元素 ||dU × dV|| -double area_integrand(double u, double v, - const Eigen::Vector3d& p, - const Eigen::Vector3d& dU, - const Eigen::Vector3d& dV) +double area_integrand(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV) { return 1.0; } - } // namespace internal \ No newline at end of file diff --git a/surface_integral/src/subface_integrator.cpp b/surface_integral/src/subface_integrator.cpp index af5bdc5..8d3ef51 100644 --- a/surface_integral/src/subface_integrator.cpp +++ b/surface_integral/src/subface_integrator.cpp @@ -8,30 +8,30 @@ namespace internal { -double integrator_t::calculate(int gauss_order, double (*func)(double u, double v, - const Eigen::Vector3d& p, - const Eigen::Vector3d& dU, - const Eigen::Vector3d& dV)) const +double integrator_t::calculate( + int gauss_order, + double (*func)(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)) const { double 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, func); + const auto& subface = m_subfaces[subface_index].object_ptr.get(); + total_integral += calculate_one_subface(subface, param_plane, gauss_order, func); } return total_integral; } - -double integrator_t::calculate_one_subface(const subface& subface, const parametric_plane_t& param_plane, int gauss_order, double (*func)(double u, double v, - const Eigen::Vector3d& p, - const Eigen::Vector3d& dU, - const Eigen::Vector3d& dV)) const +double integrator_t::calculate_one_subface( + const subface& subface, + const parametric_plane_t& param_plane, + int gauss_order, + double (*func)(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)) const { - auto solver = subface.fetch_solver_evaluator(); + auto solver = subface.fetch_solver_evaluator(); // Gaussian integration in u direction auto u_integrand = [&](double u) { // Find exact v intersections for each u - stl_vector_mp v_breaks = find_v_intersections_at_u(subface, param_plane, u);; + stl_vector_mp v_breaks = find_v_intersections_at_u(subface, param_plane, u); + ; // Gaussian integration in v direction auto v_integrand = [&](double v) { @@ -39,25 +39,25 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet if (is_point_inside_domain(subface, param_plane, u, v)) { try { // 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 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 + 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>(); // 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 + 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 // Area element: ||dU × dV|| - Eigen::Vector3d cross = dU.cross(dV); - double jacobian = cross.norm(); // Jacobian (area scaling factor) + Eigen::Vector3d cross = dU.cross(dV); + double jacobian = cross.norm(); // Jacobian (area scaling factor) return func(u, v, p, dU, dV) * jacobian; } catch (...) { - return 0.0; // Skip singular points + return 0.0; // Skip singular points } } - return 0.0; // Point not in domain + return 0.0; // Point not in domain }; double v_integral = 0.0; @@ -78,16 +78,16 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet }; // Integrate in u direction - const auto& u_breaks = compute_u_breaks(param_plane,0.0, 1.0); - double integral = 0.0; + const auto& u_breaks = compute_u_breaks(param_plane, 0.0, 1.0); + 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_v_intersections_at_u(subface, param_plane,mid_u); + auto v_intersections = find_v_intersections_at_u(subface, param_plane, mid_u); - if (!v_intersections.empty()) { + if (!v_intersections.empty()) { integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); } } @@ -98,15 +98,12 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet /** * @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 + * 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 integrator_t::compute_u_breaks( - const parametric_plane_t& param_plane, - double u_min, - double u_max) const +stl_vector_mp integrator_t::compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max) const { std::set break_set; @@ -119,41 +116,39 @@ stl_vector_mp integrator_t::compute_u_breaks( 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 + 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()); } - -stl_vector_mp integrator_t::find_v_intersections_at_u( - const subface& subface, - const parametric_plane_t& param_plane, - 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 { 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 + 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) + 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) + 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 { + 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; + if (u > u_val + eps) return 1; return 0; }; @@ -162,9 +157,7 @@ stl_vector_mp integrator_t::find_v_intersections_at_u( 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; - } + if (pos1 * pos2 > 0) { continue; } // Case 2: Both endpoints lie exactly on u = u_val → add both v-values if (pos1 == 0 && pos2 == 0) { @@ -178,12 +171,11 @@ stl_vector_mp integrator_t::find_v_intersections_at_u( intersections.push_back(v1_val); intersections.push_back(v2_val); } - } - else { + } 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 + 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); @@ -192,7 +184,7 @@ stl_vector_mp integrator_t::find_v_intersections_at_u( // 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)); + auto full_res = solver_evaluator(std::move(temp_res)); return std::get(full_res); }; @@ -209,9 +201,7 @@ stl_vector_mp integrator_t::find_v_intersections_at_u( } // Final step: sort and remove duplicates within tolerance - if (!intersections.empty()) { - sort_and_unique_with_tol(intersections, 1e-8); - } + if (!intersections.empty()) { sort_and_unique_with_tol(intersections, 1e-8); } return intersections; } @@ -223,33 +213,28 @@ stl_vector_mp integrator_t::find_v_intersections_at_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 integrator_t::is_point_inside_domain( - const subface& subface, - const parametric_plane_t& param_plane, - 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 { - auto intersections = find_v_intersections_at_u(subface, param_plane, u); - const double tol_near = 1e-8; - const double tol_above = 1e-12; + 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) { + } else if (std::abs(diff) < tol_near) { return true; // on boundary → inside } } return (count % 2) == 1; } -bool integrator_t::is_u_near_singularity(double u, double tol) const -{ - return false; -} +bool integrator_t::is_u_near_singularity(double u, double tol) const { return false; } void integrator_t::sort_and_unique_with_tol(stl_vector_mp& vec, double epsilon) const { @@ -269,17 +254,16 @@ void integrator_t::sort_and_unique_with_tol(stl_vector_mp& vec, double e } // Only accepts functions that return implicit_equation_intermediate -double newton_method( - const std::function& F, - double v_initial, - double tolerance, - int max_iterations) +double newton_method(const std::function& F, + double v_initial, + double tolerance, + int max_iterations) { double v = v_initial; for (int i = 0; i < max_iterations; ++i) { - auto res = F(v); // Known type: implicit_equation_intermediate - double f_val = res.f; + auto res = F(v); // Known type: implicit_equation_intermediate + double f_val = res.f; double df_val = res.df; if (std::abs(f_val) < tolerance) {