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

\[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 \(\setminus_r: \mathcal{P}(\mathbb{R}^2) \times \mathbb{R} \times \mathbb{R}^2 \mapsto \mathcal{P}(\mathbb{R}^2)\) to be

\[\tilde{M} \setminus_r u:=\left\{v \in \tilde{M} \mid\|u-v\|_2>r\right\}\]

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

\[\frac{1}{2} \int_{u \in \tilde{M}} \max_{t \in T \backslash p(u)} b(d(p(u), t), \hat{d})~\mathrm{d} u,\]

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

\[\frac{1}{2} \sum_{\bar{x} \in \bar{V}} w_{\bar{x}} \max _{t \in T \backslash x(\bar{x})} b(d(x(\bar{x}), t), \hat{d}),\]

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

\[\begin{split}\begin{aligned} \Psi_s(x) & =\sum_{t \in T \backslash x} b(d(x, t), \hat{d})-\sum_{e \in E_{\text {int }} \backslash x} b(d(x, e), \hat{d})+\sum_{x_2 \in V_{i n t} \backslash x} b\left(d\left(x, x_2\right), \hat{d}\right) \\ & \approx \max _{t \in T \backslash x} b(d(x, t), \hat{d}), \end{aligned}\end{split}\]

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

\[P_s(V)= \frac{1}{2} \sum_{x \in V} w_x \Psi_s(x),\]

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,

\[\text{Pa} \cdot \text{m} \cdot \text{m}^2 = \frac{\text{N}}{\text{m}^2} \cdot \text{m} \cdot \text{m}^2 = \text{N} \cdot \text{m} = \text{J}.\]

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

\[\begin{split}b(d, \hat{d})=\left\{\begin{array}{lr} -\hat{d}\left(\frac{d}{\hat{d}}-1\right)^2 \ln \left(\frac{d}{\hat{d}}\right), & 0<d<\hat{d} \\ 0 & d \geq \hat{d} \end{array}\right.\end{split}\]

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 FrictionCollisions::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.

Last update: Oct 08, 2024