Fix edge-edge parallelism threshold and add regression test

.gitignore

@@ -672,4 +672,6 @@ python/update_bindings.py
 tests/data
 notebooks/*.png
 CMakeGraphVizOptions.cmake
-.zed/*
+.zed/*
+
+CLAUDE.md

src/ipc/distance/distance_type.cpp

@@ -85,7 +85,12 @@ EdgeEdgeDistanceType edge_edge_distance_type(
     Eigen::ConstRef<Eigen::Vector3d> eb0,
     Eigen::ConstRef<Eigen::Vector3d> 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);
     }

src/ipc/potentials/normal_potential.cpp

@@ -148,8 +148,16 @@ VectorMax12d NormalPotential::gradient(
     const NormalCollision& collision,
     Eigen::ConstRef<VectorMax12d> 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)

tests/src/tests/distance/test_distance_type.cpp

@@ -6,6 +6,7 @@
 #include <catch2/generators/catch_generators_random.hpp>
 
 #include <ipc/distance/distance_type.hpp>
+#include <ipc/distance/edge_edge.hpp>
 #include <ipc/distance/point_triangle.hpp>
 
 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));
+    }
+}

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);