Individual convergent formulation flags

.github/workflows/clang-format-check.yml

@@ -29,5 +29,5 @@ jobs:
     - name: clang-format style check
       uses: jidicula/[email protected]
       with:
-        clang-format-version: '17'
+        clang-format-version: '18'
         check-path: ${{ matrix.path }}

docs/source/tutorial/convergent.rst

@@ -5,26 +5,34 @@ Convergent Formulation
 
 In addition to the original implementation of :cite:t:`Li2020IPC`, we also implement the convergent formulation of :cite:t:`Li2023Convergent`.
 
-To enable the convergent formulation, we need to set ``use_convergent_formulation`` before building ``Collisions``:
+Fully enabling the convergent formulation requires to set three flags: ``use_area_weighting`` and ``use_improved_max_approximator`` in ``Collisions`` (before calling ``build``) and ``use_physical_barrier`` in ``BarrierPotential``.
 
 .. md-tab-set::
 
     .. md-tab-item:: C++
 
         .. code-block:: c++
 
-            collisions.set_use_convergent_formulation(true);
+            collisions.set_use_area_weighting(true);
+            collisions.set_use_improved_max_approximator(true);
             collisions.build(collision_mesh, vertices, dhat);
 
+            barrier_potential.set_use_physical_barrier(true);
+            double b = barrier_potential(collisions, mesh, vertices);
+
     .. md-tab-item:: Python
 
         .. code-block:: python
 
-            collisions.use_convergent_formulation = True
+            collisions.use_area_weighting = True
+            collisions.use_improved_max_approximator = True
             collisions.build(collision_mesh, vertices, dhat)
 
+            barrier_potential.use_physical_barrier = True
+            b = barrier_potential(collisions, mesh, vertices);
+
 .. important::
-    The variable ``use_convergent_formulation`` should be set before calling ``Collisions::build`` for it to take effect. By default, it is ``false``.
+    The flags ``use_area_weighting`` and ``use_improved_max_approximator`` should be set before calling ``build`` for them to take effect. By default, they are ``false``.
 
 Technical Details
 -----------------
@@ -34,7 +42,7 @@ Technical Details
 In order to derive a convergent formulation, we first define a continuous form of our barrier potential :math:`P`. For a surface :math:`\mathcal{S}` embedded in 3D space, we parameterize the surfaces by common (possibly discontinuous) coordinates :math:`u \in \tilde{M} \subset \mathbb{R}^2`, so that :math:`\mathcal{S}(u)` traverses the material points across all surfaces contiguously. The total barrier potential is then
 
 .. math::
-   P(\mathcal{S})=\frac{1}{2} \int_{u \in \tilde{M}} \max _{v \in \tilde{M} \setminus{ }_r u} b(d(\mathcal{S}(u), \mathcal{S}(v)), \hat{d})~\mathrm{d} u,
+   P(\mathcal{S})=\frac{1}{2} \int_{u \in \tilde{M}}~\max_{v \in \tilde{M} \setminus{ }_r u} b(d(\mathcal{S}(u), \mathcal{S}(v)), \hat{d})~\mathrm{d} u,
 
 where we define the operator :math:`\setminus_r: \mathcal{P}(\mathbb{R}^2) \times \mathbb{R} \times \mathbb{R}^2 \mapsto \mathcal{P}(\mathbb{R}^2)` to be
 
@@ -57,6 +65,9 @@ Applying mesh vertices as nodes (and quadrature points), we numerically integrat
 
 where :math:`w_{\bar{x}}` are the quadrature weights, each given by one-third of the sum of the areas (in material space) of the boundary triangles incident to :math:`\bar{x}`.
 
+.. note::
+    The area weighted quadrature is enabled by setting ``use_area_weighting`` to ``true`` in ``Collisions``.
+
 We next need to smoothly approximate the max operator in the barrier potentials. However, common approaches such as an :math:`L^p`-norm or LogSumExp would decrease sparsity in subsequent numerical solves by increasing the stencil size per collision evaluation. We instead leverage the locality of our barrier function to approximate the max operator by removing duplicate distance pairs. Our resulting approximators for a triangulated surface is
 
 .. math::
@@ -67,6 +78,9 @@ We next need to smoothly approximate the max operator in the barrier potentials.
 
 where :math:`V_{\text{int}} \subseteq V` is the subset of internal surface nodes and :math:`E_{\text{int}} \subseteq E` is the subset of internal surface edges (i.e., edges incident to two triangles). For locally convex regions this estimator is tight while remaining smooth. In turn, for nonconvex regions, it improves over direct summation.
 
+.. note::
+    The improved max approximator is enabled by setting ``use_improved_max_approximator`` to ``true`` in ``Collisions``.
+
 The corresponding discrete barrier potential is then simply
 
 .. math::
@@ -101,6 +115,9 @@ The barrier stiffness (:math:`\kappa`) then has units of pressure (e.g., :math:`
 This implies we can get good solver convergence even when using a fixed :math:`\kappa` by setting it relative to the material's Young's modulus (:math:`\kappa = 0.1 E` works well in many examples).
 The intention is to treat the barrier as a thin elastic region around the mesh, and having consistent units makes it easier to pick the stiffness for this "material".
 
+.. note::
+    The physical barrier is enabled by setting ``use_physical_barrier`` to ``true`` in ``BarrierPotential``.
+
 .. _convergent-friction-formulation:
 
 Friction

python/src/barrier/barrier.cpp

@@ -12,6 +12,7 @@ class PyBarrier : public Barrier {
     double operator()(const double d, const double dhat) const override { PYBIND11_OVERRIDE_PURE_NAME(double, Barrier, "__call__", operator(), d, dhat); }
     double first_derivative(const double d, const double dhat) const override { PYBIND11_OVERRIDE_PURE(double, Barrier, first_derivative, d, dhat); }
     double second_derivative(const double d, const double dhat) const override { PYBIND11_OVERRIDE_PURE(double, Barrier, second_derivative, d, dhat); }
+    double units(const double dhat) const override { PYBIND11_OVERRIDE_PURE(double, Barrier, units, dhat); }
     // clang-format on
 };
 
@@ -25,7 +26,8 @@ void define_barrier(py::module_& m)
             py::arg("dhat"))
         .def(
             "second_derivative", &Barrier::second_derivative, py::arg("d"),
-            py::arg("dhat"));
+            py::arg("dhat"))
+        .def("units", &Barrier::units, py::arg("dhat"));
 
     py::class_<ClampedLogBarrier, Barrier, std::shared_ptr<ClampedLogBarrier>>(
         m, "ClampedLogBarrier")

python/src/collisions/collisions.cpp

@@ -138,15 +138,18 @@ void define_collisions(py::module_& m)
             "to_string", &Collisions::to_string, "", py::arg("mesh"),
             py::arg("vertices"))
         .def_property(
-            "use_convergent_formulation",
-            &Collisions::use_convergent_formulation,
-            &Collisions::set_use_convergent_formulation,
-            "If the collisions should use the convergent formulation.")
+            "use_area_weighting", &Collisions::use_area_weighting,
+            &Collisions::set_use_area_weighting,
+            "If the Collisions should use the convergent formulation.")
         .def_property(
-            "are_shape_derivatives_enabled",
-            &Collisions::are_shape_derivatives_enabled,
-            &Collisions::set_are_shape_derivatives_enabled,
-            "If the collisions are using the convergent formulation.")
+            "use_improved_max_approximator",
+            &Collisions::use_improved_max_approximator,
+            &Collisions::set_use_improved_max_approximator,
+            "If the Collisions should use the improved max approximator.")
+        .def_property(
+            "enable_shape_derivatives", &Collisions::enable_shape_derivatives,
+            &Collisions::set_enable_shape_derivatives,
+            "If the Collisions are using the convergent formulation.")
         .def(
             "to_string", &Collisions::to_string, py::arg("mesh"),
             py::arg("vertices"))

python/src/potentials/barrier_potential.cpp

@@ -9,24 +9,26 @@ void define_barrier_potential(py::module_& m)
 {
     py::class_<BarrierPotential>(m, "BarrierPotential")
         .def(
-            py::init<const double>(),
+            py::init<const double, const bool>(),
             R"ipc_Qu8mg5v7(
             Construct a barrier potential.
 
             Parameters:
                 dhat: The activation distance of the barrier.
             )ipc_Qu8mg5v7",
-            py::arg("dhat"))
+            py::arg("dhat"), py::arg("use_physical_barrier") = false)
         .def(
-            py::init<const std::shared_ptr<Barrier>, const double>(),
+            py::init<
+                const std::shared_ptr<Barrier>, const double, const bool>(),
             R"ipc_Qu8mg5v7(
             Construct a barrier potential.
 
             Parameters:
                 barrier: The barrier function.
                 dhat: The activation distance of the barrier.
             )ipc_Qu8mg5v7",
-            py::arg("barrier"), py::arg("dhat"))
+            py::arg("barrier"), py::arg("dhat"),
+            py::arg("use_physical_barrier") = false)
         .def(
             "__call__",
             py::overload_cast<

python/tests/test_ipc.py

@@ -16,12 +16,14 @@ def check_ipc_derivatives(broad_phase_method, use_convergent_formulation, mesh_n
         vertices = mesh.vertices(vertices)
 
     collisions = ipctk.Collisions()
-    collisions.use_convergent_formulation = use_convergent_formulation
+    collisions.use_area_weighting = use_convergent_formulation
+    collisions.use_improved_max_approximator = use_convergent_formulation
     collisions.build(mesh, vertices, dhat,
                      broad_phase_method=broad_phase_method)
     assert len(collisions) > 0
 
-    B = ipctk.BarrierPotential(dhat)
+    B = ipctk.BarrierPotential(
+        dhat, use_physical_barrier=use_convergent_formulation)
 
     grad_b = B.gradient(collisions, mesh, vertices)
     fgrad_b = utils.finite_gradient(

src/ipc/barrier/adaptive_stiffness.cpp

@@ -48,8 +48,7 @@ double initial_barrier_stiffness(
         assert(std::isfinite(kappa));
     }
 
-    return std::min(
-        max_barrier_stiffness, std::max(min_barrier_stiffness, kappa));
+    return std::clamp(kappa, min_barrier_stiffness, max_barrier_stiffness);
 }
 
 // Adaptive κ

src/ipc/barrier/barrier.hpp

@@ -12,11 +12,30 @@ class Barrier {
     Barrier() = default;
     virtual ~Barrier() = default;
 
+    /// @brief Evaluate the barrier function.
+    /// @param d Distance.
+    /// @param dhat Activation distance of the barrier.
+    /// @return The value of the barrier function at d.
     virtual double operator()(const double d, const double dhat) const = 0;
+
+    /// @brief Evaluate the first derivative of the barrier function wrt d.
+    /// @param d Distance.
+    /// @param dhat Activation distance of the barrier.
+    /// @retur The value of the first derivative of the barrier function at d.
     virtual double
     first_derivative(const double d, const double dhat) const = 0;
+
+    /// @brief Evaluate the second derivative of the barrier function wrt d.
+    /// @param d Distance.
+    /// @param dhat Activation distance of the barrier.
+    /// @return The value of the second derivative of the barrier function at d.
     virtual double
     second_derivative(const double d, const double dhat) const = 0;
+
+    /// @brief Get the units of the barrier function.
+    /// @param dhat The activation distance of the barrier.
+    /// @return The units of the barrier function.
+    virtual double units(const double dhat) const = 0;
 };
 
 // ============================================================================
@@ -106,6 +125,15 @@ class ClampedLogBarrier : public Barrier {
     {
         return barrier_second_derivative(d, dhat);
     }
+
+    /// @brief Get the units of the barrier function.
+    /// @param dhat The activation distance of the barrier.
+    /// @return The units of the barrier function.
+    double units(const double dhat) const override
+    {
+        // (d - d̂)² = d̂² (d/d̂ - 1)²
+        return dhat * dhat;
+    }
 };
 
 } // namespace ipc

src/ipc/broad_phase/aabb.hpp

@@ -20,8 +20,8 @@ class AABB {
 
     AABB(const AABB& aabb1, const AABB& aabb2, const AABB& aabb3)
         : AABB(
-            aabb1.min.min(aabb2.min).min(aabb3.min),
-            aabb1.max.max(aabb2.max).max(aabb3.max))
+              aabb1.min.min(aabb2.min).min(aabb3.min),
+              aabb1.max.max(aabb2.max).max(aabb3.max))
     {
     }
 

src/ipc/collision_mesh.cpp

@@ -14,11 +14,11 @@ CollisionMesh::CollisionMesh(
     const Eigen::MatrixXi& faces,
     const Eigen::SparseMatrix<double>& displacement_map)
     : CollisionMesh(
-        std::vector<bool>(rest_positions.rows(), true),
-        rest_positions,
-        edges,
-        faces,
-        displacement_map)
+          std::vector<bool>(rest_positions.rows(), true),
+          rest_positions,
+          edges,
+          faces,
+          displacement_map)
 {
 }
 

src/ipc/collisions/collisions.cpp

@@ -173,7 +173,7 @@ void Collisions::build(
     };
 
     tbb::enumerable_thread_specific<CollisionsBuilder> storage(
-        use_convergent_formulation(), are_shape_derivatives_enabled());
+        use_area_weighting(), enable_shape_derivatives());
 
     tbb::parallel_for(
         tbb::blocked_range<size_t>(size_t(0), candidates.vv_candidates.size()),
@@ -207,7 +207,7 @@ void Collisions::build(
                 r.end());
         });
 
-    if (use_convergent_formulation()) {
+    if (use_improved_max_approximator()) {
         if (candidates.ev_candidates.size() > 0) {
             // Convert edge-vertex to vertex-vertex
             const std::vector<VertexVertexCandidate> vv_candidates =
@@ -276,50 +276,41 @@ void Collisions::build(
         Collision& collision = (*this)[ci];
         collision.dmin = dmin;
     }
+}
 
-    if (use_convergent_formulation()) {
-        // NOTE: When using the convergent formulation we want the barrier to
-        // have units of Pa⋅m, so κ gets units of Pa and the barrier function
-        // should have units of m. See notebooks/physical_barrier.ipynb for more
-        // details.
-        const double barrier_to_physical_barrier_divisor =
-            dhat * std::pow(dhat + 2 * dmin, 2);
-
-        for (size_t ci = 0; ci < size(); ci++) {
-            Collision& collision = (*this)[ci];
-            collision.weight /= barrier_to_physical_barrier_divisor;
-            if (are_shape_derivatives_enabled()) {
-                collision.weight_gradient /=
-                    barrier_to_physical_barrier_divisor;
-            }
-        }
+void Collisions::set_use_area_weighting(const bool use_area_weighting)
+{
+    if (!empty() && use_area_weighting != m_use_area_weighting) {
+        logger().warn("Setting use_area_weighting after building collisions. "
+                      "Re-build collisions for this to have an effect.");
     }
+
+    m_use_area_weighting = use_area_weighting;
 }
 
-void Collisions::set_use_convergent_formulation(
-    const bool use_convergent_formulation)
+void Collisions::set_use_improved_max_approximator(
+    const bool use_improved_max_approximator)
 {
     if (!empty()
-        && use_convergent_formulation != m_use_convergent_formulation) {
+        && use_improved_max_approximator != m_use_improved_max_approximator) {
         logger().warn(
-            "Setting use_convergent_formulation after building collisions. "
+            "Setting use_improved_max_approximator after building collisions. "
             "Re-build collisions for this to have an effect.");
     }
 
-    m_use_convergent_formulation = use_convergent_formulation;
+    m_use_improved_max_approximator = use_improved_max_approximator;
 }
 
-void Collisions::set_are_shape_derivatives_enabled(
-    const bool are_shape_derivatives_enabled)
+void Collisions::set_enable_shape_derivatives(
+    const bool enable_shape_derivatives)
 {
-    if (!empty()
-        && are_shape_derivatives_enabled != m_are_shape_derivatives_enabled) {
+    if (!empty() && enable_shape_derivatives != m_enable_shape_derivatives) {
         logger().warn(
             "Setting enable_shape_derivatives after building collisions. "
             "Re-build collisions for this to have an effect.");
     }
 
-    m_are_shape_derivatives_enabled = are_shape_derivatives_enabled;
+    m_enable_shape_derivatives = enable_shape_derivatives;
 }
 
 // ============================================================================