From ee0a7961565ab89d88856c56c4cb996516ddaf74 Mon Sep 17 00:00:00 2001
From: Uday Kusupati <kusupatiuday@gmail.com>
Date: Wed, 22 Apr 2026 18:08:32 -0400
Subject: [PATCH 1/6] 2d bug fix

---
 src/ipc/smooth_contact/primitives/edge2.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/ipc/smooth_contact/primitives/edge2.cpp b/src/ipc/smooth_contact/primitives/edge2.cpp
index fd3870e02..40ffaebe9 100644
--- a/src/ipc/smooth_contact/primitives/edge2.cpp
+++ b/src/ipc/smooth_contact/primitives/edge2.cpp
@@ -15,7 +15,7 @@ Edge2::Edge2(
 
     m_is_active = (mesh.is_orient_vertex(m_vertex_ids[0])
                    && mesh.is_orient_vertex(m_vertex_ids[1]))
-        || Math<double>::cross2(
+        && Math<double>::cross2(
                d, vertices.row(m_vertex_ids[1]) - vertices.row(m_vertex_ids[0]))
             > 0;
 }

From db1e3a966970f77dc12534294f38452f0e259b0a Mon Sep 17 00:00:00 2001
From: Uday Kusupati <kusupatiuday@gmail.com>
Date: Wed, 22 Apr 2026 19:06:12 -0400
Subject: [PATCH 2/6] add rest_shape_measure for quadrature in GCP potential
 integration

---
 src/ipc/smooth_contact/common.hpp            |  2 +
 src/ipc/smooth_contact/primitives/edge2.cpp  | 62 ++++++++++++--------
 src/ipc/smooth_contact/primitives/edge2.hpp  |  3 +
 src/ipc/smooth_contact/primitives/edge3.cpp  | 31 ++++++----
 src/ipc/smooth_contact/primitives/edge3.hpp  |  1 +
 src/ipc/smooth_contact/primitives/face.cpp   | 61 ++++++++++++-------
 src/ipc/smooth_contact/primitives/face.hpp   | 12 ++--
 src/ipc/smooth_contact/primitives/point2.cpp | 57 ++++++++++++++----
 src/ipc/smooth_contact/primitives/point2.hpp |  1 +
 src/ipc/smooth_contact/primitives/point3.cpp | 52 +++++++++++-----
 src/ipc/smooth_contact/primitives/point3.hpp |  1 +
 11 files changed, 195 insertions(+), 88 deletions(-)

diff --git a/src/ipc/smooth_contact/common.hpp b/src/ipc/smooth_contact/common.hpp
index f7390f512..eac1f1115 100644
--- a/src/ipc/smooth_contact/common.hpp
+++ b/src/ipc/smooth_contact/common.hpp
@@ -99,6 +99,8 @@ struct SmoothContactParameters {
     double beta_n = 0;
     int r = 2;
 
+    bool use_rest_shape_measure = true;
+
 private:
     double m_adaptive_dhat_ratio = 0.5;
 };
diff --git a/src/ipc/smooth_contact/primitives/edge2.cpp b/src/ipc/smooth_contact/primitives/edge2.cpp
index 40ffaebe9..b56b09f7b 100644
--- a/src/ipc/smooth_contact/primitives/edge2.cpp
+++ b/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();
+
     m_is_active = (mesh.is_orient_vertex(m_vertex_ids[0])
                    && mesh.is_orient_vertex(m_vertex_ids[1]))
         && Math<double>::cross2(
@@ -27,7 +32,11 @@ 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(
@@ -35,13 +44,17 @@ Vector6d Edge2::grad(
     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(
@@ -49,23 +62,26 @@ Matrix6d Edge2::hessian(
     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
diff --git a/src/ipc/smooth_contact/primitives/edge2.hpp b/src/ipc/smooth_contact/primitives/edge2.hpp
index d9bfe5efa..7ac322229 100644
--- a/src/ipc/smooth_contact/primitives/edge2.hpp
+++ b/src/ipc/smooth_contact/primitives/edge2.hpp
@@ -30,5 +30,8 @@ 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
diff --git a/src/ipc/smooth_contact/primitives/edge3.cpp b/src/ipc/smooth_contact/primitives/edge3.cpp
index 884c355ff..35bb409f3 100644
--- a/src/ipc/smooth_contact/primitives/edge3.cpp
+++ b/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,10 @@ 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 +568,17 @@ 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);
diff --git a/src/ipc/smooth_contact/primitives/edge3.hpp b/src/ipc/smooth_contact/primitives/edge3.hpp
index a2f266e32..af1466e2f 100644
--- a/src/ipc/smooth_contact/primitives/edge3.hpp
+++ b/src/ipc/smooth_contact/primitives/edge3.hpp
@@ -177,6 +177,7 @@ 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
\ No newline at end of file
diff --git a/src/ipc/smooth_contact/primitives/face.cpp b/src/ipc/smooth_contact/primitives/face.cpp
index 9621b8f7e..03c97ecb3 100644
--- a/src/ipc/smooth_contact/primitives/face.cpp
+++ b/src/ipc/smooth_contact/primitives/face.cpp
@@ -30,6 +30,14 @@ 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 +49,42 @@ 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
diff --git a/src/ipc/smooth_contact/primitives/face.hpp b/src/ipc/smooth_contact/primitives/face.hpp
index db061b3f6..c3ab31112 100644
--- a/src/ipc/smooth_contact/primitives/face.hpp
+++ b/src/ipc/smooth_contact/primitives/face.hpp
@@ -21,12 +21,12 @@ 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
diff --git a/src/ipc/smooth_contact/primitives/point2.cpp b/src/ipc/smooth_contact/primitives/point2.cpp
index 9d7e0ba79..8d4bed30b 100644
--- a/src/ipc/smooth_contact/primitives/point2.cpp
+++ b/src/ipc/smooth_contact/primitives/point2.cpp
@@ -61,7 +61,8 @@ namespace {
         Eigen::ConstRef<Eigen::Vector2<scalar>> e0,
         Eigen::ConstRef<Eigen::Vector2<scalar>> e1,
         const SmoothContactParameters& params,
-        const bool orientable)
+        const bool orientable,
+        const scalar measure)
     {
         const Eigen::Vector2<scalar> dn = -direc.normalized();
         const Eigen::Vector2<scalar> t0 = (e0 - v).normalized(),
@@ -83,7 +84,7 @@ namespace {
                 * Math<scalar>::smooth_heaviside(tmp - 1., params.alpha_n, 0);
         }
 
-        return val * ((e0 - v).norm() + (e1 - v).norm()) / 2.;
+        return val * measure;
     }
 
     template <class scalar>
@@ -91,7 +92,8 @@ namespace {
         Eigen::ConstRef<Eigen::Vector2<scalar>> v,
         Eigen::ConstRef<Eigen::Vector2<scalar>> direc,
         Eigen::ConstRef<Eigen::Vector2<scalar>> e0,
-        const SmoothContactParameters& params)
+        const SmoothContactParameters& params,
+        const scalar measure)
     {
         const Eigen::Vector2<scalar> dn = -direc.normalized();
         const Eigen::Vector2<scalar> t0 = e0 - v;
@@ -99,7 +101,7 @@ namespace {
         const scalar tangent_term = Math<scalar>::smooth_heaviside(
             dn.dot(t0) / t0.norm(), params.alpha_t, params.beta_t);
 
-        return tangent_term * t0.norm();
+        return tangent_term * measure;
     }
 } // namespace
 
@@ -121,6 +123,13 @@ Point2::Point2(
         m_is_active = smooth_point2_term_type(
             vertices.row(id), d, vertices.row(m_vertex_ids[1]),
             vertices.row(m_vertex_ids[2]), params, orientable);
+
+        // rest-shape measure: average half-edge length from rest positions
+        const Eigen::MatrixXd& rp = mesh.rest_positions();
+        m_rest_measure =
+            ((rp.row(neighbor_verts[0]) - rp.row(id)).norm()
+             + (rp.row(neighbor_verts[1]) - rp.row(id)).norm())
+            / 2.0;
     } else if (has_neighbor_1 || has_neighbor_2) {
         m_vertex_ids = {
             { id, has_neighbor_1 ? neighbor_verts[0] : neighbor_verts[1] }
@@ -132,10 +141,15 @@ Point2::Point2(
                 .normalized();
 
         m_is_active = dn.dot(t0) > -params.alpha_t;
+        // rest-shape measure: half-edge length from rest positions
+        const Eigen::MatrixXd& rp = mesh.rest_positions();
+        m_rest_measure =
+            (rp.row(m_vertex_ids[1]) - rp.row(m_vertex_ids[0])).norm();
     } else {
         m_vertex_ids.resize(1);
         m_vertex_ids[0] = id;
         m_is_active = true;
+        m_rest_measure = 1.0; // no neighbors — unit weight
     }
 }
 
@@ -146,14 +160,21 @@ double Point2::potential(
     const VectorMax<double, MAX_SIZE>& x) const
 {
     if (has_neighbor_1 && has_neighbor_2) {
+        const double measure = m_params.use_rest_shape_measure
+            ? m_rest_measure
+            : ((x.segment<DIM>(DIM) - x.segment<DIM>(0)).norm()
+               + (x.segment<DIM>(2 * DIM) - x.segment<DIM>(0)).norm()) / 2.0;
         return smooth_point2_term<double>(
             x.segment<DIM>(0), d, x.segment<DIM>(DIM), x.segment<DIM>(2 * DIM),
-            m_params, orientable);
+            m_params, orientable, measure);
     } else if (has_neighbor_1 || has_neighbor_2) {
+        const double measure = m_params.use_rest_shape_measure
+            ? m_rest_measure
+            : (x.segment<DIM>(DIM) - x.segment<DIM>(0)).norm();
         return smooth_point2_term_one_side<double>(
-            x.segment<DIM>(0), d, x.segment<DIM>(DIM), m_params);
+            x.segment<DIM>(0), d, x.segment<DIM>(DIM), m_params, measure);
     } else {
-        return 1.;
+        return m_rest_measure;
     }
 }
 VectorMax<double, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
@@ -166,8 +187,12 @@ VectorMax<double, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
         Eigen::Vector<double, 4 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
+        const T measure = m_params.use_rest_shape_measure
+            ? T(m_rest_measure)
+            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm()) / 2.0;
         return smooth_point2_term<T>(
-                   X.row(1), X.row(0), X.row(2), X.row(3), m_params, orientable)
+                   X.row(1), X.row(0), X.row(2), X.row(3), m_params,
+                   orientable, measure)
             .grad;
     } else if (has_neighbor_1 || has_neighbor_2) {
         ScalarBase::setVariableCount(3 * DIM);
@@ -175,8 +200,11 @@ VectorMax<double, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
         Eigen::Vector<double, 3 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 3, DIM> X = slice_positions<T, 3, DIM>(tmp);
+        const T measure = m_params.use_rest_shape_measure
+            ? T(m_rest_measure)
+            : (X.row(2) - X.row(0)).norm();
         return smooth_point2_term_one_side<T>(
-                   X.row(1), X.row(0), X.row(2), m_params)
+                   X.row(1), X.row(0), X.row(2), m_params, measure)
             .grad;
     } else {
         return VectorMax<double, Point2::MAX_SIZE + Point2::DIM>::Zero(
@@ -197,8 +225,12 @@ Point2::hessian(
         Eigen::Vector<double, 4 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
+        const T measure = m_params.use_rest_shape_measure
+            ? T(m_rest_measure)
+            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm()) / 2.0;
         return smooth_point2_term<T>(
-                   X.row(1), X.row(0), X.row(2), X.row(3), m_params, orientable)
+                   X.row(1), X.row(0), X.row(2), X.row(3), m_params,
+                   orientable, measure)
             .Hess;
     } else if (has_neighbor_1 || has_neighbor_2) {
         ScalarBase::setVariableCount(3 * DIM);
@@ -206,8 +238,11 @@ Point2::hessian(
         Eigen::Vector<double, 3 * DIM> tmp;
         tmp << d, x;
         Eigen::Matrix<T, 3, DIM> X = slice_positions<T, 3, DIM>(tmp);
+        const T measure = m_params.use_rest_shape_measure
+            ? T(m_rest_measure)
+            : (X.row(2) - X.row(0)).norm();
         return smooth_point2_term_one_side<T>(
-                   X.row(1), X.row(0), X.row(2), m_params)
+                   X.row(1), X.row(0), X.row(2), m_params, measure)
             .Hess;
     } else {
         return MatrixMax<double, -1, Point2::MAX_SIZE + Point2::DIM>::Zero(
diff --git a/src/ipc/smooth_contact/primitives/point2.hpp b/src/ipc/smooth_contact/primitives/point2.hpp
index a72ac2ea2..72baa5f27 100644
--- a/src/ipc/smooth_contact/primitives/point2.hpp
+++ b/src/ipc/smooth_contact/primitives/point2.hpp
@@ -41,5 +41,6 @@ class Point2 : public Primitive {
 private:
     bool has_neighbor_1, has_neighbor_2;
     bool orientable;
+    double m_rest_measure; ///< Rest-shape vertex length measure (quadrature weight, constant)
 };
 } // namespace ipc
diff --git a/src/ipc/smooth_contact/primitives/point3.cpp b/src/ipc/smooth_contact/primitives/point3.cpp
index 2ddb6391f..c298f9613 100644
--- a/src/ipc/smooth_contact/primitives/point3.cpp
+++ b/src/ipc/smooth_contact/primitives/point3.cpp
@@ -61,6 +61,20 @@ Point3::Point3(
     n_neighbors = local_to_global_vids.size() - 1;
     m_vertex_ids = local_to_global_vids;
 
+    if (m_params.use_rest_shape_measure) {
+        // Rest-shape vertex area measure: sum of squared rest 1-ring edge lengths / 3
+        // Constant quadrature weight per the paper (Eq. 13), never differentiated through
+        const Eigen::MatrixXd& rp = mesh.rest_positions();
+        m_rest_weight = 0.0;
+        for (int a = 0; a < edges.rows(); a++) {
+            const Eigen::RowVector3d t =
+                rp.row(local_to_global_vids[edges(a, 1)])
+                - rp.row(local_to_global_vids[edges(a, 0)]);
+            m_rest_weight += t.squaredNorm();
+        }
+        m_rest_weight /= 3.0;
+    }
+
     if (m_vertex_ids.size() > N_VERT_NEIGHBORS_3D) {
         logger().error(
             "Too many neighbors for point3 primitive! {} > {}! Increase N_VERT_NEIGHBORS_3D in common.hpp",
@@ -433,10 +447,11 @@ GradientType<-1> Point3::smooth_point3_term_gradient(
     Eigen::VectorXd grad_tmp =
         tangent_grad * normal_term + normal_grad * tangent_term;
 
-    const double weight = tangents.squaredNorm() / 3.;
+    const double weight = m_params.use_rest_shape_measure ? m_rest_weight : tangents.squaredNorm() / 3.;
     grad_tmp *= weight;
-    grad_tmp.tail(n_neighbor_dofs) += (2. / 3. * val) * tangents_vec;
-
+    if (!m_params.use_rest_shape_measure) {
+        grad_tmp.tail(n_neighbor_dofs) += (2. / 3. * val) * tangents_vec;
+    }
     val *= weight;
 
     // gradient wrt. [direc, v, neighbors]
@@ -489,19 +504,21 @@ HessianType<-1> Point3::smooth_point3_term_hessian(
         + normal_hess * tangent_term + tangent_grad * normal_grad.transpose()
         + normal_grad * tangent_grad.transpose();
 
-    const double weight = tangents.squaredNorm() / 3.;
+    const double weight = m_params.use_rest_shape_measure ? m_rest_weight : tangents.squaredNorm() / 3.;
     hess_tmp *= weight;
-    hess_tmp.bottomRightCorner(n_neighbor_dofs, n_neighbor_dofs)
-        .diagonal()
-        .array() += 2. / 3. * val;
-    hess_tmp.bottomRows(n_neighbor_dofs) +=
-        2. / 3. * tangents_vec * grad_tmp.transpose();
-    hess_tmp.rightCols(n_neighbor_dofs) +=
-        2. / 3. * grad_tmp * tangents_vec.transpose();
-
+    if (!m_params.use_rest_shape_measure) {
+        hess_tmp.bottomRightCorner(n_neighbor_dofs, n_neighbor_dofs)
+            .diagonal()
+            .array() += 2. / 3. * val;
+        hess_tmp.bottomRows(n_neighbor_dofs) +=
+            2. / 3. * tangents_vec * grad_tmp.transpose();
+        hess_tmp.rightCols(n_neighbor_dofs) +=
+            2. / 3. * grad_tmp * tangents_vec.transpose();
+    }
     grad_tmp *= weight;
-    grad_tmp.tail(n_neighbor_dofs) += (2. / 3. * val) * tangents_vec;
-
+    if (!m_params.use_rest_shape_measure) {
+        grad_tmp.tail(n_neighbor_dofs) += (2. / 3. * val) * tangents_vec;
+    }
     val *= weight;
 
     // gradient wrt. [direc, v, neighbors]
@@ -560,9 +577,12 @@ scalar Point3::smooth_point3_term(
                     -dn.dot(t) / t.norm(), m_params.alpha_t, m_params.beta_t);
         }
 
-        weight = weight + t.squaredNorm();
+        if (!m_params.use_rest_shape_measure) {
+            weight = weight + t.squaredNorm();
+        }
     }
-    weight /= 3.;
+
+    weight = m_params.use_rest_shape_measure ? scalar(m_rest_weight) : weight / 3.;
 
     if (!orientable || otypes.normal_type(0) == HeavisideType::ONE) {
         normal_term = scalar(1.);
diff --git a/src/ipc/smooth_contact/primitives/point3.hpp b/src/ipc/smooth_contact/primitives/point3.hpp
index 015a70f00..d29b09d13 100644
--- a/src/ipc/smooth_contact/primitives/point3.hpp
+++ b/src/ipc/smooth_contact/primitives/point3.hpp
@@ -96,6 +96,7 @@ class Point3 : public Primitive {
     Eigen::MatrixX3i faces;
     Eigen::MatrixX2i edges;
     bool orientable;
+    double m_rest_weight; ///< Rest-shape vertex area measure (quadrature weight, constant)
 
     bool smooth_point3_term_type(
         Eigen::ConstRef<Eigen::MatrixX3d> X,

From 2ae359adf5a51138aa7d266ad6b2ec16d5ca132b Mon Sep 17 00:00:00 2001
From: Uday Kusupati <kusupatiuday@gmail.com>
Date: Thu, 23 Apr 2026 11:53:25 -0400
Subject: [PATCH 3/6] clang-tidy

---
 src/ipc/smooth_contact/primitives/edge2.cpp | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/ipc/smooth_contact/primitives/edge2.cpp b/src/ipc/smooth_contact/primitives/edge2.cpp
index b56b09f7b..93efde666 100644
--- a/src/ipc/smooth_contact/primitives/edge2.cpp
+++ b/src/ipc/smooth_contact/primitives/edge2.cpp
@@ -18,11 +18,13 @@ Edge2::Edge2(
     m_rest_length =
         (rp.row(m_vertex_ids[1]) - rp.row(m_vertex_ids[0])).norm();
 
-    m_is_active = (mesh.is_orient_vertex(m_vertex_ids[0])
-                   && mesh.is_orient_vertex(m_vertex_ids[1]))
-        && Math<double>::cross2(
-               d, vertices.row(m_vertex_ids[1]) - vertices.row(m_vertex_ids[0]))
-            > 0;
+    if (mesh.is_orient_vertex(m_vertex_ids[0])
+        && mesh.is_orient_vertex(m_vertex_ids[1])) {
+        m_is_active = Math<double>::cross2(
+            d, vertices.row(m_vertex_ids[1]) - vertices.row(m_vertex_ids[0])) > 0;
+    } else {
+        m_is_active = true;
+    }
 }
 
 int Edge2::n_vertices() const { return N_EDGE_NEIGHBORS_2D; }

From 778f0780c6a3606450412f3f3fd58697e14570b1 Mon Sep 17 00:00:00 2001
From: Uday Kusupati <kusupatiuday@gmail.com>
Date: Thu, 30 Apr 2026 20:29:15 -0400
Subject: [PATCH 4/6] tidy up clang

---
 src/ipc/smooth_contact/primitives/edge2.cpp  | 20 +++++++++++-------
 src/ipc/smooth_contact/primitives/edge2.hpp  |  3 ++-
 src/ipc/smooth_contact/primitives/edge3.cpp  | 14 ++++++-------
 src/ipc/smooth_contact/primitives/edge3.hpp  |  3 ++-
 src/ipc/smooth_contact/primitives/face.cpp   |  9 +++++---
 src/ipc/smooth_contact/primitives/face.hpp   |  9 +++++---
 src/ipc/smooth_contact/primitives/point2.cpp | 22 +++++++++++---------
 src/ipc/smooth_contact/primitives/point2.hpp |  3 ++-
 src/ipc/smooth_contact/primitives/point3.cpp | 18 ++++++++++------
 src/ipc/smooth_contact/primitives/point3.hpp |  3 ++-
 10 files changed, 63 insertions(+), 41 deletions(-)

diff --git a/src/ipc/smooth_contact/primitives/edge2.cpp b/src/ipc/smooth_contact/primitives/edge2.cpp
index 93efde666..c3cf9e26b 100644
--- a/src/ipc/smooth_contact/primitives/edge2.cpp
+++ b/src/ipc/smooth_contact/primitives/edge2.cpp
@@ -13,15 +13,18 @@ Edge2::Edge2(
 {
     m_vertex_ids = { { mesh.edges()(id, 0), mesh.edges()(id, 1) } };
 
-    // Rest-shape edge length: computed from rest positions — constant quadrature weight
+    // 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();
+    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 = Math<double>::cross2(
-            d, vertices.row(m_vertex_ids[1]) - vertices.row(m_vertex_ids[0])) > 0;
+        m_is_active =
+            Math<double>::cross2(
+                d,
+                vertices.row(m_vertex_ids[1]) - vertices.row(m_vertex_ids[0]))
+            > 0;
     } else {
         m_is_active = true;
     }
@@ -37,7 +40,8 @@ double Edge2::potential(
     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 the rest-shape edge length — constant quadrature weight per the
+    // paper (Eq. 13)
     return m_rest_length;
 }
 
@@ -75,8 +79,8 @@ Matrix6d Edge2::hessian(
 #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())
+        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);
diff --git a/src/ipc/smooth_contact/primitives/edge2.hpp b/src/ipc/smooth_contact/primitives/edge2.hpp
index 7ac322229..4e5a30a47 100644
--- a/src/ipc/smooth_contact/primitives/edge2.hpp
+++ b/src/ipc/smooth_contact/primitives/edge2.hpp
@@ -32,6 +32,7 @@ class Edge2 : public Primitive {
         Eigen::ConstRef<Eigen::Vector4d> x) const;
 
 private:
-    double m_rest_length; ///< Rest-shape edge length (quadrature weight, constant)
+    double
+        m_rest_length; ///< Rest-shape edge length (quadrature weight, constant)
 };
 } // namespace ipc
diff --git a/src/ipc/smooth_contact/primitives/edge3.cpp b/src/ipc/smooth_contact/primitives/edge3.cpp
index 35bb409f3..54962e2ba 100644
--- a/src/ipc/smooth_contact/primitives/edge3.cpp
+++ b/src/ipc/smooth_contact/primitives/edge3.cpp
@@ -82,11 +82,11 @@ Edge3::Edge3(
         faces.row(i) = face_rows[i];
     }
 
-    // Rest-shape squared edge length: from rest positions — constant quadrature weight
+    // 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();
+        m_rest_sq_length = (rp.row(e1_id) - rp.row(e0_id)).squaredNorm();
     }
 
     if (m_vertex_ids.size() > N_EDGE_NEIGHBORS_3D) {
@@ -201,9 +201,8 @@ T Edge3::smooth_edge3_term(
     }
 
     // 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();
+    const T weight = m_params.use_rest_shape_measure ? T(m_rest_sq_length)
+                                                     : (e1 - e0).squaredNorm();
 
     return weight * tangent_term * normal_term;
 }
@@ -569,7 +568,8 @@ GradientType<Eigen::Dynamic> Edge3::smooth_edge3_term_gradient(
         tangent_grad * normal_term + normal_grad * tangent_term;
 
     const Eigen::RowVector3d edge = X.row(1) - X.row(0);
-    const double weight = m_params.use_rest_shape_measure ? m_rest_sq_length : edge.squaredNorm();
+    const double weight =
+        m_params.use_rest_shape_measure ? m_rest_sq_length : edge.squaredNorm();
     grad_tmp *= weight;
     if (!m_params.use_rest_shape_measure) {
         // Derivative of weight w.r.t. e0 and e1
diff --git a/src/ipc/smooth_contact/primitives/edge3.hpp b/src/ipc/smooth_contact/primitives/edge3.hpp
index af1466e2f..5ac56e32a 100644
--- a/src/ipc/smooth_contact/primitives/edge3.hpp
+++ b/src/ipc/smooth_contact/primitives/edge3.hpp
@@ -177,7 +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)
+    double m_rest_sq_length; ///< Rest-shape squared edge length (quadrature
+                             ///< weight, constant)
 };
 
 } // namespace ipc
\ No newline at end of file
diff --git a/src/ipc/smooth_contact/primitives/face.cpp b/src/ipc/smooth_contact/primitives/face.cpp
index 03c97ecb3..ce90bf6e0 100644
--- a/src/ipc/smooth_contact/primitives/face.cpp
+++ b/src/ipc/smooth_contact/primitives/face.cpp
@@ -31,7 +31,8 @@ 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
+    // 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]);
@@ -53,9 +54,11 @@ double Face::potential(
     Eigen::ConstRef<Eigen::Vector3d> d, Eigen::ConstRef<Vector9d> x) const
 {
     if (!m_params.use_rest_shape_measure) {
-        return smooth_face_term<double>(x.head<3>(), x.segment<3>(3), x.tail<3>());
+        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 the rest-shape area — constant quadrature weight per the paper
+    // (Eq. 13)
     return m_rest_area;
 }
 Vector12d Face::grad(
diff --git a/src/ipc/smooth_contact/primitives/face.hpp b/src/ipc/smooth_contact/primitives/face.hpp
index c3ab31112..a150c24be 100644
--- a/src/ipc/smooth_contact/primitives/face.hpp
+++ b/src/ipc/smooth_contact/primitives/face.hpp
@@ -21,9 +21,12 @@ class Face : public Primitive {
     int n_vertices() const override;
     int n_dofs() const override { return n_vertices() * DIM; }
 
-    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;
+    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)
diff --git a/src/ipc/smooth_contact/primitives/point2.cpp b/src/ipc/smooth_contact/primitives/point2.cpp
index 8d4bed30b..0691f333d 100644
--- a/src/ipc/smooth_contact/primitives/point2.cpp
+++ b/src/ipc/smooth_contact/primitives/point2.cpp
@@ -126,9 +126,8 @@ Point2::Point2(
 
         // rest-shape measure: average half-edge length from rest positions
         const Eigen::MatrixXd& rp = mesh.rest_positions();
-        m_rest_measure =
-            ((rp.row(neighbor_verts[0]) - rp.row(id)).norm()
-             + (rp.row(neighbor_verts[1]) - rp.row(id)).norm())
+        m_rest_measure = ((rp.row(neighbor_verts[0]) - rp.row(id)).norm()
+                          + (rp.row(neighbor_verts[1]) - rp.row(id)).norm())
             / 2.0;
     } else if (has_neighbor_1 || has_neighbor_2) {
         m_vertex_ids = {
@@ -163,7 +162,8 @@ double Point2::potential(
         const double measure = m_params.use_rest_shape_measure
             ? m_rest_measure
             : ((x.segment<DIM>(DIM) - x.segment<DIM>(0)).norm()
-               + (x.segment<DIM>(2 * DIM) - x.segment<DIM>(0)).norm()) / 2.0;
+               + (x.segment<DIM>(2 * DIM) - x.segment<DIM>(0)).norm())
+                / 2.0;
         return smooth_point2_term<double>(
             x.segment<DIM>(0), d, x.segment<DIM>(DIM), x.segment<DIM>(2 * DIM),
             m_params, orientable, measure);
@@ -189,10 +189,11 @@ VectorMax<double, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
         const T measure = m_params.use_rest_shape_measure
             ? T(m_rest_measure)
-            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm()) / 2.0;
+            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm())
+                / 2.0;
         return smooth_point2_term<T>(
-                   X.row(1), X.row(0), X.row(2), X.row(3), m_params,
-                   orientable, measure)
+                   X.row(1), X.row(0), X.row(2), X.row(3), m_params, orientable,
+                   measure)
             .grad;
     } else if (has_neighbor_1 || has_neighbor_2) {
         ScalarBase::setVariableCount(3 * DIM);
@@ -227,10 +228,11 @@ Point2::hessian(
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
         const T measure = m_params.use_rest_shape_measure
             ? T(m_rest_measure)
-            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm()) / 2.0;
+            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm())
+                / 2.0;
         return smooth_point2_term<T>(
-                   X.row(1), X.row(0), X.row(2), X.row(3), m_params,
-                   orientable, measure)
+                   X.row(1), X.row(0), X.row(2), X.row(3), m_params, orientable,
+                   measure)
             .Hess;
     } else if (has_neighbor_1 || has_neighbor_2) {
         ScalarBase::setVariableCount(3 * DIM);
diff --git a/src/ipc/smooth_contact/primitives/point2.hpp b/src/ipc/smooth_contact/primitives/point2.hpp
index 72baa5f27..56f57b9cd 100644
--- a/src/ipc/smooth_contact/primitives/point2.hpp
+++ b/src/ipc/smooth_contact/primitives/point2.hpp
@@ -41,6 +41,7 @@ class Point2 : public Primitive {
 private:
     bool has_neighbor_1, has_neighbor_2;
     bool orientable;
-    double m_rest_measure; ///< Rest-shape vertex length measure (quadrature weight, constant)
+    double m_rest_measure; ///< Rest-shape vertex length measure (quadrature
+                           ///< weight, constant)
 };
 } // namespace ipc
diff --git a/src/ipc/smooth_contact/primitives/point3.cpp b/src/ipc/smooth_contact/primitives/point3.cpp
index c298f9613..55a74c0ba 100644
--- a/src/ipc/smooth_contact/primitives/point3.cpp
+++ b/src/ipc/smooth_contact/primitives/point3.cpp
@@ -62,8 +62,9 @@ Point3::Point3(
     m_vertex_ids = local_to_global_vids;
 
     if (m_params.use_rest_shape_measure) {
-        // Rest-shape vertex area measure: sum of squared rest 1-ring edge lengths / 3
-        // Constant quadrature weight per the paper (Eq. 13), never differentiated through
+        // Rest-shape vertex area measure: sum of squared rest 1-ring edge
+        // lengths / 3 Constant quadrature weight per the paper (Eq. 13), never
+        // differentiated through
         const Eigen::MatrixXd& rp = mesh.rest_positions();
         m_rest_weight = 0.0;
         for (int a = 0; a < edges.rows(); a++) {
@@ -447,7 +448,9 @@ GradientType<-1> Point3::smooth_point3_term_gradient(
     Eigen::VectorXd grad_tmp =
         tangent_grad * normal_term + normal_grad * tangent_term;
 
-    const double weight = m_params.use_rest_shape_measure ? m_rest_weight : tangents.squaredNorm() / 3.;
+    const double weight = m_params.use_rest_shape_measure
+        ? m_rest_weight
+        : tangents.squaredNorm() / 3.;
     grad_tmp *= weight;
     if (!m_params.use_rest_shape_measure) {
         grad_tmp.tail(n_neighbor_dofs) += (2. / 3. * val) * tangents_vec;
@@ -504,7 +507,9 @@ HessianType<-1> Point3::smooth_point3_term_hessian(
         + normal_hess * tangent_term + tangent_grad * normal_grad.transpose()
         + normal_grad * tangent_grad.transpose();
 
-    const double weight = m_params.use_rest_shape_measure ? m_rest_weight : tangents.squaredNorm() / 3.;
+    const double weight = m_params.use_rest_shape_measure
+        ? m_rest_weight
+        : tangents.squaredNorm() / 3.;
     hess_tmp *= weight;
     if (!m_params.use_rest_shape_measure) {
         hess_tmp.bottomRightCorner(n_neighbor_dofs, n_neighbor_dofs)
@@ -561,7 +566,7 @@ HessianType<-1> Point3::smooth_point3_term_hessian(
 template <typename scalar, int n_verts>
 scalar Point3::smooth_point3_term(
     const Eigen::Matrix<scalar, n_verts, 3>& X,
-    Eigen::ConstRef<Eigen::RowVector3<scalar>> direc) const
+    Eigen::ConstRef<Eigen::RowVector3<scalar>> direc) const // NOLINT(readability-named-parameter)
 {
     const Eigen::RowVector3<scalar> dn = direc.normalized();
     scalar tangent_term(1.);
@@ -582,7 +587,8 @@ scalar Point3::smooth_point3_term(
         }
     }
 
-    weight = m_params.use_rest_shape_measure ? scalar(m_rest_weight) : weight / 3.;
+    weight =
+        m_params.use_rest_shape_measure ? scalar(m_rest_weight) : weight / 3.;
 
     if (!orientable || otypes.normal_type(0) == HeavisideType::ONE) {
         normal_term = scalar(1.);
diff --git a/src/ipc/smooth_contact/primitives/point3.hpp b/src/ipc/smooth_contact/primitives/point3.hpp
index d29b09d13..d64185475 100644
--- a/src/ipc/smooth_contact/primitives/point3.hpp
+++ b/src/ipc/smooth_contact/primitives/point3.hpp
@@ -96,7 +96,8 @@ class Point3 : public Primitive {
     Eigen::MatrixX3i faces;
     Eigen::MatrixX2i edges;
     bool orientable;
-    double m_rest_weight; ///< Rest-shape vertex area measure (quadrature weight, constant)
+    double m_rest_weight; ///< Rest-shape vertex area measure (quadrature
+                          ///< weight, constant)
 
     bool smooth_point3_term_type(
         Eigen::ConstRef<Eigen::MatrixX3d> X,

From 2322b666995fbf80ce58552a9a16c60997e18655 Mon Sep 17 00:00:00 2001
From: Uday Kusupati <kusupatiuday@gmail.com>
Date: Thu, 30 Apr 2026 20:46:15 -0400
Subject: [PATCH 5/6] change default use_rest_shape_measure to false

---
 src/ipc/smooth_contact/common.hpp            | 2 +-
 src/ipc/smooth_contact/primitives/point3.cpp | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/ipc/smooth_contact/common.hpp b/src/ipc/smooth_contact/common.hpp
index eac1f1115..53df77678 100644
--- a/src/ipc/smooth_contact/common.hpp
+++ b/src/ipc/smooth_contact/common.hpp
@@ -99,7 +99,7 @@ struct SmoothContactParameters {
     double beta_n = 0;
     int r = 2;
 
-    bool use_rest_shape_measure = true;
+    bool use_rest_shape_measure = false;
 
 private:
     double m_adaptive_dhat_ratio = 0.5;
diff --git a/src/ipc/smooth_contact/primitives/point3.cpp b/src/ipc/smooth_contact/primitives/point3.cpp
index 55a74c0ba..35bba5b42 100644
--- a/src/ipc/smooth_contact/primitives/point3.cpp
+++ b/src/ipc/smooth_contact/primitives/point3.cpp
@@ -566,7 +566,8 @@ HessianType<-1> Point3::smooth_point3_term_hessian(
 template <typename scalar, int n_verts>
 scalar Point3::smooth_point3_term(
     const Eigen::Matrix<scalar, n_verts, 3>& X,
-    Eigen::ConstRef<Eigen::RowVector3<scalar>> direc) const // NOLINT(readability-named-parameter)
+    // NOLINTNEXTLINE(readability-named-parameter)
+    Eigen::ConstRef<Eigen::RowVector3<scalar>> direc) const
 {
     const Eigen::RowVector3<scalar> dn = direc.normalized();
     scalar tangent_term(1.);

From de1134724874a5c909963d6ee048fc812de81418 Mon Sep 17 00:00:00 2001
From: Uday Kusupati <kusupatiuday@gmail.com>
Date: Thu, 30 Apr 2026 21:18:43 -0400
Subject: [PATCH 6/6] bug fix

---
 src/ipc/smooth_contact/primitives/point2.cpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/ipc/smooth_contact/primitives/point2.cpp b/src/ipc/smooth_contact/primitives/point2.cpp
index 0691f333d..c818d7bcb 100644
--- a/src/ipc/smooth_contact/primitives/point2.cpp
+++ b/src/ipc/smooth_contact/primitives/point2.cpp
@@ -189,7 +189,7 @@ VectorMax<double, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
         const T measure = m_params.use_rest_shape_measure
             ? T(m_rest_measure)
-            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm())
+            : ((X.row(2) - X.row(1)).norm() + (X.row(3) - X.row(1)).norm())
                 / 2.0;
         return smooth_point2_term<T>(
                    X.row(1), X.row(0), X.row(2), X.row(3), m_params, orientable,
@@ -203,7 +203,7 @@ VectorMax<double, Point2::MAX_SIZE + Point2::DIM> Point2::grad(
         Eigen::Matrix<T, 3, DIM> X = slice_positions<T, 3, DIM>(tmp);
         const T measure = m_params.use_rest_shape_measure
             ? T(m_rest_measure)
-            : (X.row(2) - X.row(0)).norm();
+            : (X.row(2) - X.row(1)).norm();
         return smooth_point2_term_one_side<T>(
                    X.row(1), X.row(0), X.row(2), m_params, measure)
             .grad;
@@ -228,7 +228,7 @@ Point2::hessian(
         Eigen::Matrix<T, 4, DIM> X = slice_positions<T, 4, DIM>(tmp);
         const T measure = m_params.use_rest_shape_measure
             ? T(m_rest_measure)
-            : ((X.row(2) - X.row(0)).norm() + (X.row(3) - X.row(0)).norm())
+            : ((X.row(2) - X.row(1)).norm() + (X.row(3) - X.row(1)).norm())
                 / 2.0;
         return smooth_point2_term<T>(
                    X.row(1), X.row(0), X.row(2), X.row(3), m_params, orientable,
@@ -242,7 +242,7 @@ Point2::hessian(
         Eigen::Matrix<T, 3, DIM> X = slice_positions<T, 3, DIM>(tmp);
         const T measure = m_params.use_rest_shape_measure
             ? T(m_rest_measure)
-            : (X.row(2) - X.row(0)).norm();
+            : (X.row(2) - X.row(1)).norm();
         return smooth_point2_term_one_side<T>(
                    X.row(1), X.row(0), X.row(2), m_params, measure)
             .Hess;