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

fix: find_v_intersections_at_u, add bound

feat-integrator
mckay 2 weeks ago
parent
c04fba163b
commit
c6ccc29d72
2 changed files with 38 additions and 16 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 5
      surface_integral/interface/subface_integrator.hpp
  2. 49
      surface_integral/src/subface_integrator.cpp

5
surface_integral/interface/subface_integrator.hpp
View File

@ -65,7 +65,10 @@ private:
stl_vector_mp<double> find_v_intersections_at_u(const subface& subface,
const parametric_plane_t& param_plane,
double u_val) const;
double u_val,
bool if_add_boundary=true,
double v_min=0.0,
double v_max=0.0) const;
bool is_point_inside_domain(const subface& subface, const parametric_plane_t& param_plane, double u, double v) const;

49
surface_integral/src/subface_integrator.cpp
View File

@ -13,9 +13,14 @@ double integrator_t::calculate(
double (*func)(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)) const
{
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, func);
for (int subface_index = 0; subface_index < m_subfaces.size(); ++subface_index) {
if (m_uv_planes.find(subface_index) != m_uv_planes.end()) {
const auto& param_plane = m_uv_planes.at(subface_index);
const auto& subface = m_subfaces[subface_index].object_ptr.get();
total_integral += calculate_one_subface(subface, param_plane, gauss_order, func);
} else {
// the subface has no associated parametric plane, meaning it is compete face without trimming
}
}
return total_integral;
}
@ -27,11 +32,15 @@ double integrator_t::calculate_one_subface(
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();
const auto& uv_bounds = param_plane.uv_bounds;
auto u_min = uv_bounds.min().x();
auto u_max = uv_bounds.max().x();
auto v_min = uv_bounds.min().y();
auto v_max = uv_bounds.max().y();
// Gaussian integration in u direction
auto u_integrand = [&](double 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, true,v_min, v_max);
// Gaussian integration in v direction
auto v_integrand = [&](double v) {
@ -78,14 +87,14 @@ double integrator_t::calculate_one_subface(
};
// 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, u_min, u_max);
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_v_intersections_at_u(subface, param_plane, mid_u);
auto v_intersections = find_v_intersections_at_u(subface, param_plane, mid_u,true, v_min, v_max);
if (!v_intersections.empty()) {
integral += integrate_1D(a, b, u_integrand, gauss_order, is_u_near_singularity(mid_u));
@ -115,8 +124,9 @@ stl_vector_mp<double> integrator_t::compute_u_breaks(const parametric_plane_t& p
for (size_t i = 0; i < param_plane.chain_vertices.size(); ++i) {
const auto& vertices = param_plane.chain_vertices[i];
auto& vertex_flags = param_plane.vertex_special_flags[i];
auto& edge_flags = param_plane.edge_near_parallel_flags[i];
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] || edge_flags[j]) { break_set.insert(vertices[j].x()); } // Special vertex u
}
}
// Return as vector (sorted and unique due to set)
@ -125,21 +135,30 @@ stl_vector_mp<double> integrator_t::compute_u_breaks(const parametric_plane_t& p
stl_vector_mp<double> integrator_t::find_v_intersections_at_u(const subface& subface,
const parametric_plane_t& param_plane,
double u_val) const
double u_val,
bool if_add_boundary,
double v_min,
double v_max
) const
{
stl_vector_mp<double> intersections;
if (if_add_boundary) {
intersections.push_back(v_min);
intersections.push_back(v_max);
}
// Iterate over each boundary chain
for (const auto& chain : param_plane.chain_vertices) {
const size_t n_vertices = chain.size();
if (n_vertices < 2) continue; // Skip degenerate chains
const size_t n_vertices = chain.size()-1;
if (n_vertices < 2){
std::cout << "Warning: Degenerate chain with less than 2 vertices." << std::endl;
continue; // Skip degenerate chains
}
// Iterate over each edge in the chain (including closing edge if desired)
for (size_t i = 0; i < n_vertices; ++i) {
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& v2 = chain[j]; // Next vertex: (u2, v2)
const Eigen::Vector2d& v2 = chain[i + 1]; // Next vertex: (u2, v2)
double u1 = v1.x(), v1_val = v1.y();
double u2 = v2.x(), v2_val = v2.y();
@ -218,7 +237,7 @@ bool integrator_t::is_point_inside_domain(const subface& subface,
double u,
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, false);
const double tol_near = 1e-8;
const double tol_above = 1e-12;

Write Preview
Loading…
Cancel
Save