Rest shape measure

src/ipc/smooth_contact/common.hpp

@@ -99,6 +99,12 @@ struct SmoothContactParameters {
     double beta_n = 0;
     int r = 2;
 
+    // If true, use the rest-shape measure for quadrature so weights stay
+    // constant in the rest configuration and measure derivatives are omitted.
+    // Enable this when you want a rest-shape-based measure instead of one that
+    // varies with the current deformed shape.
+    bool use_rest_shape_measure = false;
+
 private:
     double m_adaptive_dhat_ratio = 0.5;
 };

src/ipc/smooth_contact/primitives/edge2.cpp

@@ -13,6 +13,11 @@ Edge2::Edge2(
 {
     m_vertex_ids = { { mesh.edges()(id, 0), mesh.edges()(id, 1) } };
 
+    // Rest-shape edge length: computed from rest positions — constant
+    // quadrature weight
+    const Eigen::MatrixXd& rp = mesh.rest_positions();
+    m_rest_length = (rp.row(m_vertex_ids[1]) - rp.row(m_vertex_ids[0])).norm();
+
     if (mesh.is_orient_vertex(m_vertex_ids[0])
         && mesh.is_orient_vertex(m_vertex_ids[1])) {
         m_is_active =
@@ -32,45 +37,57 @@ double Edge2::potential(
     Eigen::ConstRef<Eigen::Vector4d> x) const
 {
     assert(m_is_active);
-    return (x.tail<2>() - x.head<2>()).norm();
+    if (!m_params.use_rest_shape_measure) {
+        return (x.tail<2>() - x.head<2>()).norm();
+    }
+    // Return the rest-shape edge length — constant quadrature weight per the
+    // paper (Eq. 13)
+    return m_rest_length;
 }
 
 Vector6d Edge2::grad(
     Eigen::ConstRef<Eigen::Vector2d> d,
     Eigen::ConstRef<Eigen::Vector4d> x) const
 {
     assert(m_is_active);
-    const Eigen::Vector2d t = x.tail<2>() - x.head<2>();
-    const double len = t.norm();
-    Vector6d g;
-    g.setZero();
-    g.segment<2>(2) = -t / len;
-    g.segment<2>(4) = t / len;
-    return g;
+    if (!m_params.use_rest_shape_measure) {
+        const Eigen::Vector2d t = x.tail<2>() - x.head<2>();
+        const double len = t.norm();
+        Vector6d g;
+        g.setZero();
+        g.segment<2>(2) = -t / len;
+        g.segment<2>(4) = t / len;
+        return g;
+    }
+    // Rest length is constant — no gradient contribution from this primitive
+    return Vector6d::Zero();
 }
 
 Matrix6d Edge2::hessian(
     Eigen::ConstRef<Eigen::Vector2d> d,
     Eigen::ConstRef<Eigen::Vector4d> x) const
 {
     assert(m_is_active);
-    Matrix6d h;
-    h.setZero();
+    if (!m_params.use_rest_shape_measure) {
+        Matrix6d h;
+        h.setZero();
 #ifdef IPC_TOOLKIT_DEBUG_AUTODIFF
-    ScalarBase::setVariableCount(4);
-    using T = ADHessian<4>;
-    auto xAD = slice_positions<T, 2, 2>(x);
-    h.block<4, 4>(2, 2) = (xAD.row(0) - xAD.row(1)).norm().Hess;
+        ScalarBase::setVariableCount(4);
+        using T = ADHessian<4>;
+        auto xAD = slice_positions<T, 2, 2>(x);
+        h.block<4, 4>(2, 2) = (xAD.row(0) - xAD.row(1)).norm().Hess;
 #else
-    const Eigen::Vector2d t = x.tail<2>() - x.head<2>();
-    const double norm = t.norm();
-    h.block<2, 2>(2, 2) =
-        (Eigen::Matrix2d::Identity() - t * (1. / norm / norm) * t.transpose())
-        / norm;
-    h.block<2, 2>(4, 4) = h.block<2, 2>(2, 2);
-    h.block<2, 2>(2, 4) = -h.block<2, 2>(2, 2);
-    h.block<2, 2>(4, 2) = -h.block<2, 2>(2, 2);
+        const Eigen::Vector2d t = x.tail<2>() - x.head<2>();
+        const double norm = t.norm();
+        h.block<2, 2>(2, 2) = (Eigen::Matrix2d::Identity()
+                               - t * (1. / norm / norm) * t.transpose())
+            / norm;
+        h.block<2, 2>(4, 4) = h.block<2, 2>(2, 2);
+        h.block<2, 2>(2, 4) = -h.block<2, 2>(2, 2);
+        h.block<2, 2>(4, 2) = -h.block<2, 2>(2, 2);
 #endif
-    return h;
+        return h;
+    }
+    return Matrix6d::Zero();
 }
 } // namespace ipc

src/ipc/smooth_contact/primitives/edge2.hpp

@@ -30,5 +30,9 @@ class Edge2 : public Primitive {
     Matrix6d hessian(
         Eigen::ConstRef<Eigen::Vector2d> d,
         Eigen::ConstRef<Eigen::Vector4d> x) const;
+
+private:
+    double
+        m_rest_length; ///< Rest-shape edge length (quadrature weight, constant)
 };
 } // namespace ipc

src/ipc/smooth_contact/primitives/edge3.cpp

@@ -82,6 +82,13 @@ Edge3::Edge3(
         faces.row(i) = face_rows[i];
     }
 
+    // Rest-shape squared edge length: from rest positions — constant quadrature
+    // weight
+    {
+        const Eigen::MatrixXd& rp = mesh.rest_positions();
+        m_rest_sq_length = (rp.row(e1_id) - rp.row(e0_id)).squaredNorm();
+    }
+
     if (m_vertex_ids.size() > N_EDGE_NEIGHBORS_3D) {
         logger().error(
             "Too many vertex neighbors for edge3 primitive! vertex count {} exceeds N_EDGE_NEIGHBORS_3D {}. "
@@ -193,8 +200,9 @@ T Edge3::smooth_edge3_term(
         normal_term = Math<T>::smooth_heaviside(normal_term - 1, 1., 0);
     }
 
-    // Weight: squared edge length
-    const T weight = (e1 - e0).squaredNorm();
+    // Weight: squared edge length (rest-shape if enabled, deformed otherwise)
+    const T weight = m_params.use_rest_shape_measure ? T(m_rest_sq_length)
+                                                     : (e1 - e0).squaredNorm();
 
     return weight * tangent_term * normal_term;
 }
@@ -559,17 +567,18 @@ GradientType<Eigen::Dynamic> Edge3::smooth_edge3_term_gradient(
     Eigen::VectorXd grad_tmp =
         tangent_grad * normal_term + normal_grad * tangent_term;
 
-    // Weight = (e1 - e0).squaredNorm()
     const Eigen::RowVector3d edge = X.row(1) - X.row(0);
-    const double weight = edge.squaredNorm();
+    const double weight =
+        m_params.use_rest_shape_measure ? m_rest_sq_length : edge.squaredNorm();
     grad_tmp *= weight;
-    // Derivative of weight w.r.t. e0 and e1
-    // d/d(e0) (e1-e0)^2 = 2*(e0-e1), d/d(e1) (e1-e0)^2 = 2*(e1-e0)
-    const int id_e0 = 3; // offset in [direction, X]
-    const int id_e1 = 6;
-    grad_tmp.segment<3>(id_e0) += (2 * val) * (-edge).transpose();
-    grad_tmp.segment<3>(id_e1) += (2 * val) * edge.transpose();
-
+    if (!m_params.use_rest_shape_measure) {
+        // Derivative of weight w.r.t. e0 and e1
+        // d/d(e0) (e1-e0)^2 = 2*(e0-e1), d/d(e1) (e1-e0)^2 = 2*(e1-e0)
+        const int id_e0 = 3; // offset in [direction, X]
+        const int id_e1 = 6;
+        grad_tmp.segment<3>(id_e0) += (2 * val) * (-edge).transpose();
+        grad_tmp.segment<3>(id_e1) += (2 * val) * edge.transpose();
+    }
     val *= weight;
 
     return std::make_tuple(val, grad_tmp);

src/ipc/smooth_contact/primitives/edge3.hpp

@@ -177,6 +177,8 @@ class Edge3 : public Primitive {
 
     /// @brief Whether the edge is orientable. If false, the normal term is not applied.
     bool orientable;
+    double m_rest_sq_length; ///< Rest-shape squared edge length (quadrature
+                             ///< weight, constant)
 };
 
 } // namespace ipc

src/ipc/smooth_contact/primitives/face.cpp

@@ -30,6 +30,15 @@ Face::Face(
 {
     m_vertex_ids = { { mesh.faces()(id, 0), mesh.faces()(id, 1),
                        mesh.faces()(id, 2) } };
+
+    // Rest-shape area: computed from rest positions — constant quadrature
+    // weight
+    const Eigen::MatrixXd& rp = mesh.rest_positions();
+    Eigen::Vector3d ra = rp.row(m_vertex_ids[1]) - rp.row(m_vertex_ids[0]);
+    Eigen::Vector3d rb = rp.row(m_vertex_ids[2]) - rp.row(m_vertex_ids[0]);
+    m_rest_area = 0.5 * ra.cross(rb).norm();
+
+    // Active check uses deformed positions / direction
     Eigen::Vector3d a =
         vertices.row(m_vertex_ids[1]) - vertices.row(m_vertex_ids[0]);
     Eigen::Vector3d b =
@@ -41,31 +50,44 @@ Face::Face(
     m_is_active = !orientable || a.cross(b).dot(d) > 0;
 }
 int Face::n_vertices() const { return N_FACE_NEIGHBORS_3D; }
-double
-Face::potential(Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x)
+double Face::potential(
+    Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x) const
 {
-    return smooth_face_term<double>(x.head<3>(), x.segment<3>(3), x.tail<3>());
+    if (!m_params.use_rest_shape_measure) {
+        return smooth_face_term<double>(
+            x.head<3>(), x.segment<3>(3), x.tail<3>());
+    }
+    // Return the rest-shape area — constant quadrature weight per the paper
+    // (Eq. 13)
+    return m_rest_area;
 }
-Vector12d
-Face::grad(Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x)
+Vector12d Face::grad(
+    Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x) const
 {
-    Vector12d g;
-    g.setZero();
-    ScalarBase::setVariableCount(9);
-    auto X = slice_positions<ADGrad<9>, 3, 3>(x);
-    g.tail<9>() =
-        smooth_face_term<ADGrad<9>>(X.row(0), X.row(1), X.row(2)).grad;
-    return g;
+    if (!m_params.use_rest_shape_measure) {
+        Vector12d g;
+        g.setZero();
+        ScalarBase::setVariableCount(9);
+        auto X = slice_positions<ADGrad<9>, 3, 3>(x);
+        g.tail<9>() =
+            smooth_face_term<ADGrad<9>>(X.row(0), X.row(1), X.row(2)).grad;
+        return g;
+    }
+    // Rest area is constant — no gradient contribution from this primitive
+    return Vector12d::Zero();
 }
-Matrix12d
-Face::hessian(Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x)
+Matrix12d Face::hessian(
+    Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x) const
 {
-    Matrix12d h;
-    h.setZero();
-    ScalarBase::setVariableCount(9);
-    auto X = slice_positions<ADHessian<9>, 3, 3>(x);
-    h.bottomRightCorner<9, 9>() =
-        smooth_face_term<ADHessian<9>>(X.row(0), X.row(1), X.row(2)).Hess;
-    return h;
+    if (!m_params.use_rest_shape_measure) {
+        Matrix12d h;
+        h.setZero();
+        ScalarBase::setVariableCount(9);
+        auto X = slice_positions<ADHessian<9>, 3, 3>(x);
+        h.bottomRightCorner<9, 9>() =
+            smooth_face_term<ADHessian<9>>(X.row(0), X.row(1), X.row(2)).Hess;
+        return h;
+    }
+    return Matrix12d::Zero();
 }
 } // namespace ipc

src/ipc/smooth_contact/primitives/face.hpp

@@ -21,12 +21,15 @@ class Face : public Primitive {
     int n_vertices() const override;
     int n_dofs() const override { return n_vertices() * DIM; }
 
-    static double
-    potential(Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x);
-    static Vector12d
-    grad(Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x);
-    static Matrix12d
-    hessian(Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x);
+    double potential(
+        Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x) const;
+    Vector12d
+    grad(Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x) const;
+    Matrix12d hessian(
+        Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x) const;
+
+private:
+    double m_rest_area; ///< Rest-shape face area (quadrature weight, constant)
 };
 
 /// @brief d points from triangle to the point