diff --git a/surface_integral/interface/subface_integrator.hpp b/surface_integral/interface/subface_integrator.hpp index 8edcb5d..63d3809 100644 --- a/surface_integral/interface/subface_integrator.hpp +++ b/surface_integral/interface/subface_integrator.hpp @@ -65,7 +65,10 @@ private: stl_vector_mp find_v_intersections_at_u(const subface& subface, const parametric_plane_t& param_plane, - double u_val) const; + double u_val, + bool if_add_boundary=true, + double v_min=0.0, + double v_max=0.0) 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 8d3ef51..a3999b4 100644 --- a/surface_integral/src/subface_integrator.cpp +++ b/surface_integral/src/subface_integrator.cpp @@ -13,9 +13,14 @@ double integrator_t::calculate( 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); + 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); + } else { + // the subface has no associated parametric plane, meaning it is compete face without trimming + } } return total_integral; } @@ -27,11 +32,15 @@ double integrator_t::calculate_one_subface( 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(); + const auto& uv_bounds = param_plane.uv_bounds; + auto u_min = uv_bounds.min().x(); + auto u_max = uv_bounds.max().x(); + auto v_min = uv_bounds.min().y(); + auto v_max = uv_bounds.max().y(); // 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, true,v_min, v_max); // Gaussian integration in v direction auto v_integrand = [&](double v) { @@ -78,14 +87,14 @@ double integrator_t::calculate_one_subface( }; // Integrate in u direction - const auto& u_breaks = compute_u_breaks(param_plane, 0.0, 1.0); + const auto& u_breaks = compute_u_breaks(param_plane, u_min, u_max); 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,true, v_min, v_max); if (!v_intersections.empty()) { integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); @@ -115,8 +124,9 @@ stl_vector_mp integrator_t::compute_u_breaks(const parametric_plane_t& p for (size_t i = 0; i < param_plane.chain_vertices.size(); ++i) { const auto& vertices = param_plane.chain_vertices[i]; auto& vertex_flags = param_plane.vertex_special_flags[i]; + auto& edge_flags = param_plane.edge_near_parallel_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] || edge_flags[j]) { break_set.insert(vertices[j].x()); } // Special vertex u } } // Return as vector (sorted and unique due to set) @@ -125,21 +135,30 @@ stl_vector_mp integrator_t::compute_u_breaks(const parametric_plane_t& p stl_vector_mp integrator_t::find_v_intersections_at_u(const subface& subface, const parametric_plane_t& param_plane, - double u_val) const + double u_val, + bool if_add_boundary, + double v_min, + double v_max + ) const { stl_vector_mp intersections; + if (if_add_boundary) { + intersections.push_back(v_min); + intersections.push_back(v_max); + } // 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 + const size_t n_vertices = chain.size()-1; + if (n_vertices < 2){ + std::cout << "Warning: Degenerate chain with less than 2 vertices." << std::endl; + 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) - const Eigen::Vector2d& v1 = chain[i]; // Current vertex: (u1, v1) - const Eigen::Vector2d& v2 = chain[j]; // Next vertex: (u2, v2) + const Eigen::Vector2d& v2 = chain[i + 1]; // Next vertex: (u2, v2) double u1 = v1.x(), v1_val = v1.y(); double u2 = v2.x(), v2_val = v2.y(); @@ -218,7 +237,7 @@ bool integrator_t::is_point_inside_domain(const subface& subface, double u, double v) const { - auto intersections = find_v_intersections_at_u(subface, param_plane, u); + auto intersections = find_v_intersections_at_u(subface, param_plane, u, false); const double tol_near = 1e-8; const double tol_above = 1e-12;