Allow user-provided vertex and edge areas in CollisionMesh

CLAUDE.md

@@ -76,6 +76,18 @@ cd build/test && ctest --verbose -R "test_name_pattern"
 ./build/test/ipc_toolkit_tests "[tag]"
 ```
 
+### Python Tests
+
+Python binding tests live in `python/tests/` and use **nose2** (listed in `python/tests/requirements.txt`). They are plain `test_*` functions with bare `assert`s, not pytest fixtures or `unittest` classes. `find_ipctk` locates a built `ipctk` module, so build the bindings first (`cmake --build --preset=python`).
+
+```bash
+# Run all Python tests
+cd python/tests && python -m nose2 -v
+
+# Run one test file
+cd python/tests && python -m nose2 -v test_collision_mesh
+```
+
 ### Benchmarks
 
 Benchmarks are hidden Catch2 cases tagged `[!benchmark]`, so they don't run by

python/src/barrier/adaptive_stiffness.cpp

@@ -30,7 +30,7 @@ void define_adaptive_stiffness(py::module_& m)
             bbox_diagonal: Length of the diagonal of the bounding box of the scene.
             barrier: Barrier function.
             dhat: Activation distance of the barrier.
-            average_mass: Average mass of all bodies.
+            average_mass: Average mass of all bodies. With weighted collisions (e.g., area weighting), pass average_mass / mean(weight) so the adaptive range accounts for the weights.
             grad_energy: Gradient of the elasticity energy function.
             grad_barrier: Gradient of the barrier potential.
             min_barrier_stiffness_scale: Scale used to premultiply the minimum barrier stiffness.

python/src/collision_mesh.cpp

@@ -455,9 +455,14 @@ void define_collision_mesh(py::module_& m)
                 Barycentric area of vertex vi.
             )ipc_Qu8mg5v7",
             "vi"_a)
-        .def_property_readonly(
+        .def_property(
             "vertex_areas", &CollisionMesh::vertex_areas,
-            "Get the barycentric area of the vertices.")
+            &CollisionMesh::set_vertex_areas,
+            R"ipc_Qu8mg5v7(
+            The vertex areas used as collision weights (barycentric by default).
+
+            Setting custom areas (e.g., 2πR·ds for axisymmetric models) does not support shape derivatives and must be done before init_area_jacobians(). Collisions built earlier keep their old weights.
+            )ipc_Qu8mg5v7")
         .def(
             "vertex_area_gradient",
             [](const CollisionMesh& self,
@@ -486,9 +491,14 @@ void define_collision_mesh(py::module_& m)
                 Barycentric area of edge ei.
             )ipc_Qu8mg5v7",
             "ei"_a)
-        .def(
+        .def_property(
             "edge_areas", &CollisionMesh::edge_areas,
-            "Get the barycentric area of the edges.")
+            &CollisionMesh::set_edge_areas,
+            R"ipc_Qu8mg5v7(
+            The edge areas used as collision weights (barycentric by default).
+
+            See vertex_areas for the caveats on custom areas.
+            )ipc_Qu8mg5v7")
         .def(
             "edge_area_gradient",
             [](const CollisionMesh& self,

python/tests/test_collision_mesh.py

@@ -3,7 +3,12 @@
 import find_ipctk
 import numpy as np
 import scipy
-from ipctk import CollisionMesh, make_sparse_filter, make_vertex_patches_filter
+from ipctk import (
+    CollisionMesh,
+    NormalCollisions,
+    make_sparse_filter,
+    make_vertex_patches_filter,
+)
 
 
 def test_collision_mesh():
@@ -158,3 +163,103 @@ def test_face_normals_2d_raises():
         assert False
     except ValueError:
         pass
+
+
+def _two_segments():
+    V = np.array([[0, 0], [1, 0], [0, 0.1], [1, 0.1]], dtype=float)
+    E = np.array([[0, 1], [2, 3]], dtype=int)
+    return V, E
+
+
+def test_custom_areas_become_collision_weights():
+    V, E = _two_segments()
+    scale = 2 * np.pi * 3.0
+
+    default_mesh = CollisionMesh(V, E)
+    default_collisions = NormalCollisions()
+    default_collisions.use_area_weighting = True
+    default_collisions.build(default_mesh, V, dhat=0.2)
+
+    mesh = CollisionMesh(V, E)
+    mesh.vertex_areas = scale * mesh.vertex_areas
+    mesh.edge_areas = scale * mesh.edge_areas
+    assert np.allclose(mesh.vertex_areas, scale * default_mesh.vertex_areas)
+    assert np.allclose(mesh.edge_areas, scale * default_mesh.edge_areas)
+
+    collisions = NormalCollisions()
+    collisions.use_area_weighting = True
+    collisions.build(mesh, V, dhat=0.2)
+
+    assert len(collisions) > 0
+    assert len(collisions) == len(default_collisions)
+    for i in range(len(collisions)):
+        assert np.isclose(
+            collisions[i].weight, scale * default_collisions[i].weight)
+
+
+def test_invalid_custom_areas_are_rejected():
+    V, E = _two_segments()
+    mesh = CollisionMesh(V, E)
+    original = mesh.vertex_areas.copy()
+
+    for bad in [np.ones(3), np.full(4, -1.0), np.full(4, np.nan), np.full(4, np.inf)]:
+        try:
+            mesh.vertex_areas = bad
+            assert False
+        except (ValueError, RuntimeError):
+            pass
+    try:
+        mesh.edge_areas = np.ones(3)
+        assert False
+    except (ValueError, RuntimeError):
+        pass
+
+    assert np.array_equal(mesh.vertex_areas, original)
+
+
+def test_custom_areas_and_area_jacobians_are_exclusive():
+    V, E = _two_segments()
+
+    # Setting after initializing the Jacobians
+    mesh = CollisionMesh(V, E)
+    mesh.init_area_jacobians()
+    try:
+        mesh.vertex_areas = 2 * mesh.vertex_areas
+        assert False
+    except (ValueError, RuntimeError):
+        pass
+
+    # Initializing after setting
+    mesh = CollisionMesh(V, E)
+    mesh.vertex_areas = 2 * mesh.vertex_areas
+    try:
+        mesh.init_area_jacobians()
+        assert False
+    except (ValueError, RuntimeError):
+        pass
+
+
+def test_area_getters_are_read_only_views():
+    # In-place edits must fail loudly instead of bypassing the validating setter.
+    V, E = _two_segments()
+    mesh = CollisionMesh(V, E)
+    try:
+        mesh.vertex_areas[0] = 5.0
+        assert False
+    except ValueError:
+        pass
+
+
+def test_zero_custom_areas_drop_their_collisions():
+    # E.g., axisymmetric primitives on the symmetry axis have no ring area.
+    V, E = _two_segments()
+    mesh = CollisionMesh(V, E)
+    mesh.vertex_areas = np.zeros(4)
+    mesh.edge_areas = np.zeros(2)
+
+    collisions = NormalCollisions()
+    collisions.use_area_weighting = True
+    collisions.build(mesh, V, dhat=0.2)
+
+    # Zero-weight collisions are dropped, so no contact force is produced.
+    assert len(collisions) == 0

src/ipc/barrier/adaptive_stiffness.hpp

@@ -11,6 +11,7 @@
 namespace ipc {
 
 /// @brief Compute an inital barrier stiffness using the barrier potential gradient.
+/// @note This heuristic was designed for unweighted collisions. With area weighting and a physical barrier, the stiffness is a material parameter in the same units as Young's modulus, so choosing it directly relative to the elastic stiffness is usually more appropriate than adapting it. If you still use this function with weighted collisions, pass average_mass / mean(weight) as average_mass so the adaptive range accounts for the weights.
 /// @param[in] bbox_diagonal Length of the diagonal of the bounding box of the scene.
 /// @param[in] barrier Barrier function.
 /// @param[in] dhat Activation distance of the barrier.

src/ipc/collision_mesh.cpp

@@ -376,8 +376,55 @@ void CollisionMesh::init_areas()
     }
 }
 
+namespace {
+    void check_areas(
+        Eigen::ConstRef<Eigen::VectorXd> areas,
+        const Eigen::Index expected_size,
+        const char* name)
+    {
+        if (areas.size() != expected_size) {
+            throw std::invalid_argument(
+                fmt::format(
+                    "{} has size {} but expected {}", name, areas.size(),
+                    expected_size));
+        }
+        if (!areas.allFinite() || (areas.array() < 0).any()) {
+            throw std::invalid_argument(
+                fmt::format("{} must be finite and non-negative", name));
+        }
+    }
+} // namespace
+
+void CollisionMesh::set_vertex_areas(
+    Eigen::ConstRef<Eigen::VectorXd> vertex_areas)
+{
+    check_areas(vertex_areas, num_vertices(), "vertex_areas");
+    if (are_area_jacobians_initialized()) {
+        throw std::runtime_error(
+            "Cannot set custom vertex areas after init_area_jacobians(): the area Jacobians would no longer match.");
+    }
+    m_vertex_areas = vertex_areas;
+    m_has_custom_areas = true;
+}
+
+void CollisionMesh::set_edge_areas(Eigen::ConstRef<Eigen::VectorXd> edge_areas)
+{
+    check_areas(edge_areas, num_edges(), "edge_areas");
+    if (are_area_jacobians_initialized()) {
+        throw std::runtime_error(
+            "Cannot set custom edge areas after init_area_jacobians(): the area Jacobians would no longer match.");
+    }
+    m_edge_areas = edge_areas;
+    m_has_custom_areas = true;
+}
+
 void CollisionMesh::init_area_jacobians()
 {
+    if (m_has_custom_areas) {
+        throw std::runtime_error(
+            "Cannot initialize area Jacobians after setting custom areas: shape derivatives are not supported with custom areas.");
+    }
+
     std::vector<bool> was_vertex_visited(num_vertices(), false);
     std::vector<bool> was_edge_visited(num_edges(), false);
 

src/ipc/collision_mesh.hpp

@@ -288,6 +288,12 @@ class CollisionMesh {
     /// @brief Get the barycentric area of the vertices.
     const Eigen::VectorXd& vertex_areas() const { return m_vertex_areas; }
 
+    /// @brief Override the vertex areas used as collision weights.
+    /// @note Use this when the quadrature weight of a vertex is not its barycentric area (e.g., axisymmetric models where each vertex stands for a ring of contact). Collisions built before this call keep their old weights.
+    /// @note Custom areas do not support shape derivatives: this throws if the area Jacobians were already initialized, and init_area_jacobians() throws afterwards.
+    /// @param vertex_areas New vertex areas (size num_vertices()). Must be finite and non-negative (zero is allowed, e.g., for axisymmetric primitives on the symmetry axis).
+    void set_vertex_areas(Eigen::ConstRef<Eigen::VectorXd> vertex_areas);
+
     /// @brief Get the gradient of the barycentric area of a vertex wrt the rest positions of all points.
     /// @param vi Vertex ID.
     /// @return Gradient of the barycentric area of vertex vi wrt the rest positions of all points.
@@ -309,6 +315,11 @@ class CollisionMesh {
     /// @brief Get the barycentric area of the edges.
     const Eigen::VectorXd& edge_areas() const { return m_edge_areas; }
 
+    /// @brief Override the edge areas used as collision weights.
+    /// @note See set_vertex_areas() for the caveats on stale collisions and shape derivatives.
+    /// @param edge_areas New edge areas (size num_edges()). Must be finite and non-negative (zero is allowed, e.g., for axisymmetric primitives on the symmetry axis).
+    void set_edge_areas(Eigen::ConstRef<Eigen::VectorXd> edge_areas);
+
     /// @brief Get the gradient of the barycentric area of an edge wrt the rest positions of all points.
     /// @param ei Edge ID.
     /// @return Gradient of the barycentric area of edge ei wrt the rest positions of all points.
@@ -451,6 +462,8 @@ class CollisionMesh {
     /// @brief Edge areas
     /// 3D: 1/3 sum of area of connected triangles
     Eigen::VectorXd m_edge_areas;
+    /// @brief Were the areas overridden by the user (so the Jacobians no longer match)?
+    bool m_has_custom_areas = false;
 
     // Stored as a std::vector so it is easier to access the rows directly.
     /// @brief The rows of the Jacobian of the vertex areas vector.

tests/src/tests/test_collision_mesh.cpp

@@ -3,6 +3,8 @@
 #include <catch2/catch_approx.hpp>
 #include <catch2/generators/catch_generators.hpp>
 #include <ipc/collision_mesh.hpp>
+#include <ipc/collisions/normal/normal_collisions.hpp>
+#include <ipc/potentials/barrier_potential.hpp>
 #include <ipc/config.hpp>
 
 using namespace ipc;
@@ -251,3 +253,82 @@ TEST_CASE("face_normals", "[collision_mesh][face_normals]")
         }
     }
 }
+
+TEST_CASE("User-provided collision mesh areas", "[collision_mesh][areas]")
+{
+    // Two parallel horizontal segments, close enough to collide.
+    Eigen::MatrixXd V(4, 2);
+    V << 0, 0, 1, 0, 0, 0.1, 1, 0.1;
+    Eigen::MatrixXi E(2, 2);
+    E << 0, 1, 2, 3;
+
+    CollisionMesh mesh(V, E);
+    const Eigen::VectorXd default_areas = mesh.vertex_areas();
+
+    SECTION("Custom areas become the collision weights")
+    {
+        const double scale = 2 * M_PI * 3.0;
+        mesh.set_vertex_areas(scale * default_areas);
+        mesh.set_edge_areas(scale * mesh.edge_areas());
+
+        NormalCollisions collisions;
+        collisions.set_use_area_weighting(true);
+        collisions.build(mesh, V, /*dhat=*/0.2);
+        REQUIRE(collisions.size() > 0);
+
+        NormalCollisions default_collisions;
+        default_collisions.set_use_area_weighting(true);
+        const CollisionMesh default_mesh(V, E);
+        default_collisions.build(default_mesh, V, /*dhat=*/0.2);
+        REQUIRE(default_collisions.size() == collisions.size());
+
+        for (size_t i = 0; i < collisions.size(); ++i) {
+            CHECK(
+                collisions[i].weight
+                == Catch::Approx(scale * default_collisions[i].weight));
+        }
+    }
+
+    SECTION("Zero areas drop their collisions and exert no contact force")
+    {
+        // E.g., axisymmetric primitives on the symmetry axis have no ring area.
+        mesh.set_vertex_areas(Eigen::VectorXd::Zero(mesh.num_vertices()));
+        mesh.set_edge_areas(Eigen::VectorXd::Zero(mesh.num_edges()));
+
+        NormalCollisions collisions;
+        collisions.set_use_area_weighting(true);
+        collisions.build(mesh, V, /*dhat=*/0.2);
+        CHECK(collisions.empty());
+
+        const BarrierPotential potential(/*dhat=*/0.2, /*kappa=*/1.0);
+        CHECK(potential(collisions, mesh, V) == 0);
+        const Eigen::VectorXd grad = potential.gradient(collisions, mesh, V);
+        CHECK(grad.allFinite());
+        CHECK(grad.isZero());
+    }
+
+    SECTION("Invalid areas are rejected")
+    {
+        CHECK_THROWS(mesh.set_vertex_areas(Eigen::VectorXd::Ones(3)));
+        CHECK_THROWS(mesh.set_vertex_areas(Eigen::VectorXd::Constant(4, -1)));
+        CHECK_THROWS(mesh.set_edge_areas(Eigen::VectorXd::Ones(3)));
+        CHECK_THROWS(
+            mesh.set_edge_areas(Eigen::VectorXd::Constant(2, std::nan(""))));
+        CHECK(mesh.vertex_areas() == default_areas);
+    }
+
+    SECTION("Custom areas and area Jacobians are mutually exclusive")
+    {
+        SECTION("Setting after initializing Jacobians")
+        {
+            mesh.init_area_jacobians();
+            CHECK_THROWS(mesh.set_vertex_areas(default_areas));
+            CHECK_THROWS(mesh.set_edge_areas(mesh.edge_areas()));
+        }
+        SECTION("Initializing after setting")
+        {
+            mesh.set_vertex_areas(2 * default_areas);
+            CHECK_THROWS(mesh.init_area_jacobians());
+        }
+    }
+}