diff --git a/surface_integral/interface/subface_integrator.hpp b/surface_integral/interface/subface_integrator.hpp index 63d3809..4ff62a3 100644 --- a/surface_integral/interface/subface_integrator.hpp +++ b/surface_integral/interface/subface_integrator.hpp @@ -26,12 +26,7 @@ public: * throughout the lifetime of this integrator. */ 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; - std::cout << "m_uv_planes size: " << m_uv_planes.size() << std::endl; - } + const flat_hash_map_mp& uv_planes); /// Default destructor ~integrator_t() = default; @@ -61,14 +56,19 @@ private: /// Non-owning reference to the map of parametric planes (ID -> parametric_plane_t) const flat_hash_map_mp& m_uv_planes; + // Precomputed bounding boxes for each chain in each parametric plane + flat_hash_map_mp> m_chain_bboxes_hash{}; + stl_vector_mp compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max) const; - stl_vector_mp find_v_intersections_at_u(const subface& subface, + void find_v_intersections_at_u(const subface& subface, const parametric_plane_t& param_plane, double u_val, - bool if_add_boundary=true, - double v_min=0.0, - double v_max=0.0) const; + 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; diff --git a/surface_integral/src/subface_integrator.cpp b/surface_integral/src/subface_integrator.cpp index a3999b4..85b9d66 100644 --- a/surface_integral/src/subface_integrator.cpp +++ b/surface_integral/src/subface_integrator.cpp @@ -8,6 +8,37 @@ namespace internal { +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; + std::cout << "m_uv_planes size: " << m_uv_planes.size() << std::endl; + + // 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]; + chain_bboxes.clear(); + chain_bboxs.reverse(chain_vertices.size()); + for (const auto& chain: param_plane.chain_vertices){ + Eigen::AlignedBox2d bbox; + bbox.setEmpty(); + + if (!chain.empty()) { + for (const auto& pt : chain) { + bbox.extend(pt); + } + } + chain_bboxes.push_back(bbox); + } + } + } + + m_chain_bboxes_hash +} + 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 @@ -40,33 +71,32 @@ 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_breaks = find_v_intersections_at_u(subface, param_plane, u, true,v_min, v_max); + stl_vector_mp v_intersections; + stl_vector_mp intersected_chains; + find_v_intersections_at_u(subface, param_plane, u, v_min, v_max, v_intersections, intersected_chains); // Gaussian integration in v direction auto v_integrand = [&](double v) { // Check if point is inside domain - 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 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 - - // Area element: ||dU × dV|| - 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 - } + 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 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 + + // Area element: ||dU × dV|| + 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; // Point not in domain }; double v_integral = 0.0; @@ -76,7 +106,7 @@ double integrator_t::calculate_one_subface( // Check midpoint validity double mid_v = (a + b) / 2.0; - if (is_point_inside_domain(subface, param_plane, u, mid_v)) { + if (is_edge_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; @@ -94,7 +124,9 @@ double integrator_t::calculate_one_subface( 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,true, v_min, v_max); + 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); if (!v_intersections.empty()) { integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); @@ -133,28 +165,37 @@ stl_vector_mp integrator_t::compute_u_breaks(const parametric_plane_t& p return stl_vector_mp(break_set.begin(), break_set.end()); } -stl_vector_mp integrator_t::find_v_intersections_at_u(const subface& subface, +void integrator_t::find_v_intersections_at_u(const subface& subface, const parametric_plane_t& param_plane, + const stl_vector_mp& chain_bboxs, double u_val, - bool if_add_boundary, double v_min, - double v_max + double v_max, + stl_vector_mp& intersections, + stl_vector_mp& intersected_chains ) const { - stl_vector_mp intersections; - if (if_add_boundary) { - intersections.push_back(v_min); - intersections.push_back(v_max); - } + intersections.clear(); + intersected_chains.clear(); + + uint16_t chain_idx =0; // Iterate over each boundary chain for (const auto& chain : param_plane.chain_vertices) { const size_t n_vertices = chain.size()-1; if (n_vertices < 2){ std::cout << "Warning: Degenerate chain with less than 2 vertices." << std::endl; + chain_idx++; continue; // Skip degenerate chains } + // filter bbox + if (chain_bbox[chain_idx].isEmpty() + || (chain_bbox[chain_idx].x()u_val+1e-8)) { + chain_idx++; + continue; + } + // Iterate over each edge in the chain (including closing edge if desired) for (size_t i = 0; i < n_vertices; ++i) { const Eigen::Vector2d& v1 = chain[i]; // Current vertex: (u1, v1) @@ -182,6 +223,8 @@ stl_vector_mp integrator_t::find_v_intersections_at_u(const subface& if (pos1 == 0 && pos2 == 0) { intersections.push_back(v1_val); intersections.push_back(v2_val); + intersected_chains.push_back(chain_idx); + intersected_chains.push_back(chain_idx); } // Case 3: One endpoint on u = u_val or segment crosses u = u_val else if (std::abs(u1 - u2) < eps) { @@ -189,6 +232,8 @@ stl_vector_mp integrator_t::find_v_intersections_at_u(const subface& if (pos1 == 0 || pos2 == 0) { intersections.push_back(v1_val); intersections.push_back(v2_val); + intersected_chains.push_back(chain_idx); + intersected_chains.push_back(chain_idx); } } else { // General case: non-vertical segment crossing u = u_val @@ -211,48 +256,37 @@ stl_vector_mp integrator_t::find_v_intersections_at_u(const subface& try { double v_solution = newton_method(target_function, v_initial); intersections.push_back(v_solution); + intersected_chains.push_back(chain_idx); } catch (...) { // If Newton's method fails (e.g., divergence), fall back to linear approximation intersections.push_back(v_initial); + intersected_chains.push_back(chain_idx); } } } + chain_idx ++; } // 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; } -/* - 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 integrator_t::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 { - auto intersections = find_v_intersections_at_u(subface, param_plane, u, false); - 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 - } + 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 (count % 2) == 1; + 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