2 changed files with 363 additions and 0 deletions
@ -0,0 +1,69 @@ |
|||
#pragma once |
|||
|
|||
#include <primitive_process/interface/base/subface.hpp> |
|||
|
|||
namespace internal |
|||
{ |
|||
|
|||
// 每个活动子面应具有以下结构
|
|||
struct parametric_plane { |
|||
stl_vector_mp<Eigen::Vector2d> chain_vertices{}; // 链顶点
|
|||
flat_index_group chains{}; // 链组
|
|||
stl_vector_mp<uint32_t> singularity_vertices{}; // 奇异顶点，即链的交点
|
|||
stl_vector_mp<uint32_t> polar_vertices{}; // 极性顶点，即两个连接顶点周围的最小/最大 x/y
|
|||
stl_vector_mp<uint32_t> 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<double>& 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<double>& u_breaks); |
|||
|
|||
// 计算面积主函数
|
|||
template <typename Func> |
|||
double calculate(Func&& func, int gauss_order = 3) const; |
|||
|
|||
private: |
|||
const subface& m_surface; // 引用原始曲面
|
|||
stl_vector_map<double> 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 <typename Func> |
|||
double SurfaceAreaCalculator::GaussIntegrate1D(double a, double b, Func&& func, int q) const; |
|||
|
|||
// 直线u=u_val与边界的交点
|
|||
std::vector<double> FindVerticalIntersectionsOCCT(double u_val, |
|||
const std::vector<Handle(Geom2d_Curve)>& edges, |
|||
double v_min, |
|||
double v_max); |
|||
|
|||
bool IsPointInsideDomain(double u, |
|||
double v, |
|||
const std::vector<Handle(Geom2d_Curve)>& outerEdges, |
|||
const std::vector<Handle(Geom2d_Curve)>& innerEdges, |
|||
double u_min, |
|||
double u_max, |
|||
double v_min, |
|||
double v_max); |
|||
double integrate_over_uv(int num_samples) const; |
|||
}; |
|||
|
|||
} // namespace internal
|
|||
@ -0,0 +1,294 @@ |
|||
#include "surface_integral.hpp" // Replace with your actual header path
|
|||
#include <Eigen/Geometry> // For vector and cross product operations |
|||
#include <cmath> // 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<double>& 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<uint32_t> 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<double> 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<double>(unique_u_values.begin(), unique_u_values.end()); |
|||
} |
|||
|
|||
// Set u-breaks (optional trimming or partitioning lines)
|
|||
void SurfaceAreaCalculator::set_ubreaks(const stl_vector_map<double>& u_breaks) { m_u_breaks = u_breaks; } |
|||
|
|||
// Main entry point to compute surface area
|
|||
template <typename Func> |
|||
double SurfaceAreaCalculator::calculate(Func&& func, int gauss_order) const |
|||
{ |
|||
// 在u方向进行高斯积分
|
|||
auto u_integrand = [&](double u) { |
|||
// 对每个u，找到v方向的精确交点
|
|||
std::vector<double> 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<typename Func> |
|||
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<double> SurfaceAreaCalculator::FindVerticalIntersectionsOCCT( |
|||
double u_val, |
|||
const std::vector<Handle(Geom2d_Curve)>& edges, |
|||
double v_min, |
|||
double v_max) |
|||
{ |
|||
std::vector<double> 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<internal::parametric_equation_intermediate>(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
|
|||
Loading…
Reference in new issue