Fix edge_edge_mollifier_gradient_wrt_x to use ∇_rest ε_x

src/ipc/distance/edge_edge_mollifier.cpp

@@ -39,11 +39,20 @@ Matrix12d edge_edge_cross_squarednorm_hessian(
     return hess;
 }
 
+// The active branch returns M(x, eps_x) := m(x, eps_x)^2, where the original
+// base mollifier was m(x, eps_x) = 2y - y^2 (y = x/eps_x). Squaring makes the
+// blend C^1 at the parallel limit (M'(0)=0) — the original m had m'(0)=2/eps_x,
+// which combined with the C^0 EE distance at parallel edges produced a
+// spurious inflation in the EE shape derivative (see repro tests under
+// tests/src/tests/potential/ipc_test_ee_near_parallel_shape_derivative.cpp).
+// The five scalar derivatives below are the analytic chain-rule of M = m^2
+// in the partial m and its scalar partials.
 double edge_edge_mollifier(const double x, const double eps_x)
 {
     if (x < eps_x) {
-        const double x_div_eps_x = x / eps_x;
-        return (-x_div_eps_x + 2) * x_div_eps_x;
+        const double y = x / eps_x;
+        const double m = (2.0 - y) * y; // 2y - y^2
+        return m * m;
     } else {
         return 1;
     }
@@ -52,8 +61,11 @@ double edge_edge_mollifier(const double x, const double eps_x)
 double edge_edge_mollifier_gradient(const double x, const double eps_x)
 {
     if (x < eps_x) {
-        const double one_div_eps_x = 1 / eps_x;
-        return 2 * one_div_eps_x * fma(-one_div_eps_x, x, 1);
+        const double y = x / eps_x;
+        const double m = (2.0 - y) * y;
+        const double m_x = 2.0 * (1.0 - y) / eps_x; // dm/dx
+        // dM/dx = 2 m * dm/dx
+        return 2.0 * m * m_x;
     } else {
         return 0;
     }
@@ -62,7 +74,12 @@ double edge_edge_mollifier_gradient(const double x, const double eps_x)
 double edge_edge_mollifier_hessian(const double x, const double eps_x)
 {
     if (x < eps_x) {
-        return -2 / (eps_x * eps_x);
+        const double y = x / eps_x;
+        const double m = (2.0 - y) * y;
+        const double m_x = 2.0 * (1.0 - y) / eps_x; // dm/dx
+        const double m_xx = -2.0 / (eps_x * eps_x); // d2m/dx2
+        // d2M/dx2 = 2 (dm/dx)^2 + 2 m * d2m/dx2
+        return 2.0 * m_x * m_x + 2.0 * m * m_xx;
     } else {
         return 0;
     }
@@ -71,13 +88,32 @@ double edge_edge_mollifier_hessian(const double x, const double eps_x)
 double
 edge_edge_mollifier_derivative_wrt_eps_x(const double x, const double eps_x)
 {
-    return x < eps_x ? (2 * x * (-eps_x + x) / (eps_x * eps_x * eps_x)) : 0.0;
+    if (x < eps_x) {
+        const double y = x / eps_x;
+        const double m = (2.0 - y) * y;
+        const double m_eps =
+            2.0 * x * (x - eps_x) / (eps_x * eps_x * eps_x); // dm/deps
+        // dM/deps = 2 m * dm/deps
+        return 2.0 * m * m_eps;
+    }
+    return 0.0;
 }
 
 double edge_edge_mollifier_gradient_derivative_wrt_eps_x(
     const double x, const double eps_x)
 {
-    return x < eps_x ? (2 * (-eps_x + 2 * x) / (eps_x * eps_x * eps_x)) : 0.0;
+    if (x < eps_x) {
+        const double y = x / eps_x;
+        const double m = (2.0 - y) * y;
+        const double m_x = 2.0 * (1.0 - y) / eps_x;
+        const double m_eps =
+            2.0 * x * (x - eps_x) / (eps_x * eps_x * eps_x);
+        const double m_eps_x =
+            2.0 * (-eps_x + 2.0 * x) / (eps_x * eps_x * eps_x); // d2m/deps dx
+        // d2M/deps dx = 2 (dm/deps)(dm/dx) + 2 m (d2m/deps dx)
+        return 2.0 * m_eps * m_x + 2.0 * m * m_eps_x;
+    }
+    return 0.0;
 }
 
 double edge_edge_mollifier(
@@ -150,10 +186,13 @@ Vector12d edge_edge_mollifier_gradient_wrt_x(
     const double ee_cross_norm_sqr =
         edge_edge_cross_squarednorm(ea0, ea1, eb0, eb1);
     if (ee_cross_norm_sqr < eps_x) {
-        // ∇ₓ m = ∂m/∂ε ∇ₓε
+        // ∇ₓ m = ∂m/∂ε · ∇ₓε
+        // (m depends on rest positions only through eps_x, since the
+        // cross-squarednorm s is a function of POSITIONS only)
         return edge_edge_mollifier_derivative_wrt_eps_x(
                    ee_cross_norm_sqr, eps_x)
-            * edge_edge_mollifier_gradient(ea0, ea1, eb0, eb1, eps_x);
+            * edge_edge_mollifier_threshold_gradient(
+                   ea0_rest, ea1_rest, eb0_rest, eb1_rest);
     } else {
         return Vector12d::Zero();
     }

tests/src/tests/potential/CMakeLists.txt

@@ -5,6 +5,7 @@ set(SOURCES
   test_smooth_potential.cpp
   test_friction_potential.cpp
   test_distance_vector_methods.cpp
+  ipc_test_ee_near_parallel_shape_derivative.cpp
 
   # Benchmarks
 

tests/src/tests/potential/ipc_test_ee_near_parallel_shape_derivative.cpp

@@ -0,0 +1,149 @@
+// Regression test for the edge-edge parallel mollifier shape-derivative path.
+//
+// Reproduces a single mollified EE pair (near-parallel edges, distance near
+// dhat) and FD-checks the analytic shape_derivative against a centered
+// difference of `gradient` at fixed collision set — the same probe IPC's own
+// shape_derivative test uses (test_barrier_potential.cpp:377-391).
+//
+// Prior to the fix in edge_edge_mollifier.cpp, the assertion failed by
+// roughly two orders of magnitude on these geometries: the analytic
+// per-pair shape_derivative was inflated when (a) the EE pair is mollified
+// (cross_squared_norm < eps_x), AND (b) the closest distance is near dhat.
+// The root cause was `edge_edge_mollifier_gradient_wrt_x` returning
+// ∂m/∂eps_x · ∇_position m instead of ∂m/∂eps_x · ∇_rest eps_x. The fix
+// uses edge_edge_mollifier_threshold_gradient as the second factor, matching
+// what edge_edge_mollifier_gradient_jacobian_wrt_x already does.
+
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/catch_approx.hpp>
+#include <catch2/generators/catch_generators.hpp>
+
+#include <ipc/collision_mesh.hpp>
+#include <ipc/collisions/normal/normal_collisions.hpp>
+#include <ipc/potentials/barrier_potential.hpp>
+
+#include <Eigen/Core>
+
+using namespace ipc;
+
+namespace {
+
+// Build a 2-edge / 4-vertex CollisionMesh: two unit-length parallel-ish edges
+// in 3D, separated in y by `y_sep`, with a small z-twist `z_twist` on the
+// second edge to keep the cross product slightly nonzero (so the mollifier is
+// triggered but s = |A x B|^2 != 0).
+//
+//   ea0 = (0, 0,        0)        ea1 = (1, 0,        0)
+//   eb0 = (0, y_sep,    0)        eb1 = (1, y_sep,    z_twist)
+//
+// With z_twist > 0 small, |A x B|^2 = z_twist^2 > 0 but << eps_x = 1e-3
+// (since both edges have length ~1), so the mollifier branch fires. The
+// closest-point distance between the edges stays near y_sep.
+struct EEScene {
+    Eigen::MatrixXd rest;
+    Eigen::MatrixXd displaced;
+    Eigen::MatrixXi edges;
+    Eigen::MatrixXi faces;
+};
+
+EEScene make_near_parallel_ee_scene(double y_sep, double z_twist)
+{
+    EEScene s;
+    s.rest.resize(4, 3);
+    s.rest << 0.0, 0.0,     0.0,
+              1.0, 0.0,     0.0,
+              0.0, y_sep,   0.0,
+              1.0, y_sep,   z_twist;
+    // displaced = rest: the shape-derivative bug is independent of any
+    // additional displacement, so we test the simplest configuration.
+    s.displaced = s.rest;
+    s.edges.resize(2, 2);
+    s.edges << 0, 1,
+               2, 3;
+    s.faces.resize(0, 3);
+    return s;
+}
+
+} // namespace
+
+TEST_CASE(
+    "EE shape_derivative on near-parallel mollified pair",
+    "[potential][barrier_potential][shape_derivative][mollifier]")
+{
+    const double dhat = 1e-3;
+
+    // y_sep just below dhat (active contact), z_twist makes |A x B|^2 small
+    // but nonzero so the mollifier branch is entered cleanly.
+    const double y_sep   = GENERATE(9.99e-4, 9.9e-4);
+    const double z_twist = GENERATE(1e-3);
+
+    EEScene scene = make_near_parallel_ee_scene(y_sep, z_twist);
+    CollisionMesh mesh(scene.rest, scene.edges, scene.faces);
+    if (!mesh.are_area_jacobians_initialized())
+        mesh.init_area_jacobians();
+
+    NormalCollisions collisions;
+    collisions.set_enable_shape_derivatives(true);
+    collisions.build(mesh, scene.displaced, dhat);
+    REQUIRE(collisions.size() == 1);
+    REQUIRE(collisions.is_edge_edge(0));
+    REQUIRE(collisions[0].is_mollified());
+
+    // kappa = 1e8 keeps entries well above the precision floor of the
+    // centered-FD probe.
+    BarrierPotential bp(dhat, /*stiffness=*/1e8, /*use_physical_barrier=*/false);
+
+    // -------- Analytic 12x12 (collision-mesh DOF order) --------
+    Eigen::SparseMatrix<double> JF_an_sparse =
+        bp.shape_derivative(collisions, mesh, scene.displaced);
+    Eigen::MatrixXd JF_an(JF_an_sparse);
+
+    // -------- FD 12x12 column-by-column on rest DOFs, material
+    // displacement held constant, collision set held fixed. Build in
+    // ROW-MAJOR DOF order (v * dim + d) to match the layout returned by
+    // shape_derivative. --------
+    const int n_verts = static_cast<int>(scene.rest.rows());
+    const int dim     = static_cast<int>(scene.rest.cols());
+    const int ndof    = n_verts * dim;
+    const double eps  = 1e-7;
+    Eigen::MatrixXd JF_fd(ndof, ndof);
+    JF_fd.setZero();
+    for (int v = 0; v < n_verts; ++v) {
+        for (int d = 0; d < dim; ++d) {
+            Eigen::MatrixXd rest_plus  = scene.rest;
+            Eigen::MatrixXd rest_minus = scene.rest;
+            rest_plus(v, d)  += eps;
+            rest_minus(v, d) -= eps;
+            // Hold material displacement (= displaced - rest) constant.
+            const Eigen::MatrixXd disp_plus =
+                rest_plus  + (scene.displaced - scene.rest);
+            const Eigen::MatrixXd disp_minus =
+                rest_minus + (scene.displaced - scene.rest);
+            CollisionMesh mesh_plus(rest_plus, scene.edges, scene.faces);
+            if (mesh.are_area_jacobians_initialized())
+                mesh_plus.init_area_jacobians();
+            CollisionMesh mesh_minus(rest_minus, scene.edges, scene.faces);
+            if (mesh.are_area_jacobians_initialized())
+                mesh_minus.init_area_jacobians();
+            const Eigen::VectorXd F_plus =
+                bp.gradient(collisions, mesh_plus,  disp_plus);
+            const Eigen::VectorXd F_minus =
+                bp.gradient(collisions, mesh_minus, disp_minus);
+            JF_fd.col(v * dim + d) = (F_plus - F_minus) / (2.0 * eps);
+        }
+    }
+
+    const double an   = JF_an.norm();
+    const double fn   = JF_fd.norm();
+    const double diff = (JF_an - JF_fd).norm();
+    const double rel  = diff / std::max(std::max(an, fn), 1e-30);
+
+    UNSCOPED_INFO(
+        "y_sep=" << y_sep << " z_twist=" << z_twist
+        << "  ||analytic||=" << an
+        << "  ||fd||=" << fn
+        << "  ||diff||=" << diff
+        << "  rel=" << rel);
+
+    CHECK(rel < 5e-2);
+}