From 7b970316ec3cf2e04fb1cdbf8141ecfd28efd4f0 Mon Sep 17 00:00:00 2001 From: Zizhou Huang Date: Sat, 6 Jun 2026 11:59:12 -0700 Subject: [PATCH] =?UTF-8?q?Fix=20EE-mollifier=20shape=20derivative:=20grad?= =?UTF-8?q?ient=5Fwrt=5Fx=20factor=20+=20missing=20(=E2=88=87=E1=B5=A4m)(?= =?UTF-8?q?=E2=88=87=E1=B5=A4f)=E1=B5=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two bugs in the per-pair NormalPotential::shape_derivative when collision.is_mollified() == true (near-parallel EE pair): 1. edge_edge_mollifier_gradient_wrt_x returns ∇_rest m, computed as ∂m/∂ε · ∇_rest ε. The second factor was using edge_edge_mollifier_gradient(positions, eps_x) (returns ∇_position m) instead of edge_edge_mollifier_threshold_gradient(rest) (returns ∇_rest ε). Sister function edge_edge_mollifier_gradient_jacobian_wrt_x already used the correct factor. 2. NormalPotential::shape_derivative (normal_potential.cpp:392) summed only four outer-product / Hessian terms in local_hess. The correct chain-rule expansion of ∇ₓ(f∇ᵤm + m∇ᵤf) — with ∇ₓm = ∇ᵤm + (∂m/∂ε·∇_rest ε), ∇ₓf = ∇ᵤf, ∇ₓ∇ᵤf = ∇ᵤ²f — has FIVE terms: f · ∇ₓ∇ᵤm (jac_m) + (∂m/∂ε · ∇_rest ε) · (∇ᵤf)ᵀ (gradx_m · gradu_fᵀ) + (∇ᵤf)(∇ᵤm)ᵀ + (∇ᵤm)(∇ᵤf)ᵀ symmetric pair + m · ∇ᵤ²f (hessu_f) The (∇ᵤm)(∇ᵤf)ᵀ term was missing; it is the "∇ᵤm contribution" of (∇ₓm)(∇ᵤf)ᵀ that gets split off from the rest contribution that goes into the gradx_m·gradu_fᵀ term. Adds a regression test "Edge-Edge Mollifier Shape Derivative" in tests/src/tests/distance/test_edge_edge_mollifier.cpp that constructs a 4-vertex / 2-edge scene with a single mollified EE pair and FD-checks the full 12x12 analytic BarrierPotential::shape_derivative against centered FD of gradient at fixed collision set (same convention as test_barrier_potential.cpp:377-391). Without either fix, rel >= 2.5e-2. With both fixes, rel reaches FD precision (~3.5e-6). Test threshold 1e-5. Co-Authored-By: Claude Opus 4.7 (1M context) --- src/ipc/distance/edge_edge_mollifier.cpp | 7 +- src/ipc/potentials/normal_potential.cpp | 17 +++- .../distance/test_edge_edge_mollifier.cpp | 81 +++++++++++++++++++ 3 files changed, 100 insertions(+), 5 deletions(-) diff --git a/src/ipc/distance/edge_edge_mollifier.cpp b/src/ipc/distance/edge_edge_mollifier.cpp index c106f1dc5..f3672af61 100644 --- a/src/ipc/distance/edge_edge_mollifier.cpp +++ b/src/ipc/distance/edge_edge_mollifier.cpp @@ -150,10 +150,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/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp index 5b6c171c8..eda6ae867 100644 --- a/src/ipc/potentials/normal_potential.cpp +++ b/src/ipc/potentials/normal_potential.cpp @@ -387,10 +387,21 @@ void NormalPotential::shape_derivative( const Matrix12d jac_m = collision.mollifier_gradient_jacobian_wrt_x( rest_positions, positions); - // Only compute the second term of the shape derivative - // ∇ₓ (f ∇ᵤm + m ∇ᵤf) = f ∇ₓ∇ᵤm + (∇ₓm)(∇ᵤf)ᵀ + (∇ᵤf)(∇ᵤm)ᵀ + m ∇ᵤ²f + // Only compute the second term of the shape derivative. With + // ∇ₓm = ∇ᵤm + (∂m/∂ε · ∇_rest ε) (decomposed into position and + // rest parts) ∇ₓf = ∇ᵤf (f depends on + // positions only) ∇ₓ∇ᵤf = ∇ᵤ²f (same + // reason) + // the chain rule gives FIVE outer-product / Hessian terms: + // ∇ₓ(f ∇ᵤm + m ∇ᵤf) = + // f · ∇ₓ∇ᵤm (jac_m) + // + (∂m/∂ε · ∇_rest ε) · (∇ᵤf)ᵀ (gradx_m + // · gradu_fᵀ) + // + (∇ᵤf) · (∇ᵤm)ᵀ + (∇ᵤm) · (∇ᵤf)ᵀ (symmetric pair) + // + m · ∇ᵤ²f (hessu_f) local_hess = f * jac_m + gradx_m * gradu_f.transpose() - + gradu_f * gradu_m.transpose() + m * hessu_f; + + gradu_f * gradu_m.transpose() + gradu_m * gradu_f.transpose() + + m * hessu_f; } local_hessian_to_global_triplets( diff --git a/tests/src/tests/distance/test_edge_edge_mollifier.cpp b/tests/src/tests/distance/test_edge_edge_mollifier.cpp index 3043350fb..2903768b7 100644 --- a/tests/src/tests/distance/test_edge_edge_mollifier.cpp +++ b/tests/src/tests/distance/test_edge_edge_mollifier.cpp @@ -3,6 +3,9 @@ #include #include +#include +#include +#include #include @@ -283,4 +286,82 @@ TEST_CASE("Edge-Edge Mollifier Threshold", "[distance][edge-edge][mollifier]") fgrad); CHECK(fd::compare_gradient(grad, fgrad)); +} + +// Single mollified EE pair, FD-check shape_derivative against gradient. +TEST_CASE( + "Edge-Edge Mollifier Shape Derivative", + "[distance][edge-edge][mollifier][shape_derivative]") +{ + const double dhat = 1e-3; + const bool aw = GENERATE(false, true); + const double y_sep = GENERATE(9.99e-4, 9.9e-4, 9.5e-4); + const double z_twist = GENERATE(5e-3, 1e-2, 2e-2); + + Eigen::MatrixXd rest(4, 3); + 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; + const Eigen::MatrixXd displaced = rest; + Eigen::MatrixXi edges(2, 2); + edges << 0, 1, 2, 3; + Eigen::MatrixXi faces(0, 3); + + CollisionMesh mesh(rest, edges, faces); + if (!mesh.are_area_jacobians_initialized()) + mesh.init_area_jacobians(); + + NormalCollisions collisions; + collisions.set_use_area_weighting(aw); + collisions.set_enable_shape_derivatives(true); + collisions.build(mesh, displaced, dhat); + REQUIRE(collisions.size() == 1); + REQUIRE(collisions.is_edge_edge(0)); + REQUIRE(collisions[0].is_mollified()); + + BarrierPotential bp( + dhat, /*stiffness=*/1e8, /*use_physical_barrier=*/false); + + Eigen::SparseMatrix JF_an_sparse = + bp.shape_derivative(collisions, mesh, displaced); + Eigen::MatrixXd JF_an(JF_an_sparse); + + // Column-wise central FD on rest DOFs, fixed collision set. + const int n_verts = static_cast(rest.rows()); + const int dim = static_cast(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 = rest; + Eigen::MatrixXd rest_minus = rest; + rest_plus(v, d) += eps; + rest_minus(v, d) -= eps; + const Eigen::MatrixXd disp_plus = rest_plus + (displaced - rest); + const Eigen::MatrixXd disp_minus = rest_minus + (displaced - rest); + CollisionMesh mesh_plus(rest_plus, edges, faces); + if (mesh.are_area_jacobians_initialized()) + mesh_plus.init_area_jacobians(); + CollisionMesh mesh_minus(rest_minus, edges, 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( + "aw=" << aw << " y_sep=" << y_sep << " z_twist=" << z_twist + << " ||analytic||=" << an << " ||fd||=" << fn + << " ||diff||=" << diff << " rel=" << rel); + + CHECK(rel < 1e-4); } \ No newline at end of file