GIPC: Fast and stable Gauss-Newton optimization of IPC barrier energy

TL;DR

The key idea is to formulate a simplicial geometric measure of contact using mesh boundary elements, from which analytic eigensystems are derived and enhanced with filtering and stiffening terms that ensure robustness with respect to the convergence of a Project-Newton solver.

Abstract

Barrier functions are crucial for maintaining an intersection and inversion free simulation trajectory but existing methods which directly use distance can restrict implementation design and performance. We present an approach to rewriting the barrier function for arriving at an efficient and robust approximation of its Hessian. The key idea is to formulate a simplicial geometric measure of contact using mesh boundary elements, from which analytic eigensystems are derived and enhanced with filtering and stiffening terms that ensure robustness with respect to the convergence of a Project-Newton solver. A further advantage of our rewriting of the barrier function is that it naturally caters to the notorious case of nearly-parallel edge-edge contacts for which we also present a novel analytic eigensystem. Our approach is thus well suited for standard second order unconstrained optimization strategies for resolving contacts, minimizing nonlinear nonconvex functions where the Hessian may be indefinite. The efficiency of our eigensystems alone yields a 3x speedup over the standard IPC barrier formulation. We further apply our analytic proxy eigensystems to produce an entirely GPU-based implementation of IPC with significant further acceleration.

Figures (21)

• Figure 1: The ideal shape subtended by vertices $\mathbf{x}_{A}$, $\mathbf{x}_{B}$ and $\mathbf{x}_{\hat{d}}$ is deformed into the current shape (with $\mathbf{x}_d$ instead of $\mathbf{x}_{\hat{d}}$) due to proximal contact where $d < \hat{d}$. Deformation is understood to be in the direction $-\mathbf{n}(\mathbf{x})$, which is opposite to the normal vector that is computed from vertex positions $\mathbf{x}$ in the stencil of the local contact. The positions of the vertices in the stencil here are $\mathbf{x}_{A}$, $\mathbf{x}_{B}$ and $\mathbf{x}_{{d}}$, which come from the mesh boundary elements at the current step. One can always construct the ideal/undeformed shape of the simplex by duplicating a subset of the vertices $\in \mathbf{x}$ and shifting them by an offset $\mathbf{s} = (\hat{d}-d)\mathbf{n}(\mathbf{x})$ from the position of their original copy. Only $\mathbf{x}_d$ is duplicated and shifted by the offset in the example illustrated here. The linear map $\mathbf{J}(\mathbf{x})$ is then a transformation from the ideal to the current shape, which is a pure compression (without shear) due to the uniaxial change-in-shape with $\mathbf{s}$.
• Figure 2: A plot of the norms of the barrier gradient and Hessian together with their ratio for illustrating the source of vanishing displacements $\mathbf{d}$ as $g \rightarrow 0$. The scalar expressions used to measure the respective norms are described in Appendix \ref{['app:magnitude_evaluation']}. A threshold of $\hat{d} = 1$ is assumed here to aid visualization.
• Figure 3: Barrier function comparison, where a threshold of $\hat{d} = 1$ is assumed to aid visualization. Note that the form $b(d^2, \hat{d}^2)$ is actually used in the source code of 10.1145/3386569.3392425 to avoid squared roots, which results in an energy with a profile that is identical to Eq. (\ref{['eq:ivc-barrier']}).
• Figure 4: Plot of the eigenvalues of the approximate mollified Hessian in Eq. (\ref{['eq:ee-barrier-hessian']}): Eigenvalues of the first term ( Fig. \ref{['fig:lambda-alpha1']} and Fig. \ref{['fig:lambda-alpha23']}); Eigenvalues of the second term ( Fig. \ref{['fig:lambda-beta1']} and Fig. \ref{['fig:lambda-beta23']}); Eigenvalues due to the third and forth term (Fig. \ref{['fig:lambda-prime7']} and Fig. \ref{['fig:lambda-prime8']}).
• Figure 5: Evaluating the entire force-Jacobian in Eq. (\ref{['eq:hess-tens']}) versus evaluating only the first term. The latter (Fig. \ref{['fig:search-dir-ours']}) restricts the solution vector $\mathbf{d}$ (visualized by the red arrows) to the contact normal direction, and the former (Fig. \ref{['fig:search-dir-li']}) is affected by higher order variations between contact pairs due to the nonlinearity of Euclidean distance. The influence of these high-order variations on $\mathbf{d}$ are similarly observed with the distance based formulation of 10.1145/3386569.3392425. Refer to our supplementary video.