diff --git a/frontend/src/implicit_surface_network_solver.cpp b/frontend/src/implicit_surface_network_solver.cpp index 6d73e71..90bed93 100644 --- a/frontend/src/implicit_surface_network_solver.cpp +++ b/frontend/src/implicit_surface_network_solver.cpp @@ -1,8 +1,10 @@ #include -#include +#include EXTERN_C_BEGIN + + void implicit_network_solver::generate_polymesh(stl_vector_mp& output_vertices, stl_vector_mp& output_polygon_faces, stl_vector_mp& output_vertex_counts_of_face) @@ -17,7 +19,12 @@ void implicit_network_solver::generate_polymesh(stl_vector_mp& output_vertex_counts_of_face, output_parameteric_planes); - //integrator_t integrator(const subface& surface, const parametric_plane& uv_plane); + auto integrator = internal::integrator_t(m_blobtree->subfaces, output_parameteric_planes); + + double area = integrator.calculate(4, internal::area_integrand); + + std::cout << "Generated mesh area: " << area << std::endl; + m_timers.pop_timer("generate_polymesh"); } diff --git a/surface_integral/interface/SurfaceIntegrator.hpp b/surface_integral/interface/subface_integrator.hpp similarity index 69% rename from surface_integral/interface/SurfaceIntegrator.hpp rename to surface_integral/interface/subface_integrator.hpp index 3f2fc14..b8688fe 100644 --- a/surface_integral/interface/SurfaceIntegrator.hpp +++ b/surface_integral/interface/subface_integrator.hpp @@ -5,6 +5,8 @@ #include #include + + namespace internal { @@ -27,32 +29,40 @@ public: */ integrator_t(const stl_vector_mp>& surfaces, const flat_hash_map_mp& uv_planes) - : m_surfaces(surfaces), m_uv_planes(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; } /// Default destructor ~integrator_t() = default; - template - double calculate(int gauss_order, Func&& func) 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; + - template double calculate_one_subface(const subface& subface, const parametric_plane_t& param_plane, int gauss_order, - Func&& func) 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 - const stl_vector_mp>& m_surfaces; + const stl_vector_mp>& m_subfaces; /// Non-owning reference to the map of parametric planes (ID -> parametric_plane_t) const flat_hash_map_mp& m_uv_planes; - stl_vector_mp compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max); + 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, const parametric_plane_t& param_plane, @@ -71,4 +81,15 @@ 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) +{ + return 1.0; +} + + + } // namespace internal \ No newline at end of file diff --git a/surface_integral/src/SurfaceIntegrator.cpp b/surface_integral/src/subface_integrator.cpp similarity index 78% rename from surface_integral/src/SurfaceIntegrator.cpp rename to surface_integral/src/subface_integrator.cpp index a044294..af5bdc5 100644 --- a/surface_integral/src/SurfaceIntegrator.cpp +++ b/surface_integral/src/subface_integrator.cpp @@ -1,4 +1,4 @@ -#include "SurfaceIntegrator.hpp" +#include "subface_integrator.hpp" #include "quadrature.hpp" #include // For vector and cross product operations #include // For math functions like sqrt @@ -8,27 +8,30 @@ namespace internal { -// Main entry point to compute surface area -// Main entry point to compute surface area -template -double integrator_t::calculate(int gauss_order, Func&& func) 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 { - total_integral = 0.0; + 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, std::forward(func)); + total_integral += calculate_one_subface(subface, param_plane, gauss_order, func); } return total_integral; } -template -double integrator_t::calculate_one_subface(const subface& subface, const parametric_plane_t& param_plane, int gauss_order, Func&& func) 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(); // Gaussian integration in u direction auto u_integrand = [&](double u) { // Find exact v intersections for each u - std::vector v_breaks = find_vertical_intersections(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) { @@ -75,14 +78,14 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet }; // Integrate in u direction - const auto& u_breaks = compute_u_breaks(param_plane); + 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_vertical_intersections(mid_u, self.outerEdges, self.Vmin, self.Vmax); + auto v_intersections = find_v_intersections_at_u(subface, param_plane,mid_u); if (!v_intersections.empty()) { integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); @@ -92,72 +95,6 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet return integral; } -// 在 integrator_t 类中添加： -/* -double integrator_t::compute_volume(int gauss_order) const -{ - double total_volume = 0.0; - - // 外层：对 u 分段积分 - auto u_integrand = [&](double u) { - std::vector v_breaks = find_vertical_intersections(u); - double v_integral = 0.0; - - // 内层：对 v 积分 - auto v_integrand = [&](double v) -> double { - if (!is_point_inside_domain(subface, u, v)) { - return 0.0; - } - - try { - // 获取偏导数 - auto eval_du = m_surface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // ∂/∂u - auto eval_dv = m_surface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // ∂/∂v - - auto res_u = eval_du(u); - auto res_v = eval_dv(v); - - Eigen::Vector3d p = res_u.f.template head<3>(); // r(u,v) - Eigen::Vector3d dU = res_u.grad_f.template head<3>(); // r_u - Eigen::Vector3d dV = res_v.grad_f.template head<3>(); // r_v - - // 计算 r · (r_u × r_v) - Eigen::Vector3d cross = dU.cross(dV); - double mixed_product = p.dot(cross); - - return mixed_product; // 注意：不是 norm，是点积 - } catch (...) { - return 0.0; - } - }; - - for (size_t i = 0; i < v_breaks.size() - 1; ++i) { - double a = v_breaks[i]; - double b = v_breaks[i + 1]; - double mid_v = (a + b) / 2.0; - if (is_point_inside_domain(subface, u, mid_v)) { - v_integral += integrate_1D(a, b, v_integrand, gauss_order); - } - } - - return v_integral; - }; - - // 在 u 方向积分 - for (size_t i = 0; i < m_u_breaks.size() - 1; ++i) { - double a = m_u_breaks[i]; - double b = m_u_breaks[i + 1]; - double mid_u = (a + b) / 2.0; - auto v_intersections = find_vertical_intersections(mid_u); - if (!v_intersections.empty()) { - total_volume += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); - } - } - - // 乘以 1/3 - return std::abs(total_volume) / 3.0; -}*/ - /** * @brief Compute the u-parameter breakpoints for integration. * @note The function currently uses std::set for uniqueness and sorting, @@ -166,10 +103,10 @@ double integrator_t::compute_volume(int gauss_order) const * 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 compute_u_breaks( +stl_vector_mp integrator_t::compute_u_breaks( const parametric_plane_t& param_plane, double u_min, - double u_max) + double u_max) const { std::set break_set;