Implement Preconditioner from Shen et al. [2024]

python/src/candidates/collision_stencil.cpp

@@ -255,5 +255,138 @@ void define_collision_stencil(py::module_& m)
                                     vertices_t0: Stencil vertices at the start of the time step.
                     vertices_t1: Stencil vertices at the end of the time step.
                 )ipc_Qu8mg5v7",
-            "vertices_t0"_a, "vertices_t1"_a);
+            "vertices_t0"_a, "vertices_t1"_a)
+        .def(
+            "compute_distance_vector",
+            py::overload_cast<Eigen::ConstRef<VectorMax12d>>(
+                &CollisionStencil::compute_distance_vector, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the distance vector of the stencil: t = sum(c_i * x_i).
+
+            The distance vector is the vector between the closest points on the
+            collision primitives. Its squared norm equals the squared distance.
+
+            Note:
+                positions can be computed as stencil.dof(vertices, edges, faces)
+
+            Parameters:
+                positions: Stencil's vertex positions.
+
+            Returns:
+                The distance vector (dim-dimensional, i.e., 2D or 3D).
+            )ipc_Qu8mg5v7",
+            "positions"_a)
+        .def(
+            "compute_distance_vector_with_coefficients",
+            [](const CollisionStencil& self,
+               Eigen::ConstRef<VectorMax12d> positions) {
+                VectorMax4d coeffs;
+                VectorMax3d dv =
+                    self.compute_distance_vector(positions, coeffs);
+                return std::make_tuple(dv, coeffs);
+            },
+            R"ipc_Qu8mg5v7(
+            Compute the distance vector and the coefficients together.
+
+            Note:
+                positions can be computed as stencil.dof(vertices, edges, faces)
+
+            Parameters:
+                positions: Stencil's vertex positions.
+
+            Returns:
+                Tuple of:
+                The distance vector (dim-dimensional).
+                The computed coefficients c_i.
+            )ipc_Qu8mg5v7",
+            "positions"_a)
+        .def(
+            "compute_distance_vector",
+            py::overload_cast<
+                Eigen::ConstRef<Eigen::MatrixXd>,
+                Eigen::ConstRef<Eigen::MatrixXi>,
+                Eigen::ConstRef<Eigen::MatrixXi>>(
+                &CollisionStencil::compute_distance_vector, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the distance vector of the stencil.
+
+            Parameters:
+                vertices: Collision mesh vertices.
+                edges: Collision mesh edges.
+                faces: Collision mesh faces.
+
+            Returns:
+                The distance vector (dim-dimensional).
+            )ipc_Qu8mg5v7",
+            "vertices"_a, "edges"_a, "faces"_a)
+        .def(
+            "compute_distance_vector_jacobian",
+            &CollisionStencil::compute_distance_vector_jacobian,
+            R"ipc_Qu8mg5v7(
+            Compute the Jacobian of the distance vector w.r.t. positions.
+
+            J = [c_0 I, c_1 I, ..., c_n I]^T where I is the dim x dim identity.
+
+            Note:
+                positions can be computed as stencil.dof(vertices, edges, faces)
+
+            Parameters:
+                positions: Stencil's vertex positions.
+
+            Returns:
+                The Jacobian dt/dx as a matrix of shape (ndof, dim).
+            )ipc_Qu8mg5v7",
+            "positions"_a)
+        .def_static(
+            "diag_distance_vector_outer",
+            &CollisionStencil::diag_distance_vector_outer,
+            R"ipc_Qu8mg5v7(
+            Compute diag((dt/dx)(dt/dx)^T) efficiently (Eq. 11).
+
+            Result is [c_0^2, c_0^2, c_0^2, c_1^2, ...] (each c_i^2 repeated
+            dim times).
+
+            Parameters:
+                coeffs: The coefficients c_i (from compute_coefficients).
+                dim: The spatial dimension (2 or 3).
+
+            Returns:
+                The diagonal of (dt/dx)(dt/dx)^T as a vector of size ndof.
+            )ipc_Qu8mg5v7",
+            "coeffs"_a, "dim"_a)
+        .def_static(
+            "diag_distance_vector_t_outer",
+            &CollisionStencil::diag_distance_vector_t_outer,
+            R"ipc_Qu8mg5v7(
+            Compute diag((dt/dx * t)(dt/dx * t)^T) efficiently (Eq. 12).
+
+            Result is element-wise square of [c_0*t^T, c_1*t^T, ..., c_n*t^T].
+
+            Parameters:
+                coeffs: The coefficients c_i.
+                distance_vector: The distance vector t.
+
+            Returns:
+                The diagonal of (dt/dx*t)(dt/dx*t)^T as a vector of size ndof.
+            )ipc_Qu8mg5v7",
+            "coeffs"_a, "distance_vector"_a)
+        .def_static(
+            "contract_distance_vector_jacobian",
+            &CollisionStencil::contract_distance_vector_jacobian,
+            R"ipc_Qu8mg5v7(
+            Compute p^T (dt/dx) efficiently as sum(c_i * p_i) (Eqs. 13-14).
+
+            Given p = [p_0, p_1, ..., p_n]^T where p_i are dim-dimensional,
+            this computes p^T (dt/dx) = sum(c_i * p_i) which is a
+            dim-dimensional vector.
+
+            Parameters:
+                coeffs: The coefficients c_i.
+                p: A vector of size ndof (the direction for the quadratic form).
+                dim: The spatial dimension (2 or 3).
+
+            Returns:
+                p^T (dt/dx) as a dim-dimensional vector.
+            )ipc_Qu8mg5v7",
+            "coeffs"_a, "p"_a, "dim"_a);
 }

python/src/potentials/normal_potential.cpp

@@ -82,5 +82,110 @@ void define_normal_potential(py::module_& m)
             Returns:
                 The gradient of the force.
             )ipc_Qu8mg5v7",
-            "distance_squared"_a, "distance_squared_gradient"_a, "dmin"_a);
+            "distance_squared"_a, "distance_squared_gradient"_a, "dmin"_a)
+        .def(
+            "gauss_newton_hessian_diagonal",
+            py::overload_cast<
+                const NormalCollisions&, const CollisionMesh&,
+                Eigen::ConstRef<Eigen::MatrixXd>>(
+                &NormalPotential::gauss_newton_hessian_diagonal, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the diagonal of the cumulative Gauss-Newton Hessian of the potential.
+
+            Uses the distance-vector formulation to efficiently compute the
+            diagonal of the Gauss-Newton Hessian without forming full local
+            12x12 Hessian matrices. This is useful as a Jacobi preconditioner
+            for iterative solvers.
+
+            Note:
+                This is a Gauss-Newton approximation (drops derivatives of
+                closest-point coefficients), not the exact Hessian.
+
+            Parameters:
+                collisions: The set of collisions.
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+
+            Returns:
+                The diagonal of the Gauss-Newton Hessian as a vector of size vertices.size().
+            )ipc_Qu8mg5v7",
+            "collisions"_a, "mesh"_a, "vertices"_a)
+        .def(
+            "gauss_newton_hessian_diagonal",
+            py::overload_cast<
+                const NormalCollision&, Eigen::ConstRef<VectorMax12d>>(
+                &NormalPotential::gauss_newton_hessian_diagonal, py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the diagonal of the Gauss-Newton Hessian for a single collision.
+
+            Uses the distance-vector formulation (Eqs. 10-12) to efficiently
+            compute the diagonal without forming the full local Hessian.
+
+            Note:
+                This is a Gauss-Newton approximation, not the exact Hessian diagonal.
+
+            Parameters:
+                collision: The collision.
+                positions: The collision stencil's positions.
+
+            Returns:
+                The diagonal of the Gauss-Newton Hessian as a vector of size ndof.
+            )ipc_Qu8mg5v7",
+            "collision"_a, "positions"_a)
+        .def(
+            "gauss_newton_hessian_quadratic_form",
+            py::overload_cast<
+                const NormalCollisions&, const CollisionMesh&,
+                Eigen::ConstRef<Eigen::MatrixXd>,
+                Eigen::ConstRef<Eigen::VectorXd>>(
+                &NormalPotential::gauss_newton_hessian_quadratic_form,
+                py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute the product p^T H p for the cumulative Gauss-Newton Hessian.
+
+            Uses the distance-vector formulation to efficiently compute the
+            quadratic form without forming full local 12x12 Hessian matrices
+            nor the global sparse Hessian. This is useful for nonlinear
+            conjugate gradient methods.
+
+            Note:
+                This is a Gauss-Newton approximation (drops derivatives of
+                closest-point coefficients), not the exact Hessian.
+
+            Parameters:
+                collisions: The set of collisions.
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+                p: The direction vector of size vertices.size().
+
+            Returns:
+                The scalar value p^T H p (approximate).
+            )ipc_Qu8mg5v7",
+            "collisions"_a, "mesh"_a, "vertices"_a, "p"_a)
+        .def(
+            "gauss_newton_hessian_quadratic_form",
+            py::overload_cast<
+                const NormalCollision&, Eigen::ConstRef<VectorMax12d>,
+                Eigen::ConstRef<VectorMax12d>>(
+                &NormalPotential::gauss_newton_hessian_quadratic_form,
+                py::const_),
+            R"ipc_Qu8mg5v7(
+            Compute p^T H p for a single collision using the Gauss-Newton Hessian.
+
+            Uses the distance-vector formulation (Eqs. 10, 13-14) to
+            efficiently compute the quadratic form without forming the full
+            local Hessian.
+
+            Note:
+                This is a Gauss-Newton approximation, not the exact Hessian.
+
+            Parameters:
+                collision: The collision.
+                positions: The collision stencil's positions.
+                p: The local direction vector (size ndof).
+
+            Returns:
+                The scalar value p^T H p (approximate).
+            )ipc_Qu8mg5v7",
+            "collision"_a, "positions"_a, "p"_a);
 }

src/ipc/candidates/collision_stencil.cpp

@@ -4,6 +4,82 @@
 
 namespace ipc {
 
+VectorMax3d CollisionStencil::compute_distance_vector(
+    Eigen::ConstRef<VectorMax12d> positions) const
+{
+    VectorMax4d _; // Unused output
+    return compute_distance_vector(positions, _);
+}
+
+VectorMax3d CollisionStencil::compute_distance_vector(
+    Eigen::ConstRef<VectorMax12d> positions, VectorMax4d& coeffs) const
+{
+    const int n = num_vertices();
+    const int d = dim(positions.size());
+    coeffs = compute_coefficients(positions);
+    VectorMax3d dv = VectorMax3d::Zero(d);
+    for (int i = 0; i < n; i++) {
+        dv += coeffs[i] * positions.segment(d * i, d);
+    }
+    return dv;
+}
+
+MatrixMax<double, 12, 3> CollisionStencil::compute_distance_vector_jacobian(
+    Eigen::ConstRef<VectorMax12d> positions) const
+{
+    const int n = num_vertices();
+    const int d = dim(positions.size());
+    const int ndof = n * d;
+    const VectorMax4d c = compute_coefficients(positions);
+    // ∂t/∂x is ndof × dim
+    // Each block row i is cᵢ * I_{d×d}
+    MatrixMax<double, 12, 3> J;
+    J.setZero(ndof, d);
+    for (int i = 0; i < n; i++) {
+        J.block(i * d, 0, d, d).diagonal().array() = c[i];
+    }
+    return J;
+}
+
+VectorMax12d CollisionStencil::diag_distance_vector_outer(
+    Eigen::ConstRef<VectorMax4d> coeffs, const int d)
+{
+    const int n = coeffs.size();
+    VectorMax12d diag(n * d);
+    for (int i = 0; i < n; i++) {
+        const double ci2 = coeffs[i] * coeffs[i];
+        diag.segment(i * d, d).setConstant(ci2);
+    }
+    return diag;
+}
+
+VectorMax12d CollisionStencil::diag_distance_vector_t_outer(
+    Eigen::ConstRef<VectorMax4d> coeffs,
+    Eigen::ConstRef<VectorMax3d> distance_vector)
+{
+    const int n = coeffs.size();
+    const int d = distance_vector.size();
+    const VectorMax3d t2 = distance_vector.array().square();
+    VectorMax12d diag(n * d);
+    for (int i = 0; i < n; i++) {
+        diag.segment(i * d, d).array() = (coeffs[i] * coeffs[i]) * t2;
+    }
+    return diag;
+}
+
+VectorMax3d CollisionStencil::contract_distance_vector_jacobian(
+    Eigen::ConstRef<VectorMax4d> coeffs,
+    Eigen::ConstRef<VectorMax12d> p,
+    const int d)
+{
+    const int n = coeffs.size();
+    VectorMax3d result = VectorMax3d::Zero(d);
+    for (int i = 0; i < n; i++) {
+        result += coeffs[i] * p.segment(d * i, d);
+    }
+    return result;
+}
+
 VectorMax3d CollisionStencil::compute_normal(
     Eigen::ConstRef<VectorMax12d> positions,
     bool flip_if_negative,