You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
286 lines
12 KiB
286 lines
12 KiB
#include "SurfaceIntegrator.hpp"
|
|
#include "quadrature.hpp"
|
|
#include <Eigen/Geometry> // For vector and cross product operations
|
|
#include <cmath> // For math functions like sqrt
|
|
#include <set>
|
|
|
|
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_mp<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 = stl_vector_mp<double>(unique_u_values.begin(), unique_u_values.end());
|
|
}
|
|
|
|
// Set u-breaks (optional trimming or partitioning lines)
|
|
void SurfaceAreaCalculator::set_ubreaks(const stl_vector_mp<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 = find_vertical_intersections(u);
|
|
|
|
// 在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 += gauss_integrate_1D(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 = find_vertical_intersections(mid_u, self.outerEdges, self.Vmin, self.Vmax);
|
|
|
|
if (!v_intersections.empty()) { // 确保该u区间有有效区域
|
|
|
|
double integralp = gauss_integrate_1D(a, b, u_integrand, gauss_order);
|
|
integral += integralp;
|
|
std::cout << "integral " << i << ": " << integralp << std::endl;
|
|
}
|
|
}
|
|
|
|
return integral;
|
|
}
|
|
|
|
// 直线u=u_val与边界的交点
|
|
std::vector<double> SurfaceAreaCalculator::find_vertical_intersections(double u_val)
|
|
{
|
|
std::vector<double> intersections;
|
|
std::vector<int32_t> uPositionFlags;
|
|
uPositionFlags.reserve(m_uv_plane.chain_vertices.size());
|
|
|
|
std::transform(m_uv_plane.chain_vertices.begin(),
|
|
m_uv_plane.chain_vertices.end(),
|
|
std::back_inserter(uPositionFlags),
|
|
[&](const auto& currentVertex) {
|
|
double uDifference = currentVertex.x() - u_val;
|
|
if (uDifference < 0) return -1; // 在参考值左侧
|
|
if (uDifference > 0) return 1; // 在参考值右侧
|
|
return 0; // 等于参考值
|
|
});
|
|
|
|
uint32_t group_idx = 0;
|
|
for (uint32_t element_idx = 0; element_idx < m_uv_plane.chains.index_group.size() - 1; element_idx++) {
|
|
if (element_idx > m_uv_plane.chains.start_indices[group_idx + 1]) group_idx++;
|
|
if (element_idx + 1 < m_uv_plane.chains.start_indices[group_idx + 1]) {
|
|
uint32_t vertex_idx1 = m_uv_plane.chains.index_group[element_idx];
|
|
uint32_t vertex_idx2 = m_uv_plane.chains.index_group[element_idx + 1];
|
|
if (uPositionFlags[vertex_idx1] * uPositionFlags[vertex_idx2] <= 0) {
|
|
auto v1 = m_uv_plane.chain_vertices[vertex_idx1], v2 = m_uv_plane.chain_vertices[vertex_idx2];
|
|
if (v1.x() != v2.x()) {
|
|
// "The line segment is vertical (u₁ == u₂), so there is no unique v value corresponding to the given u."
|
|
double v_initial = v1.y() + (v2.y() - v1.y()) * (u_val - v1.x()) / (v2.x() - v1.x());
|
|
auto curve_evaluator = m_surface.fetch_curve_constraint_evaluator(u_val);
|
|
auto solver_evaluator = m_surface.fetch_solver_evaluator();
|
|
auto target_function = [&](double v) -> equation_intermediate_t {
|
|
constraint_curve_intermediate temp_res = curve_evaluator(v);
|
|
return solver_evaluator(std::move(temp_res));
|
|
};
|
|
double v_solution = newton_method(target_function, v_initial);
|
|
intersections.push_back(v_solution);
|
|
} else {
|
|
intersections.push_back(v1.y());
|
|
intersections.push_back(v2.y());
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// 去重排序
|
|
sort_and_unique_with_tol(intersections);
|
|
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 is_point_inside_domain(double u, double v)
|
|
{
|
|
uint32_t group_idx = 0, intersection_count = 0;
|
|
for (uint32_t element_idx = 0; element_idx < m_uv_plane.chains.index_group.size() - 1; element_idx++) {
|
|
if (element_idx > m_uv_plane.chains.start_indices[group_idx + 1]) group_idx++;
|
|
if (element_idx + 1 < m_uv_plane.chains.start_indices[group_idx + 1]) {
|
|
uint32_t vertex_idx1 = m_uv_plane.chains.index_group[element_idx];
|
|
uint32_t vertex_idx2 = m_uv_plane.chains.index_group[element_idx + 1];
|
|
auto v1 = m_uv_plane.chain_vertices[vertex_idx1], v2 = m_uv_plane.chain_vertices[vertex_idx2];
|
|
if ((v1.x() <= u && v2.x() > u) || (v2.x() < u && v1.x() >= u)) {
|
|
|
|
double v_intersected = v1.y() + (v2.y() - v1.y()) * (u_val - v1.x()) / (v2.x() - v1.x());
|
|
if (v_interdected - v >= 1e-6) { intersection_count++; }
|
|
else if (std::abs(v_intersected - v) < 1e-6) {
|
|
auto curve_evaluator = m_surface.fetch_curve_constraint_evaluator(u);
|
|
auto solver_evaluator = m_surface.fetch_solver_evaluator();
|
|
auto target_function = [&](double v) -> equation_intermediate_t {
|
|
constraint_curve_intermediate temp_res = curve_evaluator(v);
|
|
return solver_evaluator(std::move(temp_res));
|
|
};
|
|
double v_solution = newton_method(target_function, v_initial);
|
|
if (std::abs(v_solution - v) > 0) {
|
|
intersection_count++;
|
|
}
|
|
}
|
|
}
|
|
/*
|
|
case v1.x() == v2.x() == u, do not count as intersection.but will cout in next iteration
|
|
*/
|
|
}
|
|
}
|
|
return intersection_count % 2 == 1; // in domain
|
|
}
|
|
|
|
|
|
void SurfaceAreaCalculator::sort_and_unique_with_tol(std::vector<double>& vec, double epsilon)
|
|
{
|
|
if (vec.empty()) return;
|
|
|
|
std::sort(vec.begin(), vec.end());
|
|
|
|
size_t write_index = 0;
|
|
for (size_t read_index = 1; read_index < vec.size(); ++read_index) {
|
|
if (std::fabs(vec[read_index] - vec[write_index]) > epsilon) {
|
|
++write_index;
|
|
vec[write_index] = vec[read_index];
|
|
}
|
|
}
|
|
|
|
vec.resize(write_index + 1);
|
|
}
|
|
|
|
// 牛顿法求解器
|
|
double newton_method(const std::function<equation_intermediate_t(double)>& F,
|
|
double v_initial,
|
|
double tolerance,
|
|
int max_iterations)
|
|
{
|
|
double v = v_initial;
|
|
|
|
for (int i = 0; i < max_iterations; ++i) {
|
|
equation_intermediate_t res = F(v);
|
|
double f = res.f;
|
|
double df = res.df;
|
|
|
|
if (std::abs(f) < tolerance) {
|
|
std::cout << "✅ Converged in " << i + 1 << " iterations. v = " << v << std::endl;
|
|
return v;
|
|
}
|
|
|
|
if (std::abs(df) < 1e-10) {
|
|
std::cerr << "⚠️ Derivative near zero. No convergence." << std::endl;
|
|
return v;
|
|
}
|
|
|
|
v -= f / df;
|
|
|
|
std::cout << "Iteration " << i + 1 << ": v = " << v << ", f = " << f << std::endl;
|
|
}
|
|
|
|
std::cerr << "❌ Did not converge within " << max_iterations << " iterations." << std::endl;
|
|
return v;
|
|
}
|
|
|
|
|
|
} // namespace internal
|