Divide and Truncate

docs/source/references.bib

@@ -108,6 +108,15 @@ @article{Huang2024GIPC
   pages   = {1--18},
   doi     = {10.1145/3643028}
 }
+@article{Chen2026DivideAndTruncate,
+  title         = {Divide and Truncate: A Penetration and Inversion Free Framework for Coupled Multi-physics Systems},
+  author        = {Anka He Chen and Jerry Hsu and Youssef Ayman and Miles Macklin},
+  year          = 2026,
+  eprint        = {2604.15513},
+  archivePrefix = {arXiv},
+  primaryClass  = {cs.GR},
+  note          = {\url{https://arxiv.org/abs/2604.15513}},
+}
 @article{Chen2025Offset,
   title   = {Offset {Geometric} {Contact}},
   author  = {Chen, Anka He and Hsu, Jerry and Liu, Ziheng and Macklin, Miles and Yang, Yin and Yuksel, Cem},

docs/source/tutorials/ogc.rst

@@ -383,4 +383,138 @@ The ``warm_start_time_step`` function includes an adaptive heuristic. During the
 2.  If they are outside, it effectively projects them to the boundary (using the scaling method described above).
 3.  **Crucially**, if the number of vertices hitting the boundary exceeds ``update_threshold`` immediately during this warm start, it triggers a second ``update()`` call.
 
-This prevents the solver from starting with a trust region that is already too restrictive for the bulk of the mesh, potentially saving solver iterations later.
+This prevents the solver from starting with a trust region that is already too restrictive for the bulk of the mesh, potentially saving solver iterations later.
+
+Planar-DAT (Divide and Truncate)
+---------------------------------
+
+The isotropic ``filter_step`` clips every vertex displacement to a sphere, which restricts motion in *all* directions — even tangential sliding or outward motion away from a contact. In dense or cluttered contact scenarios this causes **artificial damping** and **deadlock**: vertices cannot slide past each other even when no penetration would occur.
+
+**Planar-DAT** :cite:`Chen2026DivideAndTruncate` replaces the sphere with a *direction-aware division plane* computed per collision pair. Only the displacement component that crosses the plane (i.e., moves toward the opposing primitive) is truncated; tangential and separating motion passes through unconstrained. This typically allows steps 2× larger or more, eliminating damping artifacts.
+
+How the Division Plane is Computed
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+For each active collision pair the algorithm:
+
+1. Finds the closest points between the two primitives.
+2. Constructs a unit normal :math:`\mathbf{n}` pointing from one primitive toward the other.
+3. Distributes the available gap between the two primitives in proportion to their approach speeds, yielding a **division point** :math:`\mathbf{p}` on the segment connecting the closest points.
+4. For every vertex of both primitives, finds the ray–plane intersection time :math:`t_i` such that :math:`\mathbf{x}_u + t_i \, \delta\mathbf{x}_u` lies on the plane through :math:`\mathbf{p}` with normal :math:`\mathbf{n}`. If :math:`t_i \in [0, 1/\gamma_r)`, the vertex is truncated to :math:`\gamma_r \, t_i`; otherwise it is left unchanged.
+
+For vertices not covered by an active collision pair, an isotropic fallback clamps motion to a safe-zone sphere centered at ``trust_region_centers`` with radius ``trust_region_inflation_radius``. This fallback is distinct from the isotropic clamp in ``filter_step``, which uses the per-vertex radii ``trust_region_radii(i)``, but it serves the same purpose of ensuring that no uncovered vertex can escape the safe zone.
+
+Using ``planar_filter_step``
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``planar_filter_step`` is a *drop-in replacement* for ``filter_step`` inside the solver loop. It uses the collision candidates stored in ``trust_region.candidates`` by the last call to ``update`` or ``update_if_needed``. The relaxation ratio is controlled via the ``relaxed_radius_scaling`` member of the ``TrustRegion`` object (default 0.9).
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            // Inside the solver loop, replace:
+            //   trust_region.filter_step(mesh, x, dx);
+            // with:
+            trust_region.planar_filter_step(mesh, x, dx);
+
+            // Optionally tune the relaxation ratio via the struct member:
+            // trust_region.relaxed_radius_scaling = 0.9; // default
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            # Inside the solver loop, replace:
+            #   trust_region.filter_step(mesh, x, dx)
+            # with:
+            trust_region.planar_filter_step(mesh, x, dx)
+
+            # Optionally tune the relaxation ratio via the struct member:
+            # trust_region.relaxed_radius_scaling = 0.9  # default
+
+Full Optimization Loop with Planar-DAT
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. md-tab-set::
+
+    .. md-tab-item:: C++
+
+        .. code-block:: c++
+
+            // 1. Initial setup (same as isotropic OGC)
+            ipc::NormalCollisions collisions;
+            collisions.set_collision_set_type(ipc::NormalCollisions::CollisionSetType::OGC);
+            ipc::ogc::TrustRegion trust_region(dhat);
+
+            // 2. Warm start (Predict & Initialize)
+            trust_region.warm_start_time_step(
+                mesh, x, pred_x, collisions, dhat);
+
+            // 3. Solver Loop
+            for (int i = 0; i < max_iterations; ++i) {
+                // Update trust region centers/radii if needed
+                trust_region.update_if_needed(mesh, x, collisions, dhat);
+
+                // Compute search direction (Solver specific)
+                Eigen::MatrixXd dx = compute_search_direction(x, ...);
+
+                // Filter step using Planar-DAT (direction-aware truncation)
+                trust_region.planar_filter_step(mesh, x, dx);
+
+                // Update positions
+                x += dx;
+
+                // Check convergence...
+            }
+
+    .. md-tab-item:: Python
+
+        .. code-block:: python
+
+            # 1. Initial setup (same as isotropic OGC)
+            collisions = ipctk.NormalCollisions()
+            collisions.collision_set_type = ipctk.NormalCollisions.CollisionSetType.OGC
+            trust_region = ipctk.ogc.TrustRegion(dhat)
+
+            # 2. Warm start (Predict & Initialize)
+            trust_region.warm_start_time_step(
+                mesh, x, pred_x, collisions, dhat)
+
+            # 3. Solver Loop
+            for i in range(max_iterations):
+                # Update trust region centers/radii if needed
+                trust_region.update_if_needed(mesh, x, collisions, dhat)
+
+                # Compute search direction (Solver specific)
+                dx = compute_search_direction(x, ...)
+
+                # Filter step using Planar-DAT (direction-aware truncation)
+                trust_region.planar_filter_step(mesh, x, dx)
+
+                # Update positions
+                x += dx
+
+                # Check convergence...
+
+.. note::
+    Like ``filter_step``, ``planar_filter_step`` also sets
+    ``should_update_trust_region`` when a critical mass of vertices are
+    restricted by the trust region. Call ``update_if_needed`` at the top of
+    each solver iteration to refresh the collision set when needed.
+
+When to Use Planar-DAT vs. Isotropic
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
++-----------------------------------+-----------------------------------------------+
+| Use ``filter_step`` (isotropic)   | Use ``planar_filter_step`` (Planar-DAT)       |
++===================================+===============================================+
+| Simple scenes, sparse contact     | Dense contact, many simultaneous constraints  |
++-----------------------------------+-----------------------------------------------+
+| Stepping-stone implementation     | Coupled multi-physics simulations             |
++-----------------------------------+-----------------------------------------------+
+| When implementation simplicity    | When artificial damping or deadlock           |
+| is preferred                      | is observed with isotropic filtering          |
++-----------------------------------+-----------------------------------------------+

python/src/collisions/normal/normal_collisions.cpp

@@ -6,7 +6,7 @@
 using namespace ipc;
 
 template <typename collision_type, typename parent_type>
-void define_smooth_collision_template(py::module_& m, std::string name)
+void define_smooth_collision_template(py::module_& m, const std::string& name)
 {
     py::class_<collision_type, parent_type>(m, name.c_str())
         .def("name", &collision_type::name, "Get the type name of collision")
@@ -15,7 +15,7 @@ void define_smooth_collision_template(py::module_& m, std::string name)
             "Get the number of vertices");
 }
 
-void define_smooth_collisions(py::module_& m, std::string name)
+void define_smooth_collisions(py::module_& m, const std::string& name)
 {
     py::class_<SmoothCollisions>(m, name.c_str())
         .def(py::init())

python/src/ogc/trust_region.cpp

@@ -81,6 +81,26 @@ void define_trust_region(py::module_& m)
                 Sets should_update_trust_region to true if the trust region should be updated on the next iteration.
             )ipc_Qu8mg5v7",
             "mesh"_a, "x"_a, "dx"_a)
+        .def(
+            "planar_filter_step", &ogc::TrustRegion::planar_filter_step,
+            R"ipc_Qu8mg5v7(
+            Filter the optimization step dx using Planar-DAT (Divide and Truncate).
+
+            For each collision candidate (stored in ``candidates`` by the last
+            call to ``update``), computes a direction-aware division plane and
+            truncates only the component of displacement toward that plane.
+            This eliminates the artificial damping and deadlock of the isotropic
+            ``filter_step`` while retaining the penetration-free guarantee.
+
+            See "Divide and Truncate: A Penetration and Inversion Free Framework
+            for Coupled Multi-physics Systems" [ACM SIGGRAPH 2026].
+
+            Parameters:
+                mesh: The collision mesh.
+                x: Current vertex positions.
+                dx: Proposed vertex displacements (modified in-place).
+            )ipc_Qu8mg5v7",
+            "mesh"_a, "x"_a, "dx"_a)
         .def_readwrite(
             "trust_region_centers", &ogc::TrustRegion::trust_region_centers,
             "Centers of the trust regions for each vertex.")
@@ -118,5 +138,8 @@ void define_trust_region(py::module_& m)
         .def_readwrite(
             "should_update_trust_region",
             &ogc::TrustRegion::should_update_trust_region,
-            "If true, the trust region will be updated on the next call to ``update_if_needed``.");
+            "If true, the trust region will be updated on the next call to ``update_if_needed``.")
+        .def_readwrite(
+            "candidates", &ogc::TrustRegion::candidates,
+            "Collision candidates used for Planar-DAT. Updated by ``update()``.");
 }

src/ipc/broad_phase/aabb.cpp

@@ -2,7 +2,6 @@
 
 #include <ipc/utils/profiler.hpp>
 
-#include <tbb/blocked_range.h>
 #include <tbb/parallel_for.h>
 
 #include <cfenv>
@@ -64,15 +63,10 @@ void build_vertex_boxes(
 
     vertex_boxes.resize(vertices.rows());
 
-    tbb::parallel_for(
-        tbb::blocked_range<size_t>(0, vertices.rows()),
-        [&](const tbb::blocked_range<size_t>& r) {
-            for (size_t i = r.begin(); i < r.end(); i++) {
-                vertex_boxes[i] =
-                    AABB::from_point(vertices.row(i), inflation_radius);
-                vertex_boxes[i].vertex_ids = { { index_t(i), -1, -1 } };
-            }
-        });
+    tbb::parallel_for(0, index_t(vertices.rows()), [&](index_t i) {
+        vertex_boxes[i] = AABB::from_point(vertices.row(i), inflation_radius);
+        vertex_boxes[i].vertex_ids = { { i, -1, -1 } };
+    });
 }
 
 void build_vertex_boxes(
@@ -85,15 +79,11 @@ void build_vertex_boxes(
 
     vertex_boxes.resize(vertices_t0.rows());
 
-    tbb::parallel_for(
-        tbb::blocked_range<size_t>(0, vertices_t0.rows()),
-        [&](const tbb::blocked_range<size_t>& r) {
-            for (size_t i = r.begin(); i < r.end(); i++) {
-                vertex_boxes[i] = AABB::from_point(
-                    vertices_t0.row(i), vertices_t1.row(i), inflation_radius);
-                vertex_boxes[i].vertex_ids = { { index_t(i), -1, -1 } };
-            }
-        });
+    tbb::parallel_for(0, index_t(vertices_t0.rows()), [&](index_t i) {
+        vertex_boxes[i] = AABB::from_point(
+            vertices_t0.row(i), vertices_t1.row(i), inflation_radius);
+        vertex_boxes[i].vertex_ids = { { i, -1, -1 } };
+    });
 }
 
 void build_edge_boxes(
@@ -108,15 +98,11 @@ void build_edge_boxes(
         edge_boxes.resize(edges.rows());
     }
 
-    tbb::parallel_for(
-        tbb::blocked_range<size_t>(0, edges.rows()),
-        [&](const tbb::blocked_range<size_t>& r) {
-            for (size_t i = r.begin(); i < r.end(); i++) {
-                const int e0 = edges(i, 0), e1 = edges(i, 1);
-                edge_boxes[i] = AABB(vertex_boxes[e0], vertex_boxes[e1]);
-                edge_boxes[i].vertex_ids = { { e0, e1, -1 } };
-            }
-        });
+    tbb::parallel_for(0, index_t(edges.rows()), [&](index_t i) {
+        const int e0 = edges(i, 0), e1 = edges(i, 1);
+        edge_boxes[i] = AABB(vertex_boxes[e0], vertex_boxes[e1]);
+        edge_boxes[i].vertex_ids = { { e0, e1, -1 } };
+    });
 }
 
 void build_face_boxes(
@@ -131,16 +117,12 @@ void build_face_boxes(
         face_boxes.resize(faces.rows());
     }
 
-    tbb::parallel_for(
-        tbb::blocked_range<size_t>(0, faces.rows()),
-        [&](const tbb::blocked_range<size_t>& r) {
-            for (size_t i = r.begin(); i < r.end(); i++) {
-                const int f0 = faces(i, 0), f1 = faces(i, 1), f2 = faces(i, 2);
-                face_boxes[i] =
-                    AABB(vertex_boxes[f0], vertex_boxes[f1], vertex_boxes[f2]);
-                face_boxes[i].vertex_ids = { { f0, f1, f2 } };
-            }
-        });
+    tbb::parallel_for(0, index_t(faces.rows()), [&](index_t i) {
+        const int f0 = faces(i, 0), f1 = faces(i, 1), f2 = faces(i, 2);
+        face_boxes[i] =
+            AABB(vertex_boxes[f0], vertex_boxes[f1], vertex_boxes[f2]);
+        face_boxes[i].vertex_ids = { { f0, f1, f2 } };
+    });
 }
 
 } // namespace ipc

src/ipc/broad_phase/hash_grid.cpp

@@ -64,14 +64,9 @@ void HashGrid::insert_boxes(
 {
     tbb::enumerable_thread_specific<std::vector<HashItem>> storage;
 
-    tbb::parallel_for(
-        tbb::blocked_range<long>(0L, long(boxes.size())),
-        [&](const tbb::blocked_range<long>& range) {
-            auto& local_items = storage.local();
-            for (long i = range.begin(); i != range.end(); i++) {
-                insert_box(boxes[i], i, local_items);
-            }
-        });
+    tbb::parallel_for(0L, long(boxes.size()), [&](long i) {
+        insert_box(boxes[i], i, storage.local());
+    });
 
     merge_thread_local_vectors(storage, items);