From 3154fe8fc3261770685360f5510e5e3c81a3d22b Mon Sep 17 00:00:00 2001 From: mckay Date: Thu, 17 Jul 2025 20:12:47 +0800 Subject: [PATCH] ruff init surface integral module --- .../interface/SurfaceIntegrator.hpp | 69 ++++ surface_integral/src/SurfaceIntegrator.cpp | 294 ++++++++++++++++++ 2 files changed, 363 insertions(+) create mode 100644 surface_integral/interface/SurfaceIntegrator.hpp create mode 100644 surface_integral/src/SurfaceIntegrator.cpp diff --git a/surface_integral/interface/SurfaceIntegrator.hpp b/surface_integral/interface/SurfaceIntegrator.hpp new file mode 100644 index 0000000..40f036c --- /dev/null +++ b/surface_integral/interface/SurfaceIntegrator.hpp @@ -0,0 +1,69 @@ +#pragma once + +#include + +namespace internal +{ + +// 每个活动子面应具有以下结构 +struct parametric_plane { + stl_vector_mp chain_vertices{}; // 链顶点 + flat_index_group chains{}; // 链组 + stl_vector_mp singularity_vertices{}; // 奇异顶点，即链的交点 + stl_vector_mp polar_vertices{}; // 极性顶点，即两个连接顶点周围的最小/最大 x/y + stl_vector_mp parallel_start_vertices{}; // 平行起始顶点，即边 {v, v + 1} 平行于 x/y 轴 +}; + +class SurfaceAreaCalculator +{ +public: + // 构造函数，接受对 subface 的引用 + explicit SurfaceAreaCalculator(const subface& surface); + [[deprecated("Use calculate_new() instead")]] + SurfaceAreaCalculator(const subface& surface, + const stl_vector_map& u_breaks, + double umin, + double umax, + double vmin, + double vmax); + SurfaceAreaCalculator(const subface& surface, const parametric_plane& uv_plane); + + + // 设置 u_breaks（裁剪/分割线） + void set_ubreaks(const stl_vector_map& u_breaks); + + // 计算面积主函数 + template + double calculate(Func&& func, int gauss_order = 3) const; + +private: + const subface& m_surface; // 引用原始曲面 + stl_vector_map m_u_breaks; // 分割线信息（可选） + parametric_plane m_uv_plane; + double Umin = 0.0; // 参数域范围 + double Umax = 1.0; + double Vmin = 0.0; + double Vmax = 1.0; + + // 私有辅助函数（可扩展） + template + double SurfaceAreaCalculator::GaussIntegrate1D(double a, double b, Func&& func, int q) const; + + // 直线u=u_val与边界的交点 + std::vector FindVerticalIntersectionsOCCT(double u_val, + const std::vector& edges, + double v_min, + double v_max); + + bool IsPointInsideDomain(double u, + double v, + const std::vector& outerEdges, + const std::vector& innerEdges, + double u_min, + double u_max, + double v_min, + double v_max); + double integrate_over_uv(int num_samples) const; +}; + +} // namespace internal \ No newline at end of file diff --git a/surface_integral/src/SurfaceIntegrator.cpp b/surface_integral/src/SurfaceIntegrator.cpp new file mode 100644 index 0000000..f96eb62 --- /dev/null +++ b/surface_integral/src/SurfaceIntegrator.cpp @@ -0,0 +1,294 @@ +#include "surface_integral.hpp" // Replace with your actual header path +#include // For vector and cross product operations +#include // For math functions like sqrt + +namespace internal +{ + +// Constructor 1: Initialize only with a reference to the surface +SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface) : m_surface(surface) {} + +// Constructor 2: Initialize with surface and u-breaks (e.g., trimming curves) +SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface, + const stl_vector_map& u_breaks, + double umin, + double umax, + double vmin, + double vmax) + : m_surface(surface), m_u_breaks(u_breaks), Umin(umin), Umax(umax), Vmin(vmin), Vmax(vmax) +{ +} + +SurfaceAreaCalculator::SurfaceAreaCalculator(const subface& surface, const parametric_plane& uv_plane) + : m_surface(surface), m_uv_plane(uv_plane), + Umin(0.0), Umax(0.0), Vmin(0.0), Vmax(0.0) { + + if (!uv_plane.chain_vertices.empty()) { + // 初始化为第一个点的坐标 + double min_u = uv_plane.chain_vertices[0].x(); + double max_u = uv_plane.chain_vertices[0].x(); + double min_v = uv_plane.chain_vertices[0].y(); + double max_v = uv_plane.chain_vertices[0].y(); + + // 遍历所有链顶点 + for (const auto& pt : uv_plane.chain_vertices) { + double u = pt.x(); + double v = pt.y(); + + if (u < min_u) min_u = u; + if (u > max_u) max_u = u; + if (v < min_v) min_v = v; + if (v > max_v) max_v = v; + } + + Umin = min_u; + Umax = max_u; + Vmin = min_v; + Vmax = max_v; + } else { + // 没有顶点时使用默认范围 [0, 1] × [0, 1] + Umin = 0.0; + Umax = 1.0; + Vmin = 0.0; + Vmax = 1.0; + } + + std::set unique_vertex_indices; + + // 插入所有类型的顶点索引 + unique_vertex_indices.insert(uv_plane.singularity_vertices.begin(), uv_plane.singularity_vertices.end()); + unique_vertex_indices.insert(uv_plane.polar_vertices.begin(), uv_plane.polar_vertices.end()); + unique_vertex_indices.insert(uv_plane.parallel_start_vertices.begin(), uv_plane.parallel_start_vertices.end()); + + std::set unique_u_values; + + for (uint32_t idx : unique_vertex_indices) { + if (idx < uv_plane.chain_vertices.size()) { + double u = uv_plane.chain_vertices[idx].x(); + unique_u_values.insert(u); + } + } + + // 转换为 vector 并排序（set 默认是有序的） + m_u_breaks = std::vector(unique_u_values.begin(), unique_u_values.end()); +} + +// Set u-breaks (optional trimming or partitioning lines) +void SurfaceAreaCalculator::set_ubreaks(const stl_vector_map& u_breaks) { m_u_breaks = u_breaks; } + +// Main entry point to compute surface area +template +double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const +{ + // 在u方向进行高斯积分 + auto u_integrand = [&](double u) { + // 对每个u，找到v方向的精确交点 + std::vector v_breaks = {self.Vmin, self.Vmax}; + + auto outer_intersections = FindVerticalIntersectionsOCCT(u, self.outerEdges, self.Vmin, self.Vmax); + v_breaks.insert(v_breaks.end(), outer_intersections.begin(), outer_intersections.end()); + + auto inner_intersections = FindVerticalIntersectionsOCCT(u, self.innerEdges, self.Vmin, self.Vmax); + v_breaks.insert(v_breaks.end(), inner_intersections.begin(), inner_intersections.end()); + + // 排序并去重 + sort_and_unique_with_tol(v_breaks, 1e-6); + + // 在v方向进行高斯积分 + auto v_integrand = [&](double v) { + // 判断点是否在有效域内 + if (IsPointInsideself(u, + v, + self.outerEdges, + self.innerEdges, + self.Umin, + self.Umax, + self.Vmin, + self.Vmax)) { + try { + gp_Pnt p; + gp_Vec dU, dV; + surface->D1(u, v, p, dU, dV); + + const double jacobian = dU.Crossed(dV).Magnitude(); + return func(u, v, p, dU, dV) * jacobian; + } catch (...) { + return 0.0; // 跳过奇异点 + } + } + return 0.0; // 点不在有效域内 + }; + + double v_integral = 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 (IsPointInsideself(u, + mid_v, + self.outerEdges, + self.innerEdges, + self.Umin, + self.Umax, + self.Vmin, + self.Vmax)) { + v_integral += GaussIntegrate1D(a, b, v_integrand, gauss_order); + } else { + std::cout << "uv out of domain: (" << u << "," << mid_v << ")" << std::endl; + } + } + + return v_integral; + }; + + // 在u方向积分 + 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 = FindVerticalIntersectionsOCCT(mid_u, self.outerEdges, self.Vmin, self.Vmax); + + if (!v_intersections.empty()) { // 确保该u区间有有效区域 + + double integralp = GaussIntegrate1D(a, b, u_integrand, gauss_order); + integral += integralp; + std::cout << "integral " << i << ": " << integralp << std::endl; + } + } + + return integral; +} + +// 在区间[a,b]上进行高斯积分采样 +template +double SurfaceAreaCalculator::GaussIntegrate1D(double a, double b, Func&& func, int q) const { + + double sum = 0.0; + for (int i = 0; i < q; ++i) { + double x=a+(b-a)*GaussQuad::x(q, i); + double w=GaussQuad::w(q, i); + sum += w * func(x); + } + return sum * (b-a); +} + +// 直线u=u_val与边界的交点 +std::vector SurfaceAreaCalculator::FindVerticalIntersectionsOCCT( + double u_val, + const std::vector& edges, + double v_min, + double v_max) +{ + std::vector intersections; + + // 创建垂直线(u=u_val) + gp_Lin2d vertical_line(gp_Pnt2d(u_val, v_min), gp_Dir2d(0, 1)); + Handle(Geom2d_Line) hVerticalLine = new Geom2d_Line(vertical_line); + + for (const auto& edge : edges) { + if (edge.IsNull()) { + continue; + } + + // 使用Geom2dAPI_InterCurveCurve求交 + Geom2dAPI_InterCurveCurve intersector(edge, hVerticalLine); + + if (intersector.NbPoints() > 0) { + for (int i = 1; i <= intersector.NbPoints(); ++i) { + gp_Pnt2d pt = intersector.Point(i); + double v = pt.Y(); + if (v >= v_min && v <= v_max) { + intersections.push_back(v); + } + } + } + + // 补充采样法确保可靠性 + Geom2dAdaptor_Curve adaptor(edge); + double first = adaptor.FirstParameter(); + double last = adaptor.LastParameter(); + + const int samples = 20; + gp_Pnt2d prev_p; + adaptor.D0(first, prev_p); + double prev_v = prev_p.Y(); + + for (int i = 1; i <= samples; ++i) { + double param = first + i * (last - first) / samples; + gp_Pnt2d curr_p; + adaptor.D0(param, curr_p); + double curr_v = curr_p.Y(); + + if ((prev_v - u_val) * (curr_v - u_val) <= 0) { + // 牛顿迭代法 + double t = param - (last - first)/samples; + for (int iter = 0; iter < 5; ++iter) { + gp_Pnt2d p; + gp_Vec2d deriv; + adaptor.D1(t, p, deriv); + if (std::abs(p.X() - u_val) < Precision::Confusion()) { + if (p.Y() >= v_min && p.Y() <= v_max) { + intersections.push_back(p.Y()); + } + break; + } + t -= (p.X() - u_val) / deriv.X(); + } + } + prev_v = curr_v; + } + } + + // 去重排序 + std::sort(intersections.begin(), intersections.end()); + intersections.erase(std::unique(intersections.begin(), intersections.end()), intersections.end()); + return intersections; +} + + +// Private method: Numerical integration over the UV parameter domain +double SurfaceAreaCalculator::integrate_over_uv(int num_samples) const +{ + double area = 0.0; + + double du = 1.0 / num_samples; + double dv = 1.0 / num_samples; + + // Fetch evaluators for point and derivative calculation + auto point_func = m_surface.fetch_point_by_param_evaluator(); + auto solver_func = m_surface.fetch_solver_evaluator(); + + for (int i = 0; i < num_samples; ++i) { + double u = (i + 0.5) * du; // Midpoint sampling in u-direction + + for (int j = 0; j < num_samples; ++j) { + double v = (j + 0.5) * dv; // Midpoint sampling in v-direction + + // Evaluate the surface point at (u, v) + Eigen::Vector4d pt = point_func(Eigen::Vector2d(u, v)); + + // Get constraint and solve intermediate result + auto constraint = m_surface.fetch_curve_constraint_evaluator(internal::parameter_v_t{}, v)(u); + auto result = std::get(solver_func(std::move(constraint))); + + // Extract partial derivatives ∂S/∂u and ∂S/∂v + Eigen::Vector3d dSdu = result.grad_f.col(0); + Eigen::Vector3d dSdv = result.grad_f.col(1); + + // Compute the differential area element (cross product norm) + double dA = dSdu.cross(dSdv).norm(); + + // Accumulate area using numerical integration + area += dA * du * dv; + } + } + + return area; +} + +} // namespace internal \ No newline at end of file