diff --git a/src/ipc/potentials/normal_potential.cpp b/src/ipc/potentials/normal_potential.cpp index 079cbd3c3..96b1b0e6b 100644 --- a/src/ipc/potentials/normal_potential.cpp +++ b/src/ipc/potentials/normal_potential.cpp @@ -8,6 +8,9 @@ #include #include +#include +#include + 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) diff --git a/tests/src/tests/potential/test_barrier_potential.cpp b/tests/src/tests/potential/test_barrier_potential.cpp index 2bb56f509..5a82f0403 100644 --- a/tests/src/tests/potential/test_barrier_potential.cpp +++ b/tests/src/tests/potential/test_barrier_potential.cpp @@ -15,6 +15,8 @@ #include #include +#include + 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 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 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(