ZCWang
/
ImplicitSurfaceNetwork
1
0
Fork 0
Code Issues 1 Pull Requests Projects Releases Wiki Activity
Browse Source

style(integrator): format docs of integrator hpp and cpp

feat-integrator
mckay 3 months ago
parent
edcc2f95e0
commit
b77f45ee5c
2 changed files with 74 additions and 96 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 30
      surface_integral/interface/subface_integrator.hpp
  2. 140
      surface_integral/src/subface_integrator.cpp

30
surface_integral/interface/subface_integrator.hpp
View File

@ -5,8 +5,6 @@
#include <container/stl_alias.hpp> #include <container/stl_alias.hpp>
#include <container/wrapper/flat_index_group.hpp> #include <container/wrapper/flat_index_group.hpp>
namespace internal namespace internal
{ {
@ -39,21 +37,22 @@ public:
~integrator_t() = default; ~integrator_t() = default;
double calculate(int gauss_order,
double calculate(int gauss_order, double (*func)(double u, double v, double (*func)(double u,
const Eigen::Vector3d& p, double v,
const Eigen::Vector3d& dU, const Eigen::Vector3d& p,
const Eigen::Vector3d& dV)) const; const Eigen::Vector3d& dU,
const Eigen::Vector3d& dV)) const;
double calculate_one_subface(const subface& subface, double calculate_one_subface(const subface& subface,
const parametric_plane_t& param_plane, const parametric_plane_t& param_plane,
int gauss_order, int gauss_order,
double (*func)(double u, double v, double (*func)(double u,
const Eigen::Vector3d& p, double v,
const Eigen::Vector3d& dU, const Eigen::Vector3d& p,
const Eigen::Vector3d& dV)) const; const Eigen::Vector3d& dU,
const Eigen::Vector3d& dV)) const;
private: private:
/// Non-owning reference to the list of subfaces /// Non-owning reference to the list of subfaces
@ -81,15 +80,10 @@ double newton_method(const std::function<implicit_equation_intermediate(double)>
double tolerance = 1e-8, double tolerance = 1e-8,
int max_iterations = 100); 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)
double area_integrand(double u, double v,
const Eigen::Vector3d& p,
const Eigen::Vector3d& dU,
const Eigen::Vector3d& dV)
{ {
return 1.0; return 1.0;
} }
} // namespace internal } // namespace internal

140
surface_integral/src/subface_integrator.cpp
View File

@ -8,30 +8,30 @@
namespace internal namespace internal
{ {
double integrator_t::calculate(int gauss_order, double (*func)(double u, double v, double integrator_t::calculate(
const Eigen::Vector3d& p, int gauss_order,
const Eigen::Vector3d& dU, double (*func)(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)) const
const Eigen::Vector3d& dV)) const
{ {
double total_integral = 0.0; double total_integral = 0.0;
for (const auto& [subface_index, param_plane] : m_uv_planes) { for (const auto& [subface_index, param_plane] : m_uv_planes) {
const auto& subface = m_subfaces[subface_index].object_ptr.get(); const auto& subface = m_subfaces[subface_index].object_ptr.get();
total_integral += calculate_one_subface(subface, param_plane, gauss_order, func); total_integral += calculate_one_subface(subface, param_plane, gauss_order, func);
} }
return total_integral; return total_integral;
} }
double integrator_t::calculate_one_subface(
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 subface& subface,
const Eigen::Vector3d& p, const parametric_plane_t& param_plane,
const Eigen::Vector3d& dU, int gauss_order,
const Eigen::Vector3d& dV)) const 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(); auto solver = subface.fetch_solver_evaluator();
// Gaussian integration in u direction // Gaussian integration in u direction
auto u_integrand = [&](double u) { auto u_integrand = [&](double u) {
// Find exact v intersections for each u // Find exact v intersections for each u
stl_vector_mp<double> v_breaks = find_v_intersections_at_u(subface, param_plane, u);; stl_vector_mp<double> v_breaks = find_v_intersections_at_u(subface, param_plane, u);
;
// Gaussian integration in v direction // Gaussian integration in v direction
auto v_integrand = [&](double v) { auto v_integrand = [&](double v) {
@ -39,25 +39,25 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet
if (is_point_inside_domain(subface, param_plane, u, v)) { if (is_point_inside_domain(subface, param_plane, u, v)) {
try { try {
// Get evaluators for both directions // Get evaluators for both directions
auto eval_du = subface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // Fix v, vary u → ∂/∂u auto eval_du = subface.fetch_curve_constraint_evaluator(parameter_v_t{}, v); // Fix v, vary u → ∂/∂u
auto eval_dv = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // Fix u, vary v → ∂/∂v auto eval_dv = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u); // Fix u, vary v → ∂/∂v
auto res_u = eval_du(u); // f(u,v), grad_f = ∂r/∂u auto res_u = eval_du(u); // f(u,v), grad_f = ∂r/∂u
auto res_v = eval_dv(v); // f(u,v), grad_f = ∂r/∂v auto res_v = eval_dv(v); // f(u,v), grad_f = ∂r/∂v
Eigen::Vector3d p = res_u.f.template head<3>(); // Point (x,y,z) Eigen::Vector3d p = res_u.f.template head<3>(); // Point (x,y,z)
Eigen::Vector3d dU = res_u.grad_f.template head<3>(); // ∂r/∂u Eigen::Vector3d dU = res_u.grad_f.template head<3>(); // ∂r/∂u
Eigen::Vector3d dV = res_v.grad_f.template head<3>(); // ∂r/∂v Eigen::Vector3d dV = res_v.grad_f.template head<3>(); // ∂r/∂v
// Area element: ||dU × dV|| // Area element: ||dU × dV||
Eigen::Vector3d cross = dU.cross(dV); Eigen::Vector3d cross = dU.cross(dV);
double jacobian = cross.norm(); // Jacobian (area scaling factor) double jacobian = cross.norm(); // Jacobian (area scaling factor)
return func(u, v, p, dU, dV) * jacobian; return func(u, v, p, dU, dV) * jacobian;
} catch (...) { } catch (...) {
return 0.0; // Skip singular points return 0.0; // Skip singular points
} }
} }
return 0.0; // Point not in domain return 0.0; // Point not in domain
}; };
double v_integral = 0.0; double v_integral = 0.0;
@ -78,16 +78,16 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet
}; };
// Integrate in u direction // 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, 0.0, 1.0);
double integral = 0.0; double integral = 0.0;
for (size_t i = 0; i < u_breaks.size() - 1; ++i) { for (size_t i = 0; i < u_breaks.size() - 1; ++i) {
double a = u_breaks[i]; double a = u_breaks[i];
double b = u_breaks[i + 1]; double b = u_breaks[i + 1];
double mid_u = (a + b) / 2.0; 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);
if (!v_intersections.empty()) { if (!v_intersections.empty()) {
integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u)); integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u));
} }
} }
@ -98,15 +98,12 @@ double integrator_t::calculate_one_subface(const subface& subface, const paramet
/** /**
* @brief Compute the u-parameter breakpoints for integration. * @brief Compute the u-parameter breakpoints for integration.
* @note The function currently uses std::set for uniqueness and sorting, * @note The function currently uses std::set for uniqueness and sorting,
* but floating-point precision may cause near-duplicate values to be * but floating-point precision may cause near-duplicate values to be
* treated as distinct. A tolerance-based comparison is recommended. * treated as distinct. A tolerance-based comparison is recommended.
* TODO: Use a tolerance-based approach to avoid floating-point precision issues * TODO: Use a tolerance-based approach to avoid floating-point precision issues
* when inserting u-values (e.g., merge values within 1e-12). * when inserting u-values (e.g., merge values within 1e-12).
*/ */
stl_vector_mp<double> integrator_t::compute_u_breaks( stl_vector_mp<double> integrator_t::compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max) const
const parametric_plane_t& param_plane,
double u_min,
double u_max) const
{ {
std::set<double> break_set; std::set<double> break_set;
@ -119,41 +116,39 @@ stl_vector_mp<double> integrator_t::compute_u_breaks(
const auto& vertices = param_plane.chain_vertices[i]; const auto& vertices = param_plane.chain_vertices[i];
auto& vertex_flags = param_plane.vertex_special_flags[i]; auto& vertex_flags = param_plane.vertex_special_flags[i];
for (size_t j = 0; j < vertices.size(); ++j) { 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]) { break_set.insert(vertices[j].x()); } // Special vertex u
} }
} }
// Return as vector (sorted and unique due to set) // Return as vector (sorted and unique due to set)
return stl_vector_mp<double>(break_set.begin(), break_set.end()); return stl_vector_mp<double>(break_set.begin(), break_set.end());
} }
stl_vector_mp<double> integrator_t::find_v_intersections_at_u(const subface& subface,
stl_vector_mp<double> integrator_t::find_v_intersections_at_u( const parametric_plane_t& param_plane,
const subface& subface, double u_val) const
const parametric_plane_t& param_plane,
double u_val) const
{ {
stl_vector_mp<double> intersections; stl_vector_mp<double> intersections;
// Iterate over each boundary chain // Iterate over each boundary chain
for (const auto& chain : param_plane.chain_vertices) { for (const auto& chain : param_plane.chain_vertices) {
const size_t n_vertices = chain.size(); const size_t n_vertices = chain.size();
if (n_vertices < 2) continue; // Skip degenerate chains if (n_vertices < 2) continue; // Skip degenerate chains
// Iterate over each edge in the chain (including closing edge if desired) // Iterate over each edge in the chain (including closing edge if desired)
for (size_t i = 0; i < n_vertices; ++i) { for (size_t i = 0; i < n_vertices; ++i) {
size_t j = (i + 1) % n_vertices; // Next vertex index (wraps around for closed chain) 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& v1 = chain[i]; // Current vertex: (u1, v1)
const Eigen::Vector2d& v2 = chain[j]; // Next vertex: (u2, v2) const Eigen::Vector2d& v2 = chain[j]; // Next vertex: (u2, v2)
double u1 = v1.x(), v1_val = v1.y(); double u1 = v1.x(), v1_val = v1.y();
double u2 = v2.x(), v2_val = v2.y(); double u2 = v2.x(), v2_val = v2.y();
// Classify position relative to u = u_val: -1 (left), 0 (on), +1 (right) // Classify position relative to u = u_val: -1 (left), 0 (on), +1 (right)
const double eps = 1e-12; const double eps = 1e-12;
auto sign_cmp = [u_val, eps](double u) -> int { auto sign_cmp = [u_val, eps](double u) -> int {
if (u < u_val - eps) return -1; if (u < u_val - eps) return -1;
if (u > u_val + eps) return 1; if (u > u_val + eps) return 1;
return 0; return 0;
}; };
@ -162,9 +157,7 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(
int pos2 = sign_cmp(u2); int pos2 = sign_cmp(u2);
// Case 1: Both endpoints on the same side (and not on the line) → no intersection // Case 1: Both endpoints on the same side (and not on the line) → no intersection
if (pos1 * pos2 > 0) { if (pos1 * pos2 > 0) { continue; }
continue;
}
// Case 2: Both endpoints lie exactly on u = u_val → add both v-values // Case 2: Both endpoints lie exactly on u = u_val → add both v-values
if (pos1 == 0 && pos2 == 0) { if (pos1 == 0 && pos2 == 0) {
@ -178,12 +171,11 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(
intersections.push_back(v1_val); intersections.push_back(v1_val);
intersections.push_back(v2_val); intersections.push_back(v2_val);
} }
} } else {
else {
// General case: non-vertical segment crossing u = u_val // General case: non-vertical segment crossing u = u_val
// Compute linear interpolation parameter t // Compute linear interpolation parameter t
double t = (u_val - u1) / (u2 - u1); double t = (u_val - u1) / (u2 - u1);
double v_initial = v1_val + t * (v2_val - v1_val); // Initial guess for v double v_initial = v1_val + t * (v2_val - v1_val); // Initial guess for v
// Fetch evaluators from subface for constraint solving // Fetch evaluators from subface for constraint solving
auto curve_evaluator = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u_val); auto curve_evaluator = subface.fetch_curve_constraint_evaluator(parameter_u_t{}, u_val);
@ -192,7 +184,7 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(
// Define target function: v ↦ residual of implicit equation // Define target function: v ↦ residual of implicit equation
auto target_function = [&](double v) -> internal::implicit_equation_intermediate { auto target_function = [&](double v) -> internal::implicit_equation_intermediate {
constraint_curve_intermediate temp_res = curve_evaluator(v); constraint_curve_intermediate temp_res = curve_evaluator(v);
auto full_res = solver_evaluator(std::move(temp_res)); auto full_res = solver_evaluator(std::move(temp_res));
return std::get<internal::implicit_equation_intermediate>(full_res); return std::get<internal::implicit_equation_intermediate>(full_res);
}; };
@ -209,9 +201,7 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(
} }
// Final step: sort and remove duplicates within tolerance // Final step: sort and remove duplicates within tolerance
if (!intersections.empty()) { if (!intersections.empty()) { sort_and_unique_with_tol(intersections, 1e-8); }
sort_and_unique_with_tol(intersections, 1e-8);
}
return intersections; return intersections;
} }
@ -223,33 +213,28 @@ stl_vector_mp<double> integrator_t::find_v_intersections_at_u(
NOTE: when v_intersect - v < threshold, further checks are required to accurately determine if the intersection point lies NOTE: when v_intersect - v < threshold, further checks are required to accurately determine if the intersection point lies
precisely on the boundary segment. precisely on the boundary segment.
*/ */
bool integrator_t::is_point_inside_domain( bool integrator_t::is_point_inside_domain(const subface& subface,
const subface& subface, const parametric_plane_t& param_plane,
const parametric_plane_t& param_plane, double u,
double u, double v) const
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);
const double tol_near = 1e-8; const double tol_near = 1e-8;
const double tol_above = 1e-12; const double tol_above = 1e-12;
uint32_t count = 0; uint32_t count = 0;
for (double v_int : intersections) { for (double v_int : intersections) {
double diff = v_int - v; double diff = v_int - v;
if (diff > tol_above) { if (diff > tol_above) {
count++; count++;
} } else if (std::abs(diff) < tol_near) {
else if (std::abs(diff) < tol_near) {
return true; // on boundary → inside return true; // on boundary → inside
} }
} }
return (count % 2) == 1; return (count % 2) == 1;
} }
bool integrator_t::is_u_near_singularity(double u, double tol) const bool integrator_t::is_u_near_singularity(double u, double tol) const { return false; }
{
return false;
}
void integrator_t::sort_and_unique_with_tol(stl_vector_mp<double>& vec, double epsilon) const void integrator_t::sort_and_unique_with_tol(stl_vector_mp<double>& vec, double epsilon) const
{ {
@ -269,17 +254,16 @@ void integrator_t::sort_and_unique_with_tol(stl_vector_mp<double>& vec, double e
} }
// Only accepts functions that return implicit_equation_intermediate // Only accepts functions that return implicit_equation_intermediate
double newton_method( double newton_method(const std::function<internal::implicit_equation_intermediate(double)>& F,
const std::function<internal::implicit_equation_intermediate(double)>& F, double v_initial,
double v_initial, double tolerance,
double tolerance, int max_iterations)
int max_iterations)
{ {
double v = v_initial; double v = v_initial;
for (int i = 0; i < max_iterations; ++i) { for (int i = 0; i < max_iterations; ++i) {
auto res = F(v); // Known type: implicit_equation_intermediate auto res = F(v); // Known type: implicit_equation_intermediate
double f_val = res.f; double f_val = res.f;
double df_val = res.df; double df_val = res.df;
if (std::abs(f_val) < tolerance) { if (std::abs(f_val) < tolerance) {

Write Preview
Loading…
Cancel
Save