Fix missing PSD projection of mollified m≤0 hessian blocks

src/ipc/potentials/normal_potential.cpp

@@ -8,6 +8,9 @@
 #include <tbb/parallel_for.h>
 #include <tbb/parallel_reduce.h>
 
+#include <algorithm>
+#include <cmath>
+
 namespace ipc {
 
 // -- Cumulative methods -------------------------------------------------------
@@ -183,9 +186,21 @@ MatrixMax12d NormalPotential::hessian(
     // m(x) before evaluating barrier derivatives.
     const double m = collision.mollifier(positions); // m(x)
     if (collision.is_mollified() && m <= 0) {
+        // The mollified hessian is weight * f * hess_m here (the mollifier
+        // gradient is zero when m = 0). At m == 0 the edges are exactly
+        // parallel -- a global minimum of the mollifier -- so hess_m is PSD,
+        // and f = f(d) > 0. The block is therefore a scalar multiple of a PSD
+        // matrix, so its PSD projection reduces to projecting the scalar
+        // weight * f (no eigendecomposition needed).
         const double f = (*this)(d, collision.dmin);
         const MatrixMax12d hess_m = collision.mollifier_hessian(positions);
-        return (collision.weight * f) * hess_m;
+        double scale = collision.weight * f;
+        if (project_hessian_to_psd == PSDProjectionMethod::CLAMP) {
+            scale = std::max(scale, 0.0); // NSD (w < 0) projects to zero
+        } else if (project_hessian_to_psd == PSDProjectionMethod::ABS) {
+            scale = std::abs(scale);
+        }
+        return scale * hess_m;
     }
 
     // ∇d(x)

tests/src/tests/potential/test_barrier_potential.cpp

@@ -15,6 +15,8 @@
 #include <finitediff.hpp>
 #include <igl/edges.h>
 
+#include <Eigen/Eigenvalues>
+
 using namespace ipc;
 
 TEST_CASE(
@@ -466,6 +468,71 @@ TEST_CASE(
     }
 }
 
+TEST_CASE(
+    "Mollified m<=0 hessian block is PSD",
+    "[potential][barrier_potential][hessian]")
+{
+    // The cube's edges are axis-aligned, so IMPROVED_MAX_APPROX's near-parallel
+    // negative-weight collisions have mollifier m == 0 exactly. That exercises
+    // NormalPotential::hessian()'s m <= 0 early-return branch.
+    //
+    // Before the fix, that branch returned (weight * f) * hess_m without
+    // applying the requested PSD projection, so negative-weight collisions
+    // could introduce non-PSD blocks into the assembled "PSD-projected"
+    // hessian.
+    Eigen::MatrixXd vertices;
+    Eigen::MatrixXi edges, faces;
+    REQUIRE(tests::load_mesh("cube.ply", vertices, edges, faces));
+    const CollisionMesh mesh(vertices, edges, faces);
+
+    const double dhat = std::sqrt(2.0);
+    NormalCollisions collisions;
+    collisions.set_use_area_weighting(true);
+    collisions.set_collision_set_type(
+        NormalCollisions::CollisionSetType::IMPROVED_MAX_APPROX);
+    collisions.build(mesh, vertices, dhat);
+    REQUIRE(!collisions.empty());
+
+    // Ensure the test actually exercises NormalPotential::hessian()'s m <= 0
+    // early-return branch (and the negative-weight case that motivated the
+    // fix).
+    int m_le_0_count = 0;
+    int m_le_0_negative_weight_count = 0;
+    for (size_t i = 0; i < collisions.size(); i++) {
+        const NormalCollision& c = collisions[i];
+        if (!c.is_mollified()) {
+            continue;
+        }
+        const double m =
+            c.mollifier(c.dof(vertices, mesh.edges(), mesh.faces()));
+        if (m <= 0) {
+            m_le_0_count++;
+            if (c.weight < 0) {
+                m_le_0_negative_weight_count++;
+            }
+        }
+    }
+    REQUIRE(m_le_0_count > 0);
+    REQUIRE(m_le_0_negative_weight_count > 0);
+
+    const BarrierPotential barrier_potential(dhat, /*stiffness=*/1.0);
+
+    // Triggers the m <= 0 early-return (and its debug PSD assert) for every
+    // mollified collision that reaches it.
+    const Eigen::SparseMatrix<double> H = barrier_potential.hessian(
+        collisions, mesh, vertices, PSDProjectionMethod::CLAMP);
+    REQUIRE(H.nonZeros() > 0);
+
+    // With per-block CLAMP projection (including the m <= 0 blocks), the
+    // assembled hessian must be PSD. Before the fix, the m <= 0 branch returned
+    // negative-definite blocks unprojected, making this fail.
+    const Eigen::MatrixXd H_dense = H.toDense();
+    const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(H_dense);
+    const double min_eig = eigensolver.eigenvalues().minCoeff();
+    const double scale = eigensolver.eigenvalues().cwiseAbs().maxCoeff();
+    CHECK(min_eig >= -1e-10 * std::max(scale, 1.0));
+}
+
 // -- Benchmarking ------------------------------------------------------------
 
 TEST_CASE(