# Convergent Formulation¶

In addition to the original implementation of Li *et al.* [2020], we also implement the convergent formulation of Li *et al.* [2023].

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`

.

```
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);
```

```
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 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¶

*We briefly summarize the convergent formulation here for convenience.*

In order to derive a convergent formulation, we first define a continuous form of our barrier potential \(P\). For a surface \(\mathcal{S}\) embedded in 3D space, we parameterize the surfaces by common (possibly discontinuous) coordinates \(u \in \tilde{M} \subset \mathbb{R}^2\), so that \(\mathcal{S}(u)\) traverses the material points across all surfaces contiguously. The total barrier potential is then

where we define the operator \(\setminus_r: \mathcal{P}(\mathbb{R}^2) \times \mathbb{R} \times \mathbb{R}^2 \mapsto \mathcal{P}(\mathbb{R}^2)\) to be

with \(r \rightarrow 0\).

We then define our surface discretization with a triangulated boundary mesh geometry. As in the smooth case, we can parameterize the domain across all triangles with \(u \in \tilde{M}\) so that \(p(u): \tilde{M} \mapsto \mathbb{R}^3\) traverses all material points, across all triangles \(t ∈ T\) in the triangle mesh contiguously. The corresponding surface barrier potential is then

where \(T \setminus p\) is the set of boundary faces that do not contain the point \(p\).

Applying mesh vertices as nodes (and quadrature points), we numerically integrate the surface barrier potential. For each nodal position \(x \in V\) we then have a corresponding material space coordinate \(\bar{x} \in \bar{V}\). Piecewise linear integration of the surface barrier is then

where \(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 \(\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 \(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

where \(V_{\text{int}} \subseteq V\) is the subset of internal surface nodes and \(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

where we simplify with \(w_x = w_{\bar{x}}\) defined appropriately, per domain, as covered above.

Please, see the paper for more details (including the formulation for 2D curves and edge-edge collisions) and evaluation.

The key difference between the original and the convergent formulations is that we (1) include area weights in the barrier potential and (2) include additional (negative) terms to cancel out the duplicate distance pairs. While this requires minor algorithmic changes, it produces considerably better results.

## Physical Barrier¶

We want our barrier potential to have the same units as our elastic potential (e.g., \(\text{J}\)). Together with the area weighting (discussed above), this means the barrier should have units of pressure times distance (e.g., \(\text{Pa} \cdot \text{m}\)). That is,

To achieve this, (when using the convergent formulation) we modify the barrier function to have units of distance:

Note

This is equivalent to the original barrier function of [Li *et al.*, 2020] times \(1/\hat{d}^3\) when using squared distances. Therefore, to simplify the implementation we only implement the original barrier function and multiply all barrier potentials by \(1/\hat{d}^3\).

The barrier stiffness (\(\kappa\)) then has units of pressure (e.g., \(\text{Pa}\)), the same as Young’s modulus (\(E\)) in elasticity. This implies we can get good solver convergence even when using a fixed \(\kappa\) by setting it relative to the material’s Young’s modulus (\(\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`

.

## Friction¶

Just as with the collisions, we implement both the original friction formulation of Li *et al.* [2020] and the convergent formulation of Li *et al.* [2023].

The choice of formulation is dependent on how the fixed set of `collisions`

given to `TangentialCollisions::build`

was built. If the `collisions`

were built using the convergent formulation, then the friction collisions will also use the convergent formulation. Otherwise, the original formulation will be used.