diff --git a/docs/source/tutorials/adhesion.rst b/docs/source/tutorials/adhesion.rst
index 29e4901bf..1cca5597d 100644
--- a/docs/source/tutorials/adhesion.rst
+++ b/docs/source/tutorials/adhesion.rst
@@ -47,7 +47,7 @@ We can build a normal adhesion potential object and compute the adhesion potenti
                 const double Y = 1e3;
                 const double eps_c = 0.5;
 
-                const ipc::NormalAdhesionPotential A_n(dhat_p, dhat_a, Y, eps_c)
+                const ipc::NormalAdhesionPotential A_n(dhat_p, dhat_a, Y, eps_c);
                 double adhesion_potential = A_n(normal_collisions, collision_mesh, vertices);
 
     .. md-tab-item:: Python
@@ -132,7 +132,7 @@ We can build a tangential adhesion potential object and compute the adhesion pot
 
             eps_a = 0.01
             A_t = ipctk.TangentialAdhesionPotential(eps_a)
-            adhesion_potential = A_t(tangential_collisions, collision_mesh, displacement);
+            adhesion_potential = A_t(tangential_collisions, collision_mesh, displacement)
 
 Derivatives
 ^^^^^^^^^^^
@@ -177,8 +177,7 @@ Similar to the friction model (see, `friction <advanced_friction.html#separate-c
 
             ipc::TangentialCollisions tangential_collisions;
             tangential_collisions.build(
-                collision_mesh, vertices, collisions, B, barrier_stiffness,
-                mu_s, mu_k);
+                collision_mesh, vertices, collisions, A_n, mu_s, mu_k);
 
     .. md-tab-item:: Python
 
@@ -186,8 +185,7 @@ Similar to the friction model (see, `friction <advanced_friction.html#separate-c
 
             tangential_collisions = ipctk.TangentialCollisions()
             tangential_collisions.build(
-                collision_mesh, vertices, collisions, B, barrier_stiffness,
-                mu_s, mu_k)
+                collision_mesh, vertices, collisions, A_n, mu_s, mu_k)
 
 The tangential adhesion force mollifier is then multiplied by the smooth coefficient of tangential adhesion:
 
diff --git a/docs/source/tutorials/convergent.rst b/docs/source/tutorials/convergent.rst
index 5caccfbfd..162085c5a 100644
--- a/docs/source/tutorials/convergent.rst
+++ b/docs/source/tutorials/convergent.rst
@@ -5,7 +5,7 @@ Convergent Formulation
 
 In addition to the original implementation of :cite:t:`Li2020IPC`, we also implement the convergent formulation of :cite:t:`Li2023Convergent`.
 
-Fully enabling the convergent formulation requires to set three things: ``use_area_weighting`` and ``collision_set_type`` in ``Collisions`` (before calling ``build``) and ``use_physical_barrier`` in ``BarrierPotential``.
+Fully enabling the convergent formulation requires to set three things: ``use_area_weighting`` and ``collision_set_type`` in ``NormalCollisions`` (before calling ``build``) and ``use_physical_barrier`` in ``BarrierPotential``.
 
 .. md-tab-set::
 
@@ -13,29 +13,39 @@ Fully enabling the convergent formulation requires to set three things: ``use_ar
 
         .. code-block:: c++
 
+            ipc::NormalCollisions collisions;
             collisions.set_use_area_weighting(true);
             collisions.set_collision_set_type(
                 ipc::NormalCollisions::CollisionSetType::IMPROVED_MAX_APPROX);
             collisions.build(collision_mesh, vertices, dhat);
 
+            ipc::BarrierPotential barrier_potential(dhat, stiffness);
             barrier_potential.set_use_physical_barrier(true);
-            double b = barrier_potential(collisions, mesh, vertices);
+            double b = barrier_potential(collisions, collision_mesh, vertices);
 
     .. md-tab-item:: Python
 
         .. code-block:: python
 
+            collisions = ipctk.NormalCollisions()
             collisions.use_area_weighting = True
             collisions.collision_set_type = \
                 ipctk.NormalCollisions.CollisionSetType.IMPROVED_MAX_APPROX
             collisions.build(collision_mesh, vertices, dhat)
 
+            barrier_potential = ipctk.BarrierPotential(dhat, stiffness)
             barrier_potential.use_physical_barrier = True
-            b = barrier_potential(collisions, mesh, vertices);
+            b = barrier_potential(collisions, collision_mesh, vertices)
 
 .. important::
     The members ``use_area_weighting`` and ``collision_set_type`` should be set before calling ``build`` for them to take effect. By default, they are ``false`` and ``IPC``.
 
+.. tip::
+    ``use_physical_barrier`` can also be set when constructing the potential:
+    ``ipc::BarrierPotential(dhat, stiffness, /*use_physical_barrier=*/true)`` in
+    C++ or ``ipctk.BarrierPotential(dhat, stiffness, use_physical_barrier=True)``
+    in Python.
+
 Technical Details
 -----------------
 
@@ -68,7 +78,7 @@ 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}`.
 
 .. tip::
-    The area weighted quadrature is enabled by setting ``use_area_weighting`` to ``true`` in ``Collisions``.
+    The area weighted quadrature is enabled by setting ``use_area_weighting`` to ``true`` in ``NormalCollisions``.
 
 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
 
@@ -86,7 +96,7 @@ where :math:`V_{\text{int}} \subseteq V` is the subset of internal surface nodes
     Comparison of the improved max approximator (right) to and exact max (left) and the direct summation (middle). For obtuse angles, the improved max approximator is tight, while for acute angles it might overestimate the max in concave regions.
 
 .. tip::
-    The improved max approximator is enabled by setting ``collision_set_type`` to ``IMPROVED_MAX_APPROX`` in ``Collisions``.
+    The improved max approximator is enabled by setting ``collision_set_type`` to ``IMPROVED_MAX_APPROX`` in ``NormalCollisions``.
 
 The corresponding discrete barrier potential is then simply
 
diff --git a/docs/source/tutorials/gcp.rst b/docs/source/tutorials/gcp.rst
index a77b41540..c95d1bb41 100644
--- a/docs/source/tutorials/gcp.rst
+++ b/docs/source/tutorials/gcp.rst
@@ -392,7 +392,7 @@ As :math:`\alpha + \beta` decreases, the support of the Heaviside function shrin
 
 Additional internal parameters that may affect behavior:
 
-- **Adaptive dhat ratio** (default ``0.5``): Controls the ratio :math:`\epsilon(x) / d_c(x, f_0)` in the adaptive barrier localization. Can be set via ``SmoothContactParameters::set_adaptive_dhat_ratio()``.
+- **Adaptive dhat ratio** (default ``0.5``): Controls the ratio :math:`\epsilon(x) / d_c(x, f_0)` in the adaptive barrier localization. Set it via ``SmoothContactParameters::set_adaptive_dhat_ratio()`` in C++ or the ``SmoothContactParameters.adaptive_dhat_ratio`` property in Python.
 - **Element measure** :math:`L`: For vertices, this is set to the average edge length around the vertex; for edges, to the edge length; for faces, :math:`L` is not needed. This determines the strength of the potential for low-dimensional contact (edge–edge, edge–vertex, vertex–vertex).
 
 Friction
diff --git a/docs/source/tutorials/getting_started.rst b/docs/source/tutorials/getting_started.rst
index 02702ef62..efc83a52a 100644
--- a/docs/source/tutorials/getting_started.rst
+++ b/docs/source/tutorials/getting_started.rst
@@ -301,8 +301,7 @@ To compute the adaptive barrier stiffness, we can use two functions: ``initial_b
             bbox_diagonal = ipctk.world_bbox_diagonal_length(vertices)
 
             barrier_stiffness, max_barrier_stiffness = ipctk.initial_barrier_stiffness(
-                bbox_diagonal, B.barrier, dhat, avg_mass, grad_energy, grad_barrier,
-                max_barrier_stiffness)
+                bbox_diagonal, B.barrier, dhat, avg_mass, grad_energy, grad_barrier)
 
             prev_distance = collisions.compute_minimum_distance(collision_mesh, vertices)
 
@@ -348,10 +347,10 @@ To add a collision offset, we need to set the ``dmin`` variable. For example, we
             dhat = 1e-4
             dmin = 1e-3
 
-            collisions = ipctk.Collisions()
+            collisions = ipctk.NormalCollisions()
             collisions.build(collision_mesh, vertices, dhat, dmin)
 
-This will then set the ``dmin`` field in all of the ``Collision`` objects stored in the ``collisions``.
+This will then set the ``dmin`` field in all of the ``NormalCollision`` objects stored in the ``collisions``.
 
 .. note::
     Currently, only a single thickness value is supported for the entire mesh.
@@ -371,7 +370,7 @@ Computing the friction dissipative potential is similar to the barrier potential
 
             ipc::TangentialCollisions tangential_collisions;
             tangential_collisions.build(
-                collision_mesh, vertices, collisions, B, barrier_stiffness, mu);
+                collision_mesh, vertices, collisions, B, mu);
 
     .. md-tab-item:: Python
 
@@ -379,9 +378,12 @@ Computing the friction dissipative potential is similar to the barrier potential
 
             tangential_collisions = ipctk.TangentialCollisions()
             tangential_collisions.build(
-                collision_mesh, vertices, collisions, B, barrier_stiffness, mu)
+                collision_mesh, vertices, collisions, B, mu)
+
+Here ``mu`` (:math:`\mu`) is the (global) coefficient of friction.
 
-Here ``mu`` (:math:`\mu`) is the (global) coefficient of friction, and ``barrier_stiffness`` (:math:`\kappa`) is the barrier stiffness.
+.. note::
+   The barrier stiffness (:math:`\kappa`) is retrieved from the normal potential ``B`` and used to compute the lagged normal force magnitudes.
 
 Friction Dissipative Potential
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -486,22 +488,24 @@ The following example determines the maximum step size allowable between the res
 
             Eigen::MatrixXd collision_free_vertices =
                 (vertices_t1 - vertices_t0) * max_step_size + vertices_t0;
-            assert(ipc::is_step_collision_free(mesh, vertices_t0, collision_free_vertices));
+            assert(ipc::is_step_collision_free(
+                collision_mesh, vertices_t0, collision_free_vertices));
 
     .. md-tab-item:: Python
 
         .. code-block:: python
 
-            vertices_t0 = collision_mesh.rest_positions() # vertices at t=0
-            vertices_t1 = vertices_t0.copy()              # vertices at t=1
+            vertices_t0 = collision_mesh.rest_positions # vertices at t=0
+            vertices_t1 = vertices_t0.copy()            # vertices at t=1
             vertices_t1[:, 1] *= 0.01 # squash the mesh in the y-direction
 
             max_step_size = ipctk.compute_collision_free_stepsize(
                     collision_mesh, vertices_t0, vertices_t1)
 
-            collision_free_vertices =
-                (vertices_t1 - vertices_t0) * max_step_size + vertices_t0
-            assert(ipctk.is_step_collision_free(mesh, vertices_t0, collision_free_vertices))
+            collision_free_vertices = (
+                (vertices_t1 - vertices_t0) * max_step_size + vertices_t0)
+            assert ipctk.is_step_collision_free(
+                collision_mesh, vertices_t0, collision_free_vertices)
 
 CCD is comprised of two parts (phases): broad-phase and narrow-phase.
 
@@ -523,7 +527,8 @@ The ``Candidates`` class represents the culled set of candidate pairs and is bui
             ipc::Candidates candidates;
             ipc::LBVH broad_phase;
             candidates.build(
-                mesh, vertices_t0, vertices_t1, /*inflation_radius=*/0.0, broad_phase);
+                collision_mesh, vertices_t0, vertices_t1,
+                /*inflation_radius=*/0.0, &broad_phase);
 
     .. md-tab-item:: Python
 
@@ -531,14 +536,27 @@ The ``Candidates`` class represents the culled set of candidate pairs and is bui
 
             candidates = ipctk.Candidates()
             candidates.build(
-                mesh, vertices_t0, vertices_t1, broad_phase=ipctk.LBVH())
+                collision_mesh, vertices_t0, vertices_t1,
+                broad_phase=ipctk.LBVH())
 
 Possible values for ``broad_phase`` are: ``BruteForce`` (parallel brute force culling), ``HashGrid``, ``SpatialHash`` (implementation from the original IPC codebase), ``LBVH`` (CPU implementation of :cite:t:`Karras2012HPG` using TBB), ``SweepAndPrune`` (a.k.a. Sort-and-Sweep from :cite:t:`Baraff1992PhD`), or ``SweepAndTiniestQueue`` (method of :cite:t:`Belgrod2023Time`; requires CUDA). The default is ``LBVH``.
 
 Narrow-Phase
 ^^^^^^^^^^^^
 
-The narrow phase computes the time of impact between two primitives (e.g., a point and a triangle or two edges in 3D). To do this we utilize the Tight Inclusion CCD method of :cite:t:`Wang2021TightInclusion` for the narrow phase as it is provably conservative (i.e., never misses collisions), accurate (i.e., rarely reports false positives), and efficient.
+The narrow phase computes the time of impact between two primitives (e.g., a point and a triangle or two edges in 3D). Narrow-phase algorithms are implemented as subclasses of ``NarrowPhaseCCD``. The default is ``TightInclusionCCD``, the Tight Inclusion CCD method of :cite:t:`Wang2021TightInclusion`, as it is provably conservative (i.e., never misses collisions), accurate (i.e., rarely reports false positives), and efficient. It is what every function taking a ``narrow_phase_ccd`` argument uses unless you pass something else.
+
+Two other implementations are also available:
+
+- ``AdditiveCCD``, the method of :cite:t:`Li2021CIPC`, is much faster than Tight Inclusion (>100×) and reliable in practice. It does not account for rounding error in its distance computations, so in theory it can miss collisions, but its default margin (10% of the initial separation) is large enough to avoid this. The cost of that margin is a less accurate time of impact and more false positives. Shrinking it -- pushing ``conservative_rescaling`` toward ``1.0`` -- tightens the time of impact but is what can introduce false negatives. Tight Inclusion accounts for the rounding error, so it can reduce its tolerance without that trade-off :cite:p:`Belgrod2023Time`.
+- ``InexactCCD``, the original method from the IPC codebase, is disabled by default. To use it, set the ``IPC_TOOLKIT_WITH_INEXACT_CCD`` CMake option to ``ON``.
+
+All three use the same expression for their conservative margin, stopping short of contact at a separation of
+
+.. math::
+   d_\min + (1 - r)(d_0 - d_\min),
+
+where :math:`r` is ``conservative_rescaling`` and :math:`d_0` is the distance at :math:`t=0`. Note that the margin is a fraction of the initial separation *in excess of* :math:`d_\min`, not of the raw distance. ``TightInclusionCCD`` differs only in additionally capping the second term at :math:`10^{-4}` (see the note below), which is why it usually reports a time of impact much closer to the exact one than ``AdditiveCCD`` does for the same query.
 
 The following example shows how to use the narrow phase to determine if a point is colliding with a triangle (static in this case).
 
@@ -548,7 +566,7 @@ The following example shows how to use the narrow phase to determine if a point
 
         .. code-block:: c++
 
-            #include <ipc/ccd/ccd.hpp>
+            #include <ipc/ccd/tight_inclusion_ccd.hpp>
 
             // ...
 
@@ -565,10 +583,10 @@ The following example shows how to use the narrow phase to determine if a point
             Eigen::Vector3d t2_t1 = t2_t0; // triangle vertex 2 at t=1
 
             double toi; // output time of impact
-            bool is_colliding = ipc::point_triangle_ccd(
+            bool is_colliding = ipc::TightInclusionCCD().point_triangle_ccd(
                 p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1, toi);
             assert(is_colliding);
-            assert(abs(toi - 0.5) < 1e-8);
+            assert(toi <= 0.5); // conservative estimate of the exact TOI of 0.5
 
     .. md-tab-item:: Python
 
@@ -591,12 +609,35 @@ The following example shows how to use the narrow phase to determine if a point
 
             # returns a boolean indicating if the point is colliding with the triangle
             # and the time of impact (TOI)
-            is_colliding, toi = ipctk.point_triangle_ccd(
+            is_colliding, toi = ipctk.TightInclusionCCD().point_triangle_ccd(
                 p_t0, t0_t0, t1_t0, t2_t0, p_t1, t0_t1, t1_t1, t2_t1)
-            assert(is_colliding)
-            assert(abs(toi - 0.5) < 1e-8)
+            assert is_colliding
+            assert toi <= 0.5 # conservative estimate of the exact TOI of 0.5
 
-Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function that can be used to determine if the face-vertex pairing is colliding:
+.. note::
+   The returned time of impact is a *conservative* under-estimate, so do not test
+   it for exact equality. For the query above the exact TOI is ``0.5``, but the
+   returned value is slightly less than that.
+
+   ``TightInclusionCCD`` achieves this primarily by inflating the *distance* at
+   which the query stops rather than by scaling the resulting TOI. It runs the
+   narrow-phase query with a minimum separation of
+
+   .. math::
+      d_\min + \min\left((1 - r)(d_0 - d_\min),\ 10^{-4}\right),
+
+   where :math:`r` is ``conservative_rescaling`` and :math:`d_0` is the distance
+   at :math:`t=0`. The TOI therefore comes back early by roughly that separation
+   divided by the distance travelled. Note the :math:`10^{-4}` cap: for
+   well-separated primitives it is what binds, so changing
+   ``conservative_rescaling`` has no effect on the result.
+
+   The TOI itself is multiplied by ``conservative_rescaling`` only in a fallback
+   path, when the query above returns a TOI below
+   ``TightInclusionCCD::SMALL_TOI``; the query is then rerun with the true
+   :math:`d_\min` and the result scaled to keep it away from zero.
+
+Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function that can be used to determine if the face-vertex pairing is colliding. It takes the *stencil* vertices (the four vertices of the face-vertex pair), which you can gather from the full vertex matrix using ``CollisionStencil::dof``:
 
 .. md-tab-set::
 
@@ -606,9 +647,13 @@ Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function tha
 
             ipc::FaceVertexCandidate candidate = ...; // face-vertex candidate
 
+            const Eigen::MatrixXi& edges = collision_mesh.edges();
+            const Eigen::MatrixXi& faces = collision_mesh.faces();
+
             double toi; // output time of impact
             bool is_colliding = candidate.ccd(
-                vertices_t0, vertices_t1, collision_mesh.edges(), collision_mesh.faces(), toi);
+                candidate.dof(vertices_t0, edges, faces),
+                candidate.dof(vertices_t1, edges, faces), toi);
 
     .. md-tab-item:: Python
 
@@ -616,12 +661,15 @@ Alternatively, the ``FaceVertexCandidate`` class contains a ``ccd`` function tha
 
             candidate = ... # face-vertex candidate
 
+            edges, faces = collision_mesh.edges, collision_mesh.faces
+
             # returns a boolean indicating if the point is colliding with the triangle
             # and the time of impact (TOI)
             is_colliding, toi = candidate.ccd(
-                vertices_t0, vertices_t1, collision_mesh.edges, collision_mesh.faces)
+                candidate.dof(vertices_t0, edges, faces),
+                candidate.dof(vertices_t1, edges, faces))
 
-The same can be done for point-edge collisions using the ``point_edge_ccd`` function or ``EdgeVertexCandidate`` class and for edge-edge collisions using the ``edge_edge_ccd`` function or ``EdgeEdgeCandidate`` class.
+The same can be done for point-edge collisions using the ``NarrowPhaseCCD::point_edge_ccd`` method or ``EdgeVertexCandidate`` class and for edge-edge collisions using the ``NarrowPhaseCCD::edge_edge_ccd`` method or ``EdgeEdgeCandidate`` class.
 
 .. _minimum-separation-ccd:
 
@@ -644,7 +692,8 @@ To do this, we need to set the ``min_distance`` parameter when calling ``is_step
             Eigen::MatrixXd collision_free_vertices =
                 (vertices_t1 - vertices_t0) * max_step_size + vertices_t0;
             assert(ipc::is_step_collision_free(
-                mesh, vertices_t0, collision_free_vertices, /*min_distance=*/1e-4));
+                collision_mesh, vertices_t0, collision_free_vertices,
+                /*min_distance=*/1e-4));
 
     .. md-tab-item:: Python
 
@@ -653,7 +702,8 @@ To do this, we need to set the ``min_distance`` parameter when calling ``is_step
             max_step_size = ipctk.compute_collision_free_stepsize(
                     collision_mesh, vertices_t0, vertices_t1, min_distance=1e-4)
 
-            collision_free_vertices =
-                (vertices_t1 - vertices_t0) * max_step_size + vertices_t0
-            assert(ipctk.is_step_collision_free(
-                mesh, vertices_t0, collision_free_vertices, min_distance=1e-4))
+            collision_free_vertices = (
+                (vertices_t1 - vertices_t0) * max_step_size + vertices_t0)
+            assert ipctk.is_step_collision_free(
+                collision_mesh, vertices_t0, collision_free_vertices,
+                min_distance=1e-4)
diff --git a/docs/source/tutorials/nonlinear_ccd.rst b/docs/source/tutorials/nonlinear_ccd.rst
index 347b4a7e6..2c9b9eb4a 100644
--- a/docs/source/tutorials/nonlinear_ccd.rst
+++ b/docs/source/tutorials/nonlinear_ccd.rst
@@ -5,27 +5,27 @@ We also implement CCD of nonlinear trajectories (of linear geometry) using the m
 
 The method works by transforming the nonlinear trajectories into (adaptive) piecewise linear trajectories with an envelope/minimum separation around each piece, enclosing the nonlinear trajectory. The method then performs CCD on the piecewise linear trajectories to find the earliest time of impact.
 
-We provide the following functions to perform nonlinear CCD:
+Nonlinear CCD is provided by the ``NonlinearCCD`` class, which exposes the following methods:
 
 .. md-tab-set::
 
     .. md-tab-item:: C++
 
-        * :cpp:func:`ipc::point_point_nonlinear_ccd`,
-        * :cpp:func:`ipc::point_edge_nonlinear_ccd`,
-        * :cpp:func:`ipc::edge_edge_nonlinear_ccd`, and
-        * :cpp:func:`ipc::point_triangle_nonlinear_ccd`.
+        * :cpp:func:`ipc::NonlinearCCD::point_point_ccd`,
+        * :cpp:func:`ipc::NonlinearCCD::point_edge_ccd`,
+        * :cpp:func:`ipc::NonlinearCCD::edge_edge_ccd`, and
+        * :cpp:func:`ipc::NonlinearCCD::point_triangle_ccd`.
 
-        Each of these functions take as input a :cpp:class:`ipc::NonlinearTrajectory` object for the endpoints of the linear geometry.
+        Each of these methods take as input a :cpp:class:`ipc::NonlinearTrajectory` object for the endpoints of the linear geometry.
 
     .. md-tab-item:: Python
 
-        * :py:func:`ipctk.point_point_nonlinear_ccd`,
-        * :py:func:`ipctk.point_edge_nonlinear_ccd`,
-        * :py:func:`ipctk.edge_edge_nonlinear_ccd`, and
-        * :py:func:`ipctk.point_triangle_nonlinear_ccd`.
+        * :py:meth:`ipctk.NonlinearCCD.point_point_ccd`,
+        * :py:meth:`ipctk.NonlinearCCD.point_edge_ccd`,
+        * :py:meth:`ipctk.NonlinearCCD.edge_edge_ccd`, and
+        * :py:meth:`ipctk.NonlinearCCD.point_triangle_ccd`.
 
-        Each of these functions take as input a :py:class:`ipctk.NonlinearTrajectory` object for the endpoints of the linear geometry.
+        Each of these methods take as input a :py:class:`ipctk.NonlinearTrajectory` object for the endpoints of the linear geometry.
 
 For example, the following code defines a rigid trajectory in 2D in order to perform nonlinear CCD between a point and edge:
 
@@ -207,7 +207,7 @@ The following code snippet shows an example of how to use interval arithmetic to
                 Eigen::ConstRef<VectorMax3d> center,
                 Eigen::ConstRef<VectorMax3d> point,
                 const double omega,
-                const Interval& t)
+                const filib::Interval& t)
             {
                 // 2×2 matrix of intervals representing the rotation matrix
                 Matrix2I R;
diff --git a/docs/source/tutorials/ogc.rst b/docs/source/tutorials/ogc.rst
index 995cdfa14..04cd3e1cb 100644
--- a/docs/source/tutorials/ogc.rst
+++ b/docs/source/tutorials/ogc.rst
@@ -263,18 +263,18 @@ Putting it all together, a single simulation step using OGC looks like this:
             // 2. Warm start (Predict & Initialize)
             // x is current position, pred_x is x^t + dt * v^t
             trust_region.warm_start_time_step(
-                mesh, x, pred_x, collisions, dhat);
+                collision_mesh, x, pred_x, collisions, dhat);
 
             // 3. Solver Loop
             for (int i = 0; i < max_iterations; ++i) {
                 // Update trust region if too many vertices hit the bound in previous step
-                trust_region.update_if_needed(mesh, x, collisions, dhat);
+                trust_region.update_if_needed(collision_mesh, x, collisions, dhat);
 
                 // Compute search direction (Solver specific)
                 Eigen::MatrixXd dx = compute_search_direction(x, ...);
 
                 // Filter the step to respect OGC bounds
-                trust_region.filter_step(mesh, x, dx);
+                trust_region.filter_step(collision_mesh, x, dx);
 
                 // Update positions
                 x += dx;
@@ -293,18 +293,18 @@ Putting it all together, a single simulation step using OGC looks like this:
 
             # 2. Warm start (Predict & Initialize)
             trust_region.warm_start_time_step(
-                mesh, x, pred_x, collisions, dhat)
+                collision_mesh, x, pred_x, collisions, dhat)
 
             # 3. Solver Loop
             for i in range(max_iterations):
                 # Update trust region if needed
-                trust_region.update_if_needed(mesh, x, collisions, dhat)
+                trust_region.update_if_needed(collision_mesh, x, collisions, dhat)
 
                 # Compute search direction (Solver specific)
                 dx = compute_search_direction(x, ...)
 
                 # Filter the step to respect OGC bounds
-                trust_region.filter_step(mesh, x, dx)
+                trust_region.filter_step(collision_mesh, x, dx)
 
                 # Update positions
                 x += dx
@@ -416,9 +416,9 @@ Using ``planar_filter_step``
         .. code-block:: c++
 
             // Inside the solver loop, replace:
-            //   trust_region.filter_step(mesh, x, dx);
+            //   trust_region.filter_step(collision_mesh, x, dx);
             // with:
-            trust_region.planar_filter_step(mesh, x, dx);
+            trust_region.planar_filter_step(collision_mesh, x, dx);
 
             // Optionally tune the relaxation ratio via the struct member:
             // trust_region.relaxed_radius_scaling = 0.9; // default
@@ -428,9 +428,9 @@ Using ``planar_filter_step``
         .. code-block:: python
 
             # Inside the solver loop, replace:
-            #   trust_region.filter_step(mesh, x, dx)
+            #   trust_region.filter_step(collision_mesh, x, dx)
             # with:
-            trust_region.planar_filter_step(mesh, x, dx)
+            trust_region.planar_filter_step(collision_mesh, x, dx)
 
             # Optionally tune the relaxation ratio via the struct member:
             # trust_region.relaxed_radius_scaling = 0.9  # default
@@ -451,18 +451,18 @@ Full Optimization Loop with Planar-DAT
 
             // 2. Warm start (Predict & Initialize)
             trust_region.warm_start_time_step(
-                mesh, x, pred_x, collisions, dhat);
+                collision_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);
+                trust_region.update_if_needed(collision_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);
+                trust_region.planar_filter_step(collision_mesh, x, dx);
 
                 // Update positions
                 x += dx;
@@ -481,18 +481,18 @@ Full Optimization Loop with Planar-DAT
 
             # 2. Warm start (Predict & Initialize)
             trust_region.warm_start_time_step(
-                mesh, x, pred_x, collisions, dhat)
+                collision_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)
+                trust_region.update_if_needed(collision_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)
+                trust_region.planar_filter_step(collision_mesh, x, dx)
 
                 # Update positions
                 x += dx
diff --git a/docs/source/tutorials/simulation.rst b/docs/source/tutorials/simulation.rst
index 7faf069be..ab2eb7954 100644
--- a/docs/source/tutorials/simulation.rst
+++ b/docs/source/tutorials/simulation.rst
@@ -38,7 +38,7 @@ From the full (volumetric) mesh vertices and surface edges/faces which index int
             // TODO: Show how to load a volumetric mesh from a file (e.g., using MshIO)
 
             // Faces of the surface mesh with indices into full_rest_positions
-            Eigen::MatrixXd faces;
+            Eigen::MatrixXi faces;
             igl::boundary_facets(tets, faces);
 
             // Edges of the surface mesh with indices into full_rest_positions
@@ -71,7 +71,7 @@ This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. How
         .. code-block:: c++
 
             // Convert full vertices to surface vertices
-            Eigen::VectorXd vertices = collision_mesh.vertices(full_vertices);
+            Eigen::MatrixXd vertices = collision_mesh.vertices(full_vertices);
 
             // Construct the set of collisions
             ipc::NormalCollisions collisions;
@@ -84,7 +84,7 @@ This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. How
             double b = B(collisions, collision_mesh, vertices);
 
             // Convert full velocities to surface velocities
-            Eigen::VectorXd velocities = collision_mesh.map_displacements(full_velocities);
+            Eigen::MatrixXd velocities = collision_mesh.map_displacements(full_velocities);
 
             // Construct the set of friction collisions
             ipc::TangentialCollisions tangential_collisions;
@@ -103,7 +103,7 @@ This ``CollisionMesh`` can then be used just as any other ``CollisionMesh``. How
             vertices = collision_mesh.vertices(full_vertices)
 
             # Construct the set of collisions
-            collisions = ipctk.Collisions()
+            collisions = ipctk.NormalCollisions()
             collisions.build(collision_mesh, vertices, dhat)
 
             # Construct a barrier potential
@@ -144,12 +144,12 @@ When computing the gradient and Hessian of the potentials, the derivatives will
 
         .. code-block:: python
 
-            B = BarrierPotential(dhat, stiffness)
+            B = ipctk.BarrierPotential(dhat, stiffness)
 
-            grad = B.gradient(collision, collision_mesh, vertices)
+            grad = B.gradient(collisions, collision_mesh, vertices)
             grad_full = collision_mesh.to_full_dof(grad)
 
-            hess = B.hessian(collision, collision_mesh, vertices)
+            hess = B.hessian(collisions, collision_mesh, vertices)
             hess_full = collision_mesh.to_full_dof(hess)
 
 Codimensional Vertices
@@ -170,9 +170,14 @@ In some cases, the collision mesh vertices are not the same as the surface verti
             std::vector<bool> is_on_surface = ipc::CollisionMesh::construct_is_on_surface(
                 full_rest_positions.rows(), boundary_edges, codim_vertices);
 
-            // Construct the collision mesh from the is_on_surface vector and full mesh data
+            // is_orient_vertex marks the vertices with a well-defined orientation
+            // (i.e., those usable for signed distances). Codimensional vertices
+            // have no orientation, so false is a safe default for all vertices.
+            std::vector<bool> is_orient_vertex(full_rest_positions.rows(), false);
+
+            // Construct the collision mesh from the masks and full mesh data
             ipc::CollisionMesh collision_mesh(
-                is_on_surface, full_rest_positions, edges, faces);
+                is_on_surface, is_orient_vertex, full_rest_positions, edges, faces);
 
     .. md-tab-item:: Python
 
@@ -185,9 +190,14 @@ In some cases, the collision mesh vertices are not the same as the surface verti
             is_on_surface = ipctk.CollisionMesh.construct_is_on_surface(
                 len(full_rest_positions), boundary_edges, codim_vertices)
 
-            # Construct the collision mesh from the is_on_surface vector and full mesh data
+            # is_orient_vertex marks the vertices with a well-defined orientation
+            # (i.e., those usable for signed distances). Codimensional vertices
+            # have no orientation, so False is a safe default for all vertices.
+            is_orient_vertex = [False] * len(full_rest_positions)
+
+            # Construct the collision mesh from the masks and full mesh data
             collision_mesh = ipctk.CollisionMesh(
-                is_on_surface, full_rest_positions, edges, faces)
+                is_on_surface, is_orient_vertex, full_rest_positions, edges, faces)
 
 Nonlinear Bases and Curved Meshes
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -209,8 +219,8 @@ While IPC cannot directly handle nonlinear finite element bases and/or curved me
             Eigen::MatrixXd proxy_rest_positions;
             Eigen::MatrixXi proxy_edges, proxy_faces;
             // Load the proxy mesh from a file
-            igl::read_triangle_mesh("proxy.ply", rest_positions, faces);
-            igl::edges(faces, edges);
+            igl::read_triangle_mesh("proxy.ply", proxy_rest_positions, proxy_faces);
+            igl::edges(proxy_faces, proxy_edges);
             // Or build it from the volumetric mesh
 
             // Linear map from the finite element mesh to the collision proxy
@@ -225,8 +235,8 @@ While IPC cannot directly handle nonlinear finite element bases and/or curved me
 
             # Finite element mesh
             fe_mesh = meshio.read("mesh.msh")
-            fe_rest_positions = mesh.points
-            tets = mesh.cells_dict["tetra"]
+            fe_rest_positions = fe_mesh.points
+            tets = fe_mesh.cells_dict["tetra"]
 
             # Collision proxy mesh
             # Load the proxy mesh from a file
@@ -239,10 +249,10 @@ While IPC cannot directly handle nonlinear finite element bases and/or curved me
             # Linear map from the finite element mesh to the collision proxy
             displacement_map = ... # build or load the displacement map
 
-            collision_mesh = CollisionMesh(
+            collision_mesh = ipctk.CollisionMesh(
                 proxy_rest_positions, proxy_edges, proxy_faces, displacement_map)
 
-We can then map the displacements using ``collision_mesh.map_displacement(fe_displacements)`` or directly get the displaced proxy mesh vertices using ``collision_mesh.displace_vertices(fe_displacements)``. Similarly, we can map forces/potential gradients using ``collision_mesh.to_full_dof(collision_forces)`` or force Jacobians/potential Hessians using ``collision_mesh.to_full_dof(potential_hessian)``.
+We can then map the displacements using ``collision_mesh.map_displacements(fe_displacements)`` or directly get the displaced proxy mesh vertices using ``collision_mesh.displace_vertices(fe_displacements)``. Similarly, we can map forces/potential gradients using ``collision_mesh.to_full_dof(collision_forces)`` or force Jacobians/potential Hessians using ``collision_mesh.to_full_dof(potential_hessian)``.
 
 .. warning::
     The function ``CollisionMesh::vertices(full_positions)`` should not be used in this case because the rest positions used to construct the ``CollisionMesh`` are not the same as the finite element mesh's rest positions. Instead, use ``CollisionMesh::displace_vertices(fe_displacements)`` where ``fe_displacements`` is already the solution of the PDE or can be computed as ``fe_displacements = fe_positions - fe_rest_positions`` from deformed and rest positions.
@@ -257,10 +267,24 @@ To remedy this, we can project the Hessian onto the positive semidefinite (PSD)
 
     .. md-tab-item:: C++
 
-        - ``ProjectToPSD::CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`.
-        - ``ProjectToPSD::ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
+        - ``PSDProjectionMethod::NONE``: Do not project the Hessian. This is the default.
+        - ``PSDProjectionMethod::CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`.
+        - ``PSDProjectionMethod::ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
+
+        .. code-block:: c++
+
+            Eigen::SparseMatrix<double> hess = B.hessian(
+                collisions, collision_mesh, vertices,
+                ipc::PSDProjectionMethod::CLAMP);
 
     .. md-tab-item:: Python
 
-        - ``ProjectToPSD.CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`.
-        - ``ProjectToPSD.ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
\ No newline at end of file
+        - ``PSDProjectionMethod.NONE``: Do not project the Hessian. This is the default.
+        - ``PSDProjectionMethod.CLAMP``: Clamp the negative eigenvalues of the Hessian to 0. This is the same as used by :cite:t:`Li2020IPC`.
+        - ``PSDProjectionMethod.ABS``: Set the negative eigenvalues of the Hessian to their absolute value. This is the method proposed by :cite:t:`Chen2024Stabler`.
+
+        .. code-block:: python
+
+            hess = B.hessian(
+                collisions, collision_mesh, vertices,
+                project_hessian_to_psd=ipctk.PSDProjectionMethod.CLAMP)
\ No newline at end of file
diff --git a/python/src/collisions/normal/normal_collisions.cpp b/python/src/collisions/normal/normal_collisions.cpp
index 901959ab1..d394cbc98 100644
--- a/python/src/collisions/normal/normal_collisions.cpp
+++ b/python/src/collisions/normal/normal_collisions.cpp
@@ -19,6 +19,21 @@ void define_smooth_collisions(py::module_& m, const std::string& name)
 {
     py::class_<SmoothCollisions>(m, name.c_str())
         .def(py::init())
+        .def(
+            "compute_adaptive_dhat", &SmoothCollisions::compute_adaptive_dhat,
+            R"ipc_Qu8mg5v7(
+            Compute the per-element adaptive dhat from the rest configuration.
+
+            Note:
+                Must be called before build() when using use_adaptive_dhat=True.
+
+            Parameters:
+                mesh: The collision mesh.
+                vertices: Vertices of the collision mesh.
+                params: SmoothContactParameters.
+                broad_phase: Broad phase method.
+            )ipc_Qu8mg5v7",
+            "mesh"_a, "vertices"_a, "params"_a, "broad_phase"_a = nullptr)
         .def(
             "build",
             py::overload_cast<
diff --git a/python/src/common.hpp b/python/src/common.hpp
index 2fc728f4c..719cd8816 100644
--- a/python/src/common.hpp
+++ b/python/src/common.hpp
@@ -15,6 +15,34 @@ using namespace py::literals;
 #include <Eigen/Core>
 #include <Eigen/Sparse>
 
+/// @brief Check that a parameter is strictly positive.
+/// @throws py::value_error (Python ValueError) if value is not positive.
+/// @note The C++ API only enforces this with an assert(), which is compiled
+///       out under NDEBUG. Validating here means Python users get an
+///       exception rather than undefined behavior in a release build.
+/// @note Written as !(value > 0) so that NaN is rejected too.
+inline void assert_positive(const double value, const std::string& name)
+{
+    if (!(value > 0)) {
+        throw py::value_error(
+            "Parameter " + name + " has invalid value: expected " + name
+            + " > 0 but got " + name + " = " + std::to_string(value));
+    }
+}
+
+/// @brief Check that a shared_ptr parameter is not None.
+/// @throws py::value_error (Python ValueError) if ptr is null.
+template <typename T>
+inline void
+assert_not_none(const std::shared_ptr<T>& ptr, const std::string& name)
+{
+    if (ptr == nullptr) {
+        throw py::value_error(
+            "Parameter " + name + " has invalid value: expected " + name
+            + " to be a valid object but got None");
+    }
+}
+
 template <typename Derived>
 void assert_2D_or_3D_vector(
     const Eigen::MatrixBase<Derived>& v, const std::string& name)
diff --git a/python/src/potentials/barrier_potential.cpp b/python/src/potentials/barrier_potential.cpp
index 1d8f74561..25224fc9c 100644
--- a/python/src/potentials/barrier_potential.cpp
+++ b/python/src/potentials/barrier_potential.cpp
@@ -9,7 +9,12 @@ void define_barrier_potential(py::module_& m)
 {
     py::class_<BarrierPotential, NormalPotential>(m, "BarrierPotential")
         .def(
-            py::init<const double, const double, const bool>(),
+            py::init([](const double dhat, const double stiffness,
+                        const bool use_physical_barrier) {
+                assert_positive(dhat, "dhat");
+                assert_positive(stiffness, "stiffness");
+                return BarrierPotential(dhat, stiffness, use_physical_barrier);
+            }),
             R"ipc_Qu8mg5v7(
             Construct a barrier potential.
 
@@ -17,12 +22,21 @@ void define_barrier_potential(py::module_& m)
                 dhat: The activation distance of the barrier.
                 stiffness: The stiffness of the barrier.
                 use_physical_barrier: Whether to use the physical barrier.
+
+            Raises:
+                ValueError: If dhat or stiffness is not positive.
             )ipc_Qu8mg5v7",
             "dhat"_a, "stiffness"_a, "use_physical_barrier"_a = false)
         .def(
-            py::init<
-                const std::shared_ptr<Barrier>, const double, const double,
-                const bool>(),
+            py::init([](std::shared_ptr<Barrier> barrier, const double dhat,
+                        const double stiffness,
+                        const bool use_physical_barrier) {
+                assert_not_none(barrier, "barrier");
+                assert_positive(dhat, "dhat");
+                assert_positive(stiffness, "stiffness");
+                return BarrierPotential(
+                    std::move(barrier), dhat, stiffness, use_physical_barrier);
+            }),
             R"ipc_Qu8mg5v7(
             Construct a barrier potential.
 
@@ -31,18 +45,40 @@ void define_barrier_potential(py::module_& m)
                 dhat: The activation distance of the barrier.
                 stiffness: The stiffness of the barrier.
                 use_physical_barrier: Whether to use the physical barrier.
+
+            Raises:
+                ValueError: If barrier is None, or dhat or stiffness is not positive.
             )ipc_Qu8mg5v7",
             "barrier"_a, "dhat"_a, "stiffness"_a,
             "use_physical_barrier"_a = false)
         .def_property(
-            "dhat", &BarrierPotential::dhat, &BarrierPotential::set_dhat,
-            "Barrier activation distance.")
+            "dhat", &BarrierPotential::dhat,
+            [](BarrierPotential& self, const double dhat) {
+                assert_positive(dhat, "dhat");
+                self.set_dhat(dhat);
+            },
+            "Barrier activation distance. Must be positive.")
+        .def_property(
+            "stiffness", &BarrierPotential::stiffness,
+            [](BarrierPotential& self, const double stiffness) {
+                assert_positive(stiffness, "stiffness");
+                self.set_stiffness(stiffness);
+            },
+            "Barrier stiffness. Must be positive.")
+        .def_property(
+            "use_physical_barrier", &BarrierPotential::use_physical_barrier,
+            &BarrierPotential::set_use_physical_barrier,
+            "Whether to use the physical barrier.")
         .def_property(
             "barrier",
             py::cpp_function(
                 &BarrierPotential::barrier, py::return_value_policy::reference),
-            &BarrierPotential::set_barrier,
-            "Barrier function used to compute the potential.");
+            [](BarrierPotential& self,
+               const std::shared_ptr<Barrier>& barrier) {
+                assert_not_none(barrier, "barrier");
+                self.set_barrier(barrier);
+            },
+            "Barrier function used to compute the potential. Must not be None.");
 }
 
 void define_smooth_potential(py::module_& m)
@@ -73,9 +109,15 @@ void define_smooth_potential(py::module_& m)
         .def_readonly("beta_t", &SmoothContactParameters::beta_t)
         .def_readonly("alpha_n", &SmoothContactParameters::alpha_n)
         .def_readonly("beta_n", &SmoothContactParameters::beta_n)
-        .def_readonly("r", &SmoothContactParameters::r);
+        .def_readonly("r", &SmoothContactParameters::r)
+        .def_property(
+            "adaptive_dhat_ratio",
+            &SmoothContactParameters::adaptive_dhat_ratio,
+            &SmoothContactParameters::set_adaptive_dhat_ratio,
+            "Ratio of the distance to the interaction set in the rest "
+            "configuration used as the per-element adaptive dhat.");
 
-    py::class_<SmoothContactPotential>(m, "SmoothPotential")
+    py::class_<SmoothContactPotential>(m, "SmoothContactPotential")
         .def(
             py::init<const SmoothContactParameters&>(),
             R"ipc_Qu8mg5v7(
diff --git a/python/tests/test_potentials.py b/python/tests/test_potentials.py
new file mode 100644
index 000000000..4c13e4900
--- /dev/null
+++ b/python/tests/test_potentials.py
@@ -0,0 +1,264 @@
+"""Tests for the potential bindings.
+
+Covers the Python-side input validation (which the C++ API only enforces with
+assert(), compiled out under NDEBUG) and the smooth-contact/GCP bindings.
+
+Uses unittest.TestCase so that assertRaises is available under both nose2 (the
+runner used in CI) and pytest, without adding a pytest dependency.
+"""
+
+import unittest
+
+import numpy as np
+from find_ipctk import ipctk
+from utils import load_mesh
+
+# two-cubes-close.ply has a ~0.069 gap between the cubes, so dhat=0.1 activates
+# a few hundred collisions while staying fast (~1 ms to build).
+DHAT = 0.1
+STIFFNESS = 1.0
+
+
+def two_cubes():
+    V, E, F = load_mesh("two-cubes-close.ply")
+    mesh = ipctk.CollisionMesh(V, E, F)
+    return mesh, mesh.rest_positions
+
+
+def normal_collisions(mesh, vertices, dhat=DHAT):
+    collisions = ipctk.NormalCollisions()
+    collisions.build(mesh, vertices, dhat)
+    return collisions
+
+
+class TestBarrierPotentialValidation(unittest.TestCase):
+    """The C++ ctor/setters assert dhat > 0, stiffness > 0, barrier != nullptr.
+
+    assert() is compiled out under NDEBUG, so the bindings must validate and
+    raise ValueError instead of admitting undefined behavior in a release build.
+    """
+
+    INVALID = [0.0, -1.0, -1e-3, float("nan")]
+
+    def test_ctor_rejects_invalid_dhat(self):
+        for dhat in self.INVALID:
+            with self.subTest(dhat=dhat):
+                with self.assertRaises(ValueError):
+                    ipctk.BarrierPotential(dhat, STIFFNESS)
+
+    def test_ctor_rejects_invalid_stiffness(self):
+        for stiffness in self.INVALID:
+            with self.subTest(stiffness=stiffness):
+                with self.assertRaises(ValueError):
+                    ipctk.BarrierPotential(DHAT, stiffness)
+
+    def test_ctor_rejects_none_barrier(self):
+        with self.assertRaises(ValueError):
+            ipctk.BarrierPotential(None, DHAT, STIFFNESS)
+
+    def test_barrier_ctor_rejects_invalid_dhat(self):
+        for dhat in self.INVALID:
+            with self.subTest(dhat=dhat):
+                with self.assertRaises(ValueError):
+                    ipctk.BarrierPotential(
+                        ipctk.ClampedLogBarrier(), dhat, STIFFNESS)
+
+    def test_dhat_setter_rejects_invalid(self):
+        for dhat in self.INVALID:
+            with self.subTest(dhat=dhat):
+                B = ipctk.BarrierPotential(DHAT, STIFFNESS)
+                with self.assertRaises(ValueError):
+                    B.dhat = dhat
+
+    def test_stiffness_setter_rejects_invalid(self):
+        for stiffness in self.INVALID:
+            with self.subTest(stiffness=stiffness):
+                B = ipctk.BarrierPotential(DHAT, STIFFNESS)
+                with self.assertRaises(ValueError):
+                    B.stiffness = stiffness
+
+    def test_barrier_setter_rejects_none(self):
+        B = ipctk.BarrierPotential(DHAT, STIFFNESS)
+        with self.assertRaises(ValueError):
+            B.barrier = None
+
+    def test_state_unchanged_after_rejected_set(self):
+        """A rejected assignment must not partially apply."""
+        B = ipctk.BarrierPotential(DHAT, STIFFNESS)
+        for attr, bad in (("dhat", 0.0), ("stiffness", -1.0)):
+            with self.subTest(attr=attr):
+                with self.assertRaises(ValueError):
+                    setattr(B, attr, bad)
+        self.assertEqual(B.dhat, DHAT)
+        self.assertEqual(B.stiffness, STIFFNESS)
+
+    def test_valid_values_accepted(self):
+        B = ipctk.BarrierPotential(DHAT, STIFFNESS)
+        B.dhat = 2e-3
+        self.assertEqual(B.dhat, 2e-3)
+        B.stiffness = 5.0
+        self.assertEqual(B.stiffness, 5.0)
+        B.barrier = ipctk.ClampedLogBarrier()
+        # use_physical_barrier has no precondition
+        B.use_physical_barrier = True
+        self.assertTrue(B.use_physical_barrier)
+        B.use_physical_barrier = False
+        self.assertFalse(B.use_physical_barrier)
+        # Tiny but positive must be allowed; only <= 0 and NaN are rejected.
+        ipctk.BarrierPotential(1e-300, 1e-300)
+        ipctk.BarrierPotential(ipctk.ClampedLogBarrier(), DHAT, STIFFNESS)
+
+
+class TestBarrierPotentialProperties(unittest.TestCase):
+    """The stiffness/use_physical_barrier properties must reach the evaluation
+    path, not merely round-trip through a stored field."""
+
+    @classmethod
+    def setUpClass(cls):
+        cls.mesh, cls.vertices = two_cubes()
+        cls.collisions = normal_collisions(cls.mesh, cls.vertices)
+        assert len(cls.collisions) > 0, "fixture produced no collisions"
+
+    def test_ctor_roundtrip(self):
+        B = ipctk.BarrierPotential(DHAT, 2.5, use_physical_barrier=True)
+        self.assertEqual(B.dhat, DHAT)
+        self.assertEqual(B.stiffness, 2.5)
+        self.assertTrue(B.use_physical_barrier)
+
+    def test_stiffness_scales_potential_linearly(self):
+        """kappa multiplies the barrier potential, so tripling it must triple
+        the value. Guards against a setter that stores but is never read."""
+        B1 = ipctk.BarrierPotential(DHAT, 1.0)
+        B3 = ipctk.BarrierPotential(DHAT, 1.0)
+        B3.stiffness = 3.0
+
+        p1 = B1(self.collisions, self.mesh, self.vertices)
+        p3 = B3(self.collisions, self.mesh, self.vertices)
+        self.assertGreater(p1, 0.0)
+        # rtol, not exact: the sum is a parallel reduction, so the operand
+        # order (and thus rounding) is not guaranteed to be reproducible.
+        np.testing.assert_allclose(p3, 3.0 * p1, rtol=1e-9)
+
+        g1 = B1.gradient(self.collisions, self.mesh, self.vertices)
+        g3 = B3.gradient(self.collisions, self.mesh, self.vertices)
+        np.testing.assert_allclose(g3, 3.0 * g1, rtol=1e-9)
+
+    def test_use_physical_barrier_setter_matches_ctor(self):
+        """Setting the property must be equivalent to passing the ctor kwarg."""
+        for flag in (False, True):
+            with self.subTest(use_physical_barrier=flag):
+                via_ctor = ipctk.BarrierPotential(
+                    DHAT, 2.5, use_physical_barrier=flag)
+                via_setter = ipctk.BarrierPotential(DHAT, 1.0)
+                via_setter.stiffness = 2.5
+                via_setter.use_physical_barrier = flag
+
+                args = (self.collisions, self.mesh, self.vertices)
+                np.testing.assert_allclose(
+                    via_setter(*args), via_ctor(*args), rtol=1e-9)
+                np.testing.assert_allclose(
+                    via_setter.gradient(*args), via_ctor.gradient(*args),
+                    rtol=1e-9)
+                np.testing.assert_allclose(
+                    via_setter.hessian(*args).todense(),
+                    via_ctor.hessian(*args).todense(), rtol=1e-9)
+
+    def test_use_physical_barrier_changes_result(self):
+        """The flag must actually alter the potential, otherwise the two
+        branches of the test above would agree trivially."""
+        off = ipctk.BarrierPotential(DHAT, STIFFNESS, False)
+        on = ipctk.BarrierPotential(DHAT, STIFFNESS, True)
+        args = (self.collisions, self.mesh, self.vertices)
+        self.assertNotEqual(off(*args), on(*args))
+
+
+class TestSmoothContactParameters(unittest.TestCase):
+    def test_adaptive_dhat_ratio_default(self):
+        params = ipctk.SmoothContactParameters(DHAT, 0.5, 0.0, 0.1, 0.0, 2)
+        self.assertEqual(params.adaptive_dhat_ratio, 0.5)
+
+    def test_adaptive_dhat_ratio_roundtrip(self):
+        params = ipctk.SmoothContactParameters(DHAT, 0.5, 0.0, 0.1, 0.0, 2)
+        for ratio in (0.1, 0.25, 0.9):
+            params.adaptive_dhat_ratio = ratio
+            self.assertEqual(params.adaptive_dhat_ratio, ratio)
+
+    def test_adaptive_dhat_ratio_affects_adaptive_dhat(self):
+        """A larger ratio yields larger per-element dhat, so a deformed
+        configuration activates more collisions. Guards against the property
+        being stored but never reaching compute_adaptive_dhat()."""
+        mesh, rest = two_cubes()
+
+        # Move the right cube toward the left one so the deformed state is
+        # close enough to activate, but only for large enough adaptive dhat.
+        deformed = rest.copy()
+        deformed[rest[:, 0] > 0.96, 0] -= 0.04
+
+        counts = []
+        for ratio in (0.1, 0.5, 0.9):
+            params = ipctk.SmoothContactParameters(DHAT, 0.5, 0.0, 0.1, 0.0, 2)
+            params.adaptive_dhat_ratio = ratio
+            collisions = ipctk.SmoothCollisions()
+            collisions.compute_adaptive_dhat(mesh, rest, params)
+            collisions.build(mesh, deformed, params, True)
+            counts.append(len(collisions))
+
+        self.assertEqual(counts[0], 0, f"expected no activation, got {counts}")
+        self.assertLess(counts[0], counts[1], f"not monotonic: {counts}")
+        self.assertLess(counts[1], counts[2], f"not monotonic: {counts}")
+
+
+class TestSmoothCollisionsAdaptiveDhat(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls.mesh, cls.rest = two_cubes()
+        cls.params = ipctk.SmoothContactParameters(
+            DHAT, 0.5, 0.0, 0.1, 0.0, 2)
+        cls.potential = ipctk.SmoothContactPotential(cls.params)
+
+    def _build(self, use_adaptive_dhat):
+        collisions = ipctk.SmoothCollisions()
+        if use_adaptive_dhat:
+            collisions.compute_adaptive_dhat(self.mesh, self.rest, self.params)
+        collisions.build(self.mesh, self.rest, self.params, use_adaptive_dhat)
+        return collisions
+
+    def test_non_adaptive_has_spurious_rest_forces(self):
+        """Baseline: without adaptive dhat, this mesh/dhat pair produces a
+        nonzero potential at rest. This is what adaptive dhat exists to fix, so
+        if it ever becomes zero the test below stops proving anything."""
+        collisions = self._build(use_adaptive_dhat=False)
+        potential = self.potential(collisions, self.mesh, self.rest)
+        gradient = self.potential.gradient(collisions, self.mesh, self.rest)
+        self.assertGreater(potential, 0.0)
+        self.assertGreater(np.abs(gradient).max(), 0.0)
+
+    def test_adaptive_dhat_eliminates_spurious_rest_forces(self):
+        """GCP's 'no spurious forces' guarantee: with adaptive dhat computed
+        from the rest configuration, the potential and its gradient are exactly
+        zero in that configuration."""
+        collisions = self._build(use_adaptive_dhat=True)
+        potential = self.potential(collisions, self.mesh, self.rest)
+        gradient = self.potential.gradient(collisions, self.mesh, self.rest)
+        self.assertEqual(potential, 0.0)
+        np.testing.assert_array_equal(gradient, np.zeros_like(gradient))
+
+    def test_broad_phase_argument_accepted(self):
+        collisions = ipctk.SmoothCollisions()
+        collisions.compute_adaptive_dhat(
+            self.mesh, self.rest, self.params, ipctk.LBVH())
+
+
+class TestSmoothContactPotentialNaming(unittest.TestCase):
+    """The Python class was previously exposed as "SmoothPotential", which did
+    not match the C++ name. Guard the rename in both directions."""
+
+    def test_matches_cpp_name(self):
+        self.assertTrue(hasattr(ipctk, "SmoothContactPotential"))
+
+    def test_old_name_removed(self):
+        self.assertFalse(hasattr(ipctk, "SmoothPotential"))
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/src/ipc/ccd/additive_ccd.cpp b/src/ipc/ccd/additive_ccd.cpp
index b0e536ccb..0cd0abf21 100644
--- a/src/ipc/ccd/additive_ccd.cpp
+++ b/src/ipc/ccd/additive_ccd.cpp
@@ -68,11 +68,11 @@ AdditiveCCD::AdditiveCCD(
     conservative_rescaling = _conservative_rescaling;
 }
 
+template <typename DistanceSqrFunc>
 bool AdditiveCCD::additive_ccd(
     VectorMax12d x, // mutable copy
     Eigen::ConstRef<VectorMax12d> dx,
-    const std::function<double(Eigen::ConstRef<VectorMax12d>)>&
-        distance_squared,
+    const DistanceSqrFunc& distance_squared,
     const double max_disp_mag,
     double& toi,
     const double min_distance,
diff --git a/src/ipc/ccd/additive_ccd.hpp b/src/ipc/ccd/additive_ccd.hpp
index 875544c22..412758671 100644
--- a/src/ipc/ccd/additive_ccd.hpp
+++ b/src/ipc/ccd/additive_ccd.hpp
@@ -137,11 +137,11 @@ class AdditiveCCD : public NarrowPhaseCCD {
     /// @param min_distance The minimum distance between the objects.
     /// @param tmax The maximum time to check for collisions.
     /// @return True if a collision was detected, false otherwise.
+    template <typename DistanceSqrFunc>
     bool additive_ccd(
         VectorMax12d x, // mutable copy
         Eigen::ConstRef<VectorMax12d> dx,
-        const std::function<double(Eigen::ConstRef<VectorMax12d>)>&
-            distance_squared,
+        const DistanceSqrFunc& distance_squared,
         const double max_disp_mag,
         double& toi,
         const double min_distance = 0.0,