From 97320c95c7d65bb948d6571b08c29aae19de5381 Mon Sep 17 00:00:00 2001 From: mckay Date: Wed, 22 Oct 2025 20:11:29 +0800 Subject: [PATCH] fix(integrator): fix small bug in integrator --- .../interface/subface_integrator.hpp | 27 +++++---- surface_integral/src/subface_integrator.cpp | 57 ++++++++++--------- 2 files changed, 47 insertions(+), 37 deletions(-) diff --git a/surface_integral/interface/subface_integrator.hpp b/surface_integral/interface/subface_integrator.hpp index 8bb3b9c..646ebb7 100644 --- a/surface_integral/interface/subface_integrator.hpp +++ b/surface_integral/interface/subface_integrator.hpp @@ -40,9 +40,10 @@ public: const Eigen::Vector3d& dV)) const; - double calculate_one_subface(const subface& subface, - const parametric_plane_t& param_plane, - int gauss_order, + double calculate_one_subface(const subface& subface, + const parametric_plane_t& param_plane, + const stl_vector_mp& chain_bboxes, + int gauss_order, double (*func)(double u, double v, const Eigen::Vector3d& p, @@ -61,15 +62,19 @@ private: stl_vector_mp compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max) const; - void find_v_intersections_at_u(const subface& subface, - const parametric_plane_t& param_plane, - double u_val, - double v_min, - double v_max, - stl_vector_mp& intersections, - stl_vector_mp& intersected_chains) const; + void find_v_intersections_at_u(const subface& subface, + const parametric_plane_t& param_plane, + const stl_vector_mp& chain_bboxes, + double u_val, + double v_min, + double v_max, + stl_vector_mp& intersections, + stl_vector_mp& intersected_chains) const; - bool is_point_inside_domain(const subface& subface, const parametric_plane_t& param_plane, double u, double v) const; + + bool is_edge_inside_domain(const stl_vector_mp>& chain_group_indices, + uint16_t chain_idx1, + uint16_t chain_idx2) const; bool is_u_near_singularity(double u, double tol = 1e-6) const; diff --git a/surface_integral/src/subface_integrator.cpp b/surface_integral/src/subface_integrator.cpp index fc3d60c..21c9132 100644 --- a/surface_integral/src/subface_integrator.cpp +++ b/surface_integral/src/subface_integrator.cpp @@ -4,12 +4,13 @@ #include // For math functions like sqrt #include #include +#include namespace internal { -integrator_t(const stl_vector_mp>& surfaces, - const flat_hash_map_mp& uv_planes) +integrator_t::integrator_t(const stl_vector_mp>& surfaces, + const flat_hash_map_mp& uv_planes) : m_subfaces(surfaces), m_uv_planes(uv_planes) { std::cout << "Integrator initialized with " << surfaces.size() << " subfaces." << std::endl; @@ -18,10 +19,10 @@ integrator_t(const stl_vector_mp>& surfaces, // calculate chain bbox for (int subface_index = 0; subface_index < m_subfaces.size(); ++subface_index) { if (m_uv_planes.find(subface_index) != m_uv_planes.end()) { - const auto& param_plane = m_uv_planes.at(subface_index); - auto& chain_bboxs = m_chain_bboxes_hash[subface_index]; + const auto& param_plane = m_uv_planes.at(subface_index); + auto& chain_bboxes = m_chain_bboxes_hash[subface_index]; chain_bboxes.clear(); - chain_bboxs.reverse(chain_vertices.size()); + chain_bboxes.reserve(param_plane.chain_vertices.size()); for (const auto& chain : param_plane.chain_vertices) { Eigen::AlignedBox2d bbox; bbox.setEmpty(); @@ -33,8 +34,6 @@ integrator_t(const stl_vector_mp>& surfaces, } } } - - m_chain_bboxes_hash } double integrator_t::calculate( @@ -44,9 +43,10 @@ double integrator_t::calculate( double total_integral = 0.0; for (int subface_index = 0; subface_index < m_subfaces.size(); ++subface_index) { if (m_uv_planes.find(subface_index) != m_uv_planes.end()) { - const auto& param_plane = m_uv_planes.at(subface_index); - const auto& subface = m_subfaces[subface_index].object_ptr.get(); - total_integral += calculate_one_subface(subface, param_plane, gauss_order, func); + const auto& param_plane = m_uv_planes.at(subface_index); + const auto& subface = m_subfaces[subface_index].object_ptr.get(); + const auto& chain_bboxes = m_chain_bboxes_hash.at(subface_index); + total_integral += calculate_one_subface(subface, param_plane, chain_bboxes, gauss_order, func); } else { // the subface has no associated parametric plane, meaning it is compete face without trimming } @@ -55,9 +55,10 @@ double integrator_t::calculate( } double integrator_t::calculate_one_subface( - const subface& subface, - const parametric_plane_t& param_plane, - int gauss_order, + const subface& subface, + const parametric_plane_t& param_plane, + const stl_vector_mp& chain_bboxes, + 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(); @@ -69,9 +70,9 @@ double integrator_t::calculate_one_subface( // Gaussian integration in u direction auto u_integrand = [&](double u) { // Find exact v intersections for each u - stl_vector_mp v_intersections; + stl_vector_mp v_breaks; stl_vector_mp intersected_chains; - find_v_intersections_at_u(subface, param_plane, u, v_min, v_max, v_intersections, intersected_chains); + find_v_intersections_at_u(subface, param_plane, chain_bboxes, u, v_min, v_max, v_breaks, intersected_chains); // Gaussian integration in v direction auto v_integrand = [&](double v) { @@ -104,7 +105,7 @@ double integrator_t::calculate_one_subface( // Check midpoint validity double mid_v = (a + b) / 2.0; - if (is_edge_inside_domain(subface, param_plane, u, mid_v)) { + if (is_edge_inside_domain(param_plane.chain_group_indices, intersected_chains[i], intersected_chains[i + 1])) { v_integral += integrate_1D(a, b, v_integrand, gauss_order); } else { std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl; @@ -124,7 +125,7 @@ double integrator_t::calculate_one_subface( double mid_u = (a + b) / 2.0; stl_vector_mp v_intersections; stl_vector_mp intersected_chains; - find_v_intersections_at_u(subface, param_plane, mid_u, v_min, v_max, v_intersections, intersected_chains); + find_v_intersections_at_u(subface, param_plane, chain_bboxes, mid_u, v_min, v_max, v_intersections, intersected_chains); if (!v_intersections.empty()) { integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); @@ -165,7 +166,7 @@ stl_vector_mp integrator_t::compute_u_breaks(const parametric_plane_t& p void integrator_t::find_v_intersections_at_u(const subface& subface, const parametric_plane_t& param_plane, - const stl_vector_mp& chain_bboxs, + const stl_vector_mp& chain_bboxes, double u_val, double v_min, double v_max, @@ -187,9 +188,10 @@ void integrator_t::find_v_intersections_at_u(const subface& } // filter bbox - if (chain_bbox[chain_idx].isEmpty() - || (chain_bbox[chain_idx].x() < u_val - 1e-8 || chain_bbox[chain_idx].x() > u_val + 1e-8)) { - chain_idx++; + const auto& box = chain_bboxes[chain_idx]; + if (box.isEmpty() + || box.max()(0) < u_val - 1e-8 // x_max < u_val - eps + || box.min()(0) > u_val + 1e-8) { // x_min > u_val + eps continue; } @@ -268,14 +270,17 @@ void integrator_t::find_v_intersections_at_u(const subface& // if (!intersections.empty()) { sort_and_unique_with_tol(intersections, 1e-8); } } -bool is_edge_inside_domain(const stl_vector_mp>& chain_group_indices, - uint16_t chain_idx1, - uint16_t chain_idx2) const +bool integrator_t::is_edge_inside_domain(const stl_vector_mp>& chain_group_indices, + uint16_t chain_idx1, + uint16_t chain_idx2) const { + if (chain_idx1 >= chain_group_indices.size() || chain_idx2 >= chain_group_indices.size()) { return false; } + const auto& group1 = chain_group_indices[chain_idx1]; const auto& group2 = chain_group_indices[chain_idx2]; - if (group1.find(chain_idx2) != group1.end() && group2.find(chain_idx1) != group2.end()) { return true; } - return false; + + return (std::find(group1.begin(), group1.end(), chain_idx2) != group1.end()) + && (std::find(group2.begin(), group2.end(), chain_idx1) != group2.end()); } bool integrator_t::is_u_near_singularity(double u, double tol) const { return false; }