diff --git a/.gitignore b/.gitignore index f337d76b8..fe03e6d0d 100644 --- a/.gitignore +++ b/.gitignore @@ -672,4 +672,6 @@ python/update_bindings.py tests/data notebooks/*.png CMakeGraphVizOptions.cmake -.zed/* \ No newline at end of file +.zed/* + +CLAUDE.md \ No newline at end of file diff --git a/src/ipc/distance/distance_type.cpp b/src/ipc/distance/distance_type.cpp index 000efc348..aa9b145b9 100644 --- a/src/ipc/distance/distance_type.cpp +++ b/src/ipc/distance/distance_type.cpp @@ -85,7 +85,12 @@ EdgeEdgeDistanceType edge_edge_distance_type( Eigen::ConstRef eb0, Eigen::ConstRef eb1) { - constexpr double PARALLEL_THRESHOLD = 1.0e-20; + // Relative sin² threshold for parallelism: treat edges as parallel when + // sin²(θ) < PARALLEL_THRESHOLD, scaled by a*c. This avoids misclassifying + // short nearly-collinear coplanar segments (which the old absolute + // threshold could not handle) and routes such cases to + // edge_edge_parallel_distance_type. + constexpr double PARALLEL_THRESHOLD = 2.5e-16; const Eigen::Vector3d u = ea1 - ea0; const Eigen::Vector3d v = eb1 - eb0; @@ -108,7 +113,7 @@ EdgeEdgeDistanceType edge_edge_distance_type( } // Special handling for parallel edges - const double parallel_tolerance = PARALLEL_THRESHOLD * std::max(1.0, a * c); + const double parallel_tolerance = PARALLEL_THRESHOLD * a * c; if (u.cross(v).squaredNorm() < parallel_tolerance) { return edge_edge_parallel_distance_type(ea0, ea1, eb0, eb1); } diff --git a/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp index 5b6c171c8..d079c0a9a 100644 --- a/src/ipc/potentials/normal_potential.cpp +++ b/src/ipc/potentials/normal_potential.cpp @@ -148,8 +148,16 @@ VectorMax12d NormalPotential::gradient( const NormalCollision& collision, Eigen::ConstRef positions) const { + // When m=0 the mollifier also satisfies ∇m=0, so the full gradient + // f(d)·∇m + m·f'(d)·∇d = 0 without needing to evaluate f'(d). + const double m = collision.mollifier(positions); // m(x) + if (m <= 0) { + return VectorMax12d::Zero(positions.size()); + } + // d(x) const double d = collision.compute_distance(positions); + // ∇d(x) const VectorMax12d grad_d = collision.compute_distance_gradient(positions); @@ -163,7 +171,7 @@ VectorMax12d NormalPotential::gradient( return (collision.weight * grad_f) * grad_d; } - const double m = collision.mollifier(positions); // m(x) + // Mollified (edge-edge) path: need to apply product rule to m(x) * f(d(x)). const VectorMax12d grad_m = collision.mollifier_gradient(positions); // ∇m(x) @@ -179,6 +187,16 @@ MatrixMax12d NormalPotential::hessian( { // d(x) const double d = collision.compute_distance(positions); + + // Mollified (edge-edge) path: same reasoning as gradient() — check + // m(x) before evaluating barrier derivatives. + const double m = collision.mollifier(positions); // m(x) + if (collision.is_mollified() && m <= 0) { + const double f = (*this)(d, collision.dmin); + const MatrixMax12d hess_m = collision.mollifier_hessian(positions); + return (collision.weight * f) * hess_m; + } + // ∇d(x) const VectorMax12d grad_d = collision.compute_distance_gradient(positions); // ∇²d(x) @@ -198,8 +216,6 @@ MatrixMax12d NormalPotential::hessian( } else { const double f = (*this)(d, collision.dmin); // f(d(x)) - // m(x) - const double m = collision.mollifier(positions); // ∇ m(x) const VectorMax12d grad_m = collision.mollifier_gradient(positions); // ∇² m(x) diff --git a/tests/src/tests/distance/test_distance_type.cpp b/tests/src/tests/distance/test_distance_type.cpp index e5b76de80..23f993195 100644 --- a/tests/src/tests/distance/test_distance_type.cpp +++ b/tests/src/tests/distance/test_distance_type.cpp @@ -6,6 +6,7 @@ #include #include +#include #include using namespace ipc; @@ -421,3 +422,90 @@ TEST_CASE( || dtype == EdgeEdgeDistanceType::EA1_EB)); } } + +// Regression test for a bug where two coplanar nearly-collinear non-overlapping +// edges were assigned EdgeEdgeDistanceType::EA_EB by edge_edge_distance_type(). +// EA_EB delegates to line_line_distance(), whose formula +// d² = ((eb0−ea0)·(u×v))² / |u×v|² +// returns exactly 0 for coplanar edges (the vector eb0−ea0 lies in the same +// plane as u and v, so it is perpendicular to the cross product). This caused +// the IPC barrier gradient to assert-fail on isfinite(f'(0)). +// +// The correct type for non-overlapping parallel/collinear segments is an +// endpoint–endpoint type (here EA1_EB0), which gives a positive distance equal +// to point_point_distance(ea1, eb0). +TEST_CASE( + "Edge-edge coplanar near-collinear non-overlapping distance type", + "[distance][distance-type][edge-edge][regression]") +{ + // Exact vertex coordinates captured from the cloth-ball crash. + // All four vertices share the same x-coordinate, making the two edges + // coplanar in the plane x = -0.81818…. Their directions are nearly + // identical (differ only by ~3 nm in y), so |u×v|² ≈ 1.63e-19. In the + // current edge_edge_distance_type() implementation, parallelism is tested + // against the relative threshold PARALLEL_THRESHOLD * a * c with + // PARALLEL_THRESHOLD = 2.5e-16; here a*c ≈ 0.001, so the effective + // threshold is ≈ 2.7e-19, which is slightly larger than |u×v|². + SECTION("exact crash coordinates") + { + const Eigen::Vector3d ea0( + -0.81818181276321411, 0.073941159961546266, 0.090909108519554152); + const Eigen::Vector3d ea1( + -0.81818181276321411, 0.073941161500775773, 0.272727280855178830); + const Eigen::Vector3d eb0( + -0.81818181276321411, 0.073941163540152718, 0.454545468091964780); + const Eigen::Vector3d eb1( + -0.81818181276321411, 0.073941167300585323, 0.636363625526428220); + + const EdgeEdgeDistanceType dtype = + edge_edge_distance_type(ea0, ea1, eb0, eb1); + + // Must not be EA_EB — that would give line_line_distance = 0. + CHECK(dtype != EdgeEdgeDistanceType::EA_EB); + + // The correct type for non-overlapping collinear segments separated + // along their shared direction is EA1_EB0 (closest endpoints facing + // each other). + CHECK(dtype == EdgeEdgeDistanceType::EA1_EB0); + + // Distance must be strictly positive. + const double dist = edge_edge_distance(ea0, ea1, eb0, eb1, dtype); + CHECK(dist > 0); + CHECK(std::isfinite(dist)); + } + + // Parametric version: build two parallel segments in the same plane + // (x = 0) with a variable z-gap between them, separated by a small + // in-plane perpendicular offset dy. The closest features are always the + // facing endpoints ea1 and eb0. + SECTION("parametric coplanar parallel non-overlapping") + { + // Length of each segment (same for both) + const double L = GENERATE(0.1, 0.5, 1.0, 2.0); + // Gap between ea1.z and eb0.z + const double gap = GENERATE(0.05, 0.1, 0.2, 0.5); + // Tiny in-plane y-offset (mimics the ~3 nm offset in the crash case) + const double dy = GENERATE(0.0, 1e-9, 1e-6); + + const Eigen::Vector3d ea0(0.0, 0.0, 0.0); + const Eigen::Vector3d ea1(0.0, dy, L); + const Eigen::Vector3d eb0(0.0, 2.0 * dy, L + gap); + const Eigen::Vector3d eb1(0.0, 3.0 * dy, L + gap + L); + + CAPTURE(L, gap, dy); + + const EdgeEdgeDistanceType dtype = + edge_edge_distance_type(ea0, ea1, eb0, eb1); + + CHECK(dtype != EdgeEdgeDistanceType::EA_EB); + + const double dist = edge_edge_distance(ea0, ea1, eb0, eb1, dtype); + CHECK(dist > 0); + CHECK(std::isfinite(dist)); + + // The distance must equal point_point_distance(ea1, eb0) since those + // are the closest points. + const double expected = (ea1 - eb0).squaredNorm(); + CHECK(dist == Catch::Approx(expected).margin(1e-10)); + } +} \ No newline at end of file diff --git a/tests/src/tests/distance/test_edge_edge.cpp b/tests/src/tests/distance/test_edge_edge.cpp index 0aad5816f..a3dadcb3e 100644 --- a/tests/src/tests/distance/test_edge_edge.cpp +++ b/tests/src/tests/distance/test_edge_edge.cpp @@ -67,8 +67,9 @@ TEST_CASE("Edge-edge distance !EA_EB", "[distance][edge-edge]") { double alpha = GENERATE(range(-1.0, 2.0, 0.1)); double s = GENERATE(range(-10.0, 10.0, 1.0)); - if (s == 0) + if (s == 0) { return; + } const bool swap_ea = GENERATE(false, true); const bool swap_eb = GENERATE(false, true); const bool swap_edges = GENERATE(false, true); @@ -342,8 +343,9 @@ TEST_CASE( { double alpha = GENERATE(range(-1.0, 2.0, 0.1)); double s = GENERATE(range(-10.0, 10.0, 1.0)); - if (s == 0) + if (s == 0) { return; + } const bool swap_ea = GENERATE(false, true); const bool swap_eb = GENERATE(false, true); const bool swap_edges = GENERATE(false, true); @@ -850,8 +852,9 @@ TEMPLATE_TEST_CASE_SIG( { double alpha = GENERATE(range(-1.0, 2.0, 0.4)); double s = GENERATE(range(-10.0, 10.0, 5.0)); - if (s == 0) + if (s == 0) { return; + } const bool swap_ea = GENERATE(false, true); const bool swap_eb = GENERATE(false, true); const bool swap_edges = GENERATE(false, true);