Convergent Formulation

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

To enable the convergent formulation, we need to set use_convergent_formulation before building Constraints:

collision_constraints.set_use_convergent_formulation(true);
collision_constraints.build(collision_mesh, vertices, dhat);

collision_constraints.use_convergent_formulation = True
collision_constraints.build(collision_mesh, vertices, dhat)

Important

The variable use_convergent_formulation should be set before calling CollisionConstraints::build for it to take effect. By default, it is 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 contact 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 contact 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 contact 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}\).

We next need to smoothly approximate the max operator in the contact 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 contact 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.

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 contact) 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 constraints 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”.

Friction

Just as with the collision constraints, 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 contact_constraints given to FrictionConstraints::build was built. If the contact_constraints were built using the convergent formulation, then the friction constraints will also use the convergent formulation. Otherwise, the original formulation will be used.