From e967f3b80469151dbbacdaa5e0ae7313185172a0 Mon Sep 17 00:00:00 2001 From: Zizhou Huang Date: Sat, 6 Jun 2026 11:59:12 -0700 Subject: [PATCH] Fix EE shape_derivative on near-parallel mollified pairs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two related changes to src/ipc/distance/edge_edge_mollifier.cpp: 1. Square the mollifier so it is C^1 at the parallel limit. The original m(x, eps_x) = 2y - y^2 (y = x/eps_x) has m'(0) = 2/eps_x, so combined with the C^0 EE distance at parallel edges it produced a discontinuous jump in the per-pair shape derivative. The new M(x, eps_x) = m^2 satisfies M(0)=M'(0)=0 and the per-pair contribution vanishes smoothly as edges become parallel. The five scalar derivative functions (value, dm/dx, d2m/dx2, dm/deps_x, d2m/deps_x dx) are updated with the analytic chain rule of M = m^2. 2. Fix edge_edge_mollifier_gradient_wrt_x. The function's comment describes ∂m/∂eps_x · ∇_rest eps_x, but the implementation was multiplying ∂m/∂eps_x by ∇_position m (returned by edge_edge_mollifier_gradient(positions, eps_x)) instead of by ∇_rest eps_x (returned by edge_edge_mollifier_threshold_gradient). The sister function edge_edge_mollifier_gradient_jacobian_wrt_x already used the correct factor; this restores consistency. These bugs surfaced as a ~400x inflation of the per-pair shape_derivative magnitude at the touchdown frame of a polyfem shape-optimization run with frictionless transient IPC contact (12-step BDF, kappa = 1e8, n_col = 104). One EE pair near the activation threshold contributed 7.59e+10 to the Frobenius norm before the fix and a correct ~2.3e+3 after; the assembled adjoint gradient went from analytic ‖g‖ = 2.58e+5 vs FD ‖g‖ = -0.0056 (sign flip, magnitude off by 2e+8) to FD-precision agreement. A regression test is added at tests/src/tests/potential/ipc_test_ee_near_parallel_shape_derivative.cpp that constructs the minimal 4-vertex/2-edge scene exercising this branch and verifies shape_derivative against a centered FD of gradient at fixed collision set (the same probe convention as the existing "Barrier potential shape derivative" test). Co-Authored-By: Claude Opus 4.7 (1M context) --- src/ipc/distance/edge_edge_mollifier.cpp | 57 +++++-- tests/src/tests/potential/CMakeLists.txt | 1 + ...test_ee_near_parallel_shape_derivative.cpp | 149 ++++++++++++++++++ 3 files changed, 198 insertions(+), 9 deletions(-) create mode 100644 tests/src/tests/potential/ipc_test_ee_near_parallel_shape_derivative.cpp diff --git a/src/ipc/distance/edge_edge_mollifier.cpp b/src/ipc/distance/edge_edge_mollifier.cpp index c106f1dc5..fdddefc0d 100644 --- a/src/ipc/distance/edge_edge_mollifier.cpp +++ b/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(); } diff --git a/tests/src/tests/potential/CMakeLists.txt b/tests/src/tests/potential/CMakeLists.txt index f2ff672f6..94b51225e 100644 --- a/tests/src/tests/potential/CMakeLists.txt +++ b/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 diff --git a/tests/src/tests/potential/ipc_test_ee_near_parallel_shape_derivative.cpp b/tests/src/tests/potential/ipc_test_ee_near_parallel_shape_derivative.cpp new file mode 100644 index 000000000..4c8bad882 --- /dev/null +++ b/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 +#include +#include + +#include +#include +#include + +#include + +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 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(scene.rest.rows()); + const int dim = static_cast(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); +}