The Incremental Potential Contact (IPC) method enables robust complex simulations of deformable objects with contact and friction.  The key to IPC's robustness is its strict adherence to geometric constraints, avoiding intersections, which are a common cause of robustness issues in contact mechanics. A key element of the IPC approach to contact is a geometric barrier function, which is defined directly in the discrete setting.
While IPC achieves its main goal of providing guarantees for contact constraints,  its parameters need to be chosen carefully to avoid significant simulation artifacts and inaccuracies.
We present a systematic derivation of an IPC-like continuum potential defined for smooth and piecewise smooth surfaces, starting from identifying a set of  natural requirements for contact potentials, including the barrier property, locality, differentiable dependence of shape, and absence of forces in rest configurations, based on the idea of candidate sets.  Our potential is formulated in a way independent of surface discretization.

This new potential is suitable for piecewise-linear surfaces and its efficiency is similar to standard IPC. We demonstrate its behavior and compare it to IPC on a range of challenging contact examples.

Spurious stresses: nonuniform 3D mesh Incremental potential contact <cit.> introduces spurious contact forces in many configurations. We show a 3D block with rounded corners, maximal edge length $\qty{0.42}{\m}$, and minimum edge length $\qty{0.03}{\m}$ Left: With $\hat{d}$ chosen suitably small for the coarsely discretized regions (in this example $\hat{d} = \qty{0.07}{\m}$,), IPC introduces spurious forces at finely meshed corners, leading to jagged deformation artifacts. Right: our method maintains the rest configuration in the absence of other outside forces by deriving a new IPC-like contact formulation from the continous case.

§ INTRODUCTION

Contact modeling is a critical component of simulation tools for many domains, including physical simulation for computer graphics, robotics,  mechanical design, and biomedical engineering.  Recent advances, in particular, IPC <cit.>, make it possible to perform, in a reliable way, complex, large-deformation, and long-duration simulations of deformable objects with self-contact, with minimal parameter tuning.

The fundamental principle of the IPC solvers is the strict adherence to geometric constraints over the entire trajectory (no intersections including self-intersections), even when approximations are made in other aspects of the problem.  This feature is the key to robustness: overwhelmingly, robustness issues in contact mechanics are geometric in nature; once a contact constraint is violated, it may be very difficult for the simulation to recover.

Geometric constraint satisfaction in IPC is accomplished by a combination of using a geometric barrier function, i.e. an interaction potential that approaches infinity for contact points, and conservative CCD. IPC formulation excels at its main goal, qualitatively improving the robustness of contact simulations.

At the same time, the approach has limitations, and certain conditions may lead to inaccurate results and simulation artifacts.

The primary reason for this behavior is several specific features of the design of the barrier potentials used in the simulation.

The potential was designed directly for discretized models, rather than obtained from a discretization of a continuum potential; as we discuss in more detail in <ref>. In particular, there is a direct, and strong limitation on how local the potential is (and, as a consequence, how stiff it is) and how the object is discretized.

In this paper, we undertake a systematic derivation of a continuum barrier potential and its discretization from a set of desiderata, elaborated in <ref>. Similarly to Convergent IPC <cit.>, we assume that the total contact potential is obtained by integration over the surface of a pointwise potential $\psi(x)$. The following requirements, which we make more precise in <ref> include:

* : The total potential over the surface is well-defined for any piecewise smooth surface not in contact.

* : The potential $\psi(x)$ grows to infinity as the objects in the simulation approach (self-)contact at $x$, but is finite for any configuration not in contact.

* : In the undeformed configuration the potential is zero.

* : The potential has a parameter $\hat{d}$ and vanishes if the objects are further away than $\hat{d}$ from contact.

* : The potential depends smoothly on the surface configuration (e.g., for piecewise surfaces, on mesh vertex positions).

* : The localization of the contact potential can be controlled independently from mesh resolution.

Note that two key requirements (and, less directly, depend on the notion of distance to contact: how we measure this distance for self-contact is critical for choosing the right potential <ref>.

The original IPC formulation <cit.> meets these requirements only partially for piecewise linear surfaces.

Namely, it can violate and depending on the choice of $\dhat$.
To address this we introduce an IPC-like potential designed to satisfy these requirements. The key features of our approach include:

* The concept of contact candidate sets, that, in turn, supports a definition of distance to contact suitable for contact barrier potentials.

* An approach for defining a contact potential for piecewise smooth surfaces
that meets the requirements enumerated above,

independently of the choice of discretization.

* A discrete version of this potential for piecewise-linear surfaces retaining the same properties and efficiency of standard IPC.

We demonstrate in a set of benchmarks that satisfying the properties above is crucial for accurate elastodynamic: our formulation avoids spurious, discretization-dependent forces introduced by the original IPC formulation <cit.> (<ref>). The decoupling of the locality from mesh resolution further reduces the parameter tuning required for scenes with large deformations (<ref>).

§ RELATED WORK

Traditional methods for contact response mostly formulate non-interpenetration constraints based on the normal displacements of the contacting regions with locally approximated geometries. Then the constrained time integration problem is solved with Lagrange multipliers <cit.>, penalty methods <cit.>, barrier methods <cit.>, etc. To guarantee global convergence for the solution of frictionless large deformation contact problems for hyperelastic materials, <cit.> proposed a filter-trust-region method to solve the constrained minimization problem formulated from implicit time integration.

Mortar methods <cit.> are commonly used in engineering <cit.> and biomechanics <cit.> but require to (a priori) mark all possible pairs of contacting mesh boundary. A clear limitation of these methods is then that they cannot handle, without extensions, self-collisions or collisions in regions with more than two contacting surface regions. <cit.> provide a didactic comparison of IPC and one such mortar method <cit.>. There it can be seen that such methods enforce contact constraints weakly and therefore allow intersections (especially at large timesteps and/or velocities).

Instead of directly treating contact as unilateral constraints, <cit.> and <cit.> formulate frictionless contact based on distance-dependent surface interaction potentials, providing a unified framework for various contact formulations including Lagrange multiplier, penalty, and barrier. Starting with the first and second laws of thermodynamics, <cit.> recently extended <cit.> to account for Coulomb friction under isothermal conditions, and demonstrated the convergence of the method under refinement.

Gap functions Gap functions compute a signed distance for pairs of points or primitives, where negative distance measures the interpenetration depth. Collisions are then resolved by enforcing that all gap functions are non-negative. Commonly used gap functions include projecting between pairwise surface primitives onto a fixed geometric normal <cit.> or using CCD to determine the point of contact and compute a projected distance along the contact normal <cit.>. <cit.> introduced the method of space-time interference volumes which was in part motivated by the need to define a version of the gap function in a fully general setting; the complexity of their construction shows that even rigorously formulating the self-contact problem in inequality-constrained is nontrivial. We discuss the connection of our method to one type of gap function  in Section <ref> in more detail.

Barrier-based methods The use of barrier functions for handling contact has been utilized in graphics for over a decade now. Most notably, <cit.> and <cit.> utilize a set of layered discrete penalty barriers that grow unbounded as the configuration approaches contact. However, this incremental construction makes it unsuitable for optimization-based implicit time integration and therefore requires small time-steps for stability.

<cit.> simulate knitted cloth at the yarn level

where yarn-yarn collisions are handled through a continuous potential integrating a barrier function over two disjoint spline segments. This approach, however,  ignores self-collisions within a single segment and between neighboring segments. To discretize their model, <cit.> use “Simpson’s quadrature at fixed positions in parameter space” which has the limitation of potentially missing the nearly singular portion of the integral and therefore introducing intersections.

Our method largely follows the IPC method of <cit.>, which transforms discrete gap constraints into a sum of log barrier functions. The method of <cit.> has been directly utilized and extended by numerous follow-up works to include support for codimensional elements <cit.>, rigid/affine body dynamics <cit.>, medial elastics <cit.>, solid-fluid interactions <cit.>, higher-order finite element analysis <cit.>, etc.

The most significant limitation of IPC <cit.> is its lack of convergence under mesh refinement due to its purely discrete formulation (see <cit.> and <cit.> for details).

To address this, <cit.> form a continuous formulation for the original IPC model which requires minimal adjustments to existing implementations. The  limitation of their approach is the requirement for joint refinement of $\hat{d}$ and mesh resolution. We provide further discussion of the limitations below.

<cit.> show the tendency of “Node-to-Segment” methods (e.g., IPC <cit.>) to produce spurious tangential contact forces. To remedy this they propose to use a fully $C^1$-continuous surface representation; computing distances used in the contact potential between the nodes and smooth surface (extending the method of <cit.>). The main limitation of <cit.>, however, is the inability to handle self-collisions.

Tangent-point energies
A possible continuous approach to avoiding intersections are the tangent-point energies <cit.> which repel surface points by integrating over all surface pairs $(x,y) \in \boundary \times \boundary$:

\[
\int_{x\in\boundary} \int_{y\in\boundary} \frac{1}{r(f(x),f(y))^\alpha}~\mathrm{d}x\mathrm{d}y,
\]

where $r$ is the radius of the smallest sphere tangent to $f(x)$ and passing through $f(y)$ and $\alpha$ controls the repulsion strength. This avoids having infinite energy for geodesic neighbors as the radius is large, but grows to infinity for points approaching zero distance in Euclidean space.
Recently, <cit.> utilize tangent-point energies to prohibit intersections while optimizing surface geometry. In pure form, this method is unsuitable for collision response due to its global coupling of all points, but it is closely related to our approach as discussed in Section <ref> in the case of smooth surfaces.  This not only leads to a dense system that needs to be factorized in the optimization but also spurious contact forces between distant pairs.

§ CONTACT POTENTIAL REQUIREMENTS

We consider a collection of deformable objects defined on a material domain $\domain$, whose boundary is $\boundary$. $\domain$ may have multiple connected components (<ref>).

An admissible deformation
$f: \domain \mapsto \mathbb{R}^n$, $n=2,3$ is noninjective on the boundary only, i.e., $f$ is injective in
the interior of the domain, but may have contact points
$x \neq y$ on the boundary $\boundary$, for which $f(x) = f(y)$.

A collection of objects is transformed using a deformation $f$ (left), which can cause the objects to intersect at a contact point (right).

We say that an admissible $f$  is a contact deformation, if it has contact points, otherwise, we call it contact-free.

We assume that for a point $x \in \boundary$, we have a measure $d_c(x,f)$ of how far it is from being a contact point. Defining this measure is a key aspect of our construction, and is necessary to formulate our contact requirements.

We note that the problem of defining $d_c(x,f)$ is far from trivial. While for rigid objects the distance to the closest point on another object is an adequate measure, this is not the case for deformable objects. Due to the possibility of self-contact, we cannot exclude points on the same object and, in this case, there are always points on $\boundary$ arbitrarily close to $x$ in Euclidean distance that need to be considered to be far from being in contact with $x$.  At the same time,

an arbitrarily small, if measured by Euclidean distance, perturbation of a surface may create an actual contact, so the distance from $x$ to a real contact point $y$ can be arbitrarily small.

Before describing our approach to the distance to contact, we make the properties we stated in the introduction more precise.

A total contact potential
\[\Psi(f) = \int_\boundary \potential(x;f)~\mathrm{d}x\]
is well-defined, i.e., $\potential(x;f)$ is integrable on $\boundary$ for a sufficiently broad class of surfaces.  Our construction applies to piecewise-smooth surfaces, although in our implementation we focus on piecewise linear.

$\potential(x;f)$ increases to infinity for $f$  if the distance to contact $d_c(x,f)$ goes to zero.  Combined with incremental potential time-stepping and CCD, this can be used to guarantee that all configurations remain contact-free.

For the undeformed configuration $f = \operatorname{Id}$,  $\potential(x;\operatorname{Id})$ is zero everywhere.  This is necessary to ensure that the potential does not create artificial forces that would cause motion/deformations even if no external forces are acting on the object.

$\potential(x;f)$ has a locality parameter $\epsilon > 0$ ($\dhat$ in the notation of <cit.>), with no restrictions on its magnitude, on which it depends at least continuously, and the potential vanishes if $d_c(x,f) > \epsilon$. If we solve a sequence of problems with decreasing $\epsilon$, we approach a solution of the standard inequality-constrained formulation of contact problems.

If $f$ is defined by a finite number of parameters (in the simplest case, vertex positions of a mesh) then $\potential(x;f)$  depends differentiably, and piecewise twice differentiably on the parameters of $f$.  Then the potential leads to a force with piecewise continuous Jacobian, allowing for second-order methods for implicit time-stepping.

Importantly, these properties are independent of whether we consider a smooth problem or its discretization: below, we require that both a smooth potential and its discretization satisfy these properties.

In the discrete case, the total potential is a sum of discrete potentials defined on pairs of discrete primitives. We note that it is important that all requirements above are satisfied by the discretized total contact potential directly, not only in its continuous formulation.

§ FORMULATION

Our approach to defining contact potentials is conveniently formulated in terms of contact candidate sets $C(x; f)$. These sets are defined for every point $x$ as a subset of $\boundary$, and the potential we construct is restricted to these.

To be able to use them in the definition of a contact potential, we choose the contact sets $C(x)$ with the following properties

* For $f$ in contact, $C(x)$ contains the contact point of $x$;
* If is not in contact, the (Euclidean) distance from $f(x)$ to $C(x)$ is positive, and it decreases to zero if we interpolate between an $f$ in contact and $\tilde{f}$

in contact.

We define the distance to contact $d_C(x,f)$ of a point $x$ under deformation $f$ as the distance to the contact candidate set of $x$:

\[
d_C(x,f) = \min_{y \in C(x)}\|f(x)-f(y)\|.
\]

If a potential $\psi(x;f)$ is obtained by integrating a positive function $g(x,y)$ supported on the contact set, depends on the distance between $x$ and $y$, and increases to infinity if the distance is zero, then we can observe that using such candidate sets allows us to satisfy <ref>.

In particular, the <ref> () is satisfied if $g(x,y)$  is chosen to vanish for a small $d_C(x,f)$, and <ref> () is satisfied if it vanishes for $d_C(x,Id)$.

§.§ Rigid objects

In the rigid case, for a point $x$, the contact candidate set (highlighted in blue) is the set of points in all connected components of $\boundary$ that do not contain $x$.

To make the contact set definition less abstract, we briefly consider the case of rigid objects, for which $f$ restricted to each connected component of $\boundary$ is rigid motion, as a guiding example for defining a contact candidate set and a corresponding potential satisfying the contact potential requirements (<ref>).

In this case, the choice is simple and intuitive: for $x$, the contact candidate set $C(x)$ is the set of points in all connected components of $\boundary$ that do not contain $x$  (<ref>). Then the distance $d_c(x)$ is simply the distance to the closest object. The requirements above for the contact set are satisfied, and one can define the potential as

\begin{equation}
\psi^R(x; f)  = \int_{y \in C(x)} \pex(\|f(x)-f(y)\|)~\mathrm{d}y
\end{equation}

where $\pe(z) \rightarrow \infty$ as $z \rightarrow 0$, fast enough for the integral to be unbounded if $x$ is in contact and

there are points in $C(x)$ arbitrarily close to $x$ (<ref>).  We choose

\[\pe(z) = \frac{\he(z)}{z^{n-1}}\]

for $n=2,3$,

\[
\he(z) = \frac{3}{2} B^3\left(2z/\epsilon\right)
\]

where $B^3\in C^2(\mathbb{R})$ is a cubic spline (<ref>)

\[
B^3(v) = \begin{cases}
{\frac{2}{3}}-{{v}^{2}+\frac {1}{2}|v|^{3}} & |v|<1\\
{\frac {1}{6} \left( 2-|v| \right)^3} & 1 \leq |v|<2 \\
0 & 2 \leq |v| \\
\end{cases}.
\]
Therefore, $\he \in C^2(\mathbb{R})$ satisfies

\[\he(\epsilon) = \he'(0) = \he''(\epsilon) = 0.\]

Plot of $B^3(z)$.

We note that somewhat counter-intuitively, it is not strictly necessary for the potential $\psi(z)$ to be infinite to guarantee that no contact happens (see <ref>).
The choice of the power for the potential is important: on the one hand, e.g., the logarithmic potential used in IPC, in the continuum case may result in a finite integral, i.e., requirement <ref> not satisfied; on the other hand, if the power is too high, that the pointwise potential may lead to divergent integrals for piecewise-smooth surfaces folding along sharp features (see Section <ref>).

Barriers Left: plot of the cubic spline $\he(z)$ and barrier function $\barrier_{\epsilon}(z)=\he(z)/z^{n-1}$ for $\epsilon=1$. We also include a plot of the log barrier $\pipc(z)$ of <cit.> for comparison. Right: our barrier improves approximation to the discontinuous function as $\epsilon$ goes to 0.

Then $\he$ vanishes if $d_C(x;f)$ exceeds $\epsilon$, and as a consequence, $\psi^R$ is zero. By construction, $\psi^R(x; f)$ becomes infinite if $d_C(x;f) \rightarrow 0$, i.e., this potential satisfies <ref>, as long as the potential grows fast enough. To satisfy <ref> (), we choose a variable $\epsilon(x)$ for each point $x$ to be equal to

\[\min( d_C(x,f_0)/2, \epsilon_\text{trg}),\]

where $\epsilon_\text{trg}$ is the global parameter determining the maximal distance at which the potential may be nonzero.

Note that $\epsilon(x)$ is not necessarily smooth or even continuous as a function of $x$.  This may affect the convergence of the outer integral in computing $\Psi$; however, it does not affect the required properties, in particular, <ref> () of the total potential's dependence on the shape parameters of $f$.

§.§ Deformable smooth surfaces

Contact candidate sets for self-contact
For an arbitrary deformation $f$, self-contact needs to be taken into account since points on the same objects cannot be excluded. We need a different approach for defining $C(x)$.

Our key observation is:
for a fixed point $x$ and $y \neq x$ in contact, if we consider  $d(y) = \|f(x) - f(y)\|$ with $y$ varying over the surface, $d(y)$ has a local minimum at $y$ (<ref>).

Clearly, $d(y)$ has also a local minimum at $x$ which we exclude, but not at any point in a neighborhood of $x$.  This leads to the idea of including in $C(x)$ points that are close to being local minima of the distance function $\|f(x) - f(y)\|$ in the candidate set, and, in addition, the vector $f(x)-f(y)$ points towards the exterior of the surface at $f(x)$ and towards the interior at $f(y)$.  The latter is critical for handling thin objects.

We first consider smooth surfaces, to describe the main ideas without considering many cases needed for piecewise smooth surfaces, and then extend the formulation to a general case.

In this case, the local minima of the distance function from $f(x)$ to $f(y)$ for a contact-free $f$  satisfy the condition

\[\Phi(x,y)^2 = \|(f(y)-f(x))_+ \times n(x)\|^2  =  1- ((f(y)-f(x))_+ \cdot n(x) )^2 = 0,\]

where $(\cdot)_+$ denotes normalization to a unit vector.
The resulting distance to contact and potentials are related to one of the commonly considered gap function types in the FEM literature <cit.> and the potential to tangent-point kernels <cit.>.

Additionally, define

\[
\Phi^n(x,y) = n(x) \cdot ( f(y) -f(x))_+
\]

Then $f(y)$ is on the side of the outward pointing normal from $x$, if $\Phi^n(x,y) > 0$.

This suggests the following definition of a contact candidate set for a smooth surface:

\begin{equation}\label{eq:csmooth}
C(x) = \{ y \,|\;  \Phi^n(x,y) \geq -\alpha, \Phi(y,x)  \leq \alpha \}
\end{equation}

(<ref>). Observe that $C(x)$ for smooth surfaces is isolated from $x$:
as $y$ approaches $x$, $f(y)-f(x)$ becomes tangent to the surface and, as a consequence, perpendicular to the normal $n(y)$, so the distance from $C(x)$ to $x$ is positive for any $\alpha < 1$.

In the case of smooth surfaces $\Phi(x,y)^2 + \Phi^n(x,y)^2 = 1$, so restricting $\Phi(x,y)^2$ to be close to zero is equivalent to restricting $\Phi^n(x,y)^2$ to be 1. However, the relationship is more complex for piecewise smooth surfaces, so we treat these separately.

At the same time, as $f$ approaches contact at $x$, the set will eventually contain the contact point $y$.

Potential distribution on the surface with respect to the red dot. From left to right are potential distributions. (A) Without constraint on surface tangent, the potential distributes spherically around the red dot regardless of the surface shape. (B) With tangent constraints, only surfaces that are close to the local minima of distance from the red dot have high values. (C) With both tangent and normal constraints, only surfaces on the closer side of the volumetric object have high values.

In this case, the property <ref>, i.e., the fact that contact potential is well-defined for a contact-free $f$: while $\psi(x)$ may be defined for any point $x$,
it may have no upper bound, so the integral over $x$ is not defined. However,
there is a global positive lower bound on the separation distance for non-empty $C(x)$: the set of such points is closed, and the bound is positive for all points.  This ensures that <ref> () is satisfied for smooth surfaces.

Contact potential. Using a $C^2$ smooth function

* has a maximal value at zero,

* decays to zero for $|z| = \alpha$ (i.e., on the boundary of the contact candidate set), and

* integrates to 1,

we also use the smoothed Heaviside  step function $H^\alpha(z) = H(z/\alpha)\in C^2(\mathbb{R})$, with

\[
H(z) = \begin{cases}
0 &  z < -3\\
\frac{1}{6}(3+z)^3 & -3 \leq z<-2\\
\frac{1}{6}(3-9z-9z^2-2z^3) & -2 \leq z<-1\\
1+\frac{1}{6}z^3 & -1 \leq z<0\\
1 & 0 \leq z
\end{cases};
\]

Smoothed Heavyside function $H^{\alpha}(z)$ with $\alpha=\tfrac{1}{2}$.

Then we define the angle factor as

\begin{equation}\label{eq:smooth-angle-factor}
\gamma^S(x,y) = \delta_\alpha(\Phi(x,y)) H^{\alpha}(\Phi^n(x,y))
\end{equation}

which ensures that the potential is supported within the contact set $C(x)$.

We can define a potential for deformable contact as

\begin{equation}
\Psi(f) =
\int_{x \in \boundary} \int_{y \in \boundary} \gamma^S(x,y)\gamma^S(y,x)  p_{\epsilon(x)}(\|f(x)-f(y))\|).
\label{eq:potential-deform}
\end{equation}

Here, similar to the rigid motion case, we choose

\[\epsilon(x) = \min( d_C(x; f_0)/2, \epsilon_\text{trg}),\]

where $\epsilon_\text{trg}$ is the maximal distance at which the potential is non-zero so that in the undeformed state, the potential vanishes everywhere.

Relation to gap functions and tangent kernels. One of the gap functions described in FEM literature (e.g., <cit.>)  is obtained by shooting a ray along the normal direction $n(x)$ to obtain a point $y(x)$, defining in this way the correspondence between surfaces, and defining the gap function as $n(x)\cdot (f(y(x))-f(x))$.  This is similar to applying $\Phi(x,y)$  for a specific choice of $y$. We note that in general geometry, $y(x)$ may not exist, or move discontinuously.  Similarly, repulsive surfaces <cit.> use potentials of the form $(n(x)\cdot (f(y)-f(x)))^p/\|x-y\|^{2p}$. Our potential can be considered a localization of repulsive potentials of this type in space and angle.

§.§ Piecewise smooth contact

We now extend the potential to piecewise smooth surfaces. There are two reasons for this extension: (1) piecewise smooth surfaces are ubiquitous in graphics and scientific computing as they can represent shapes with sharp features, and (2) smooth surfaces are often approximated by piecewise-linear for which efficient robust continuous collision detection is available.

We consider surfaces consisting of patches $\boundary_i$, which form a curved manifold mesh satisfying the standard definition <cit.>. We assume that each patch has continuously varying normals defined everywhere including boundary (i.e. no cones are allowed) and the edge curves meeting at a vertex have distinct tangents.
We refer to the curves separating patches as edge curves, or simply edges, and points shared by more than 2 patches
as vertices.

In this case, there are six possibilities for
a contact point, corresponding to the possible pairs of contacts between any two element types: Face-Face, Face-Edge, Face-Vertex, etc.

The smooth-surface definition can be applied in cases involving a face, as there is a well-defined normal on one side, and the definition of the potential requires only the normal at $y$. However, cases that do not involve a face (for example, Edge-Edge) require special handling and a generalization of our definition of $C(x)$ and distance to contact. Observe that we cannot
assume that, e.g., nearby face-vertex contacts are sufficient: two edge points may be close, but all points may be far from contact in the sense of $d_C(x; f)$ defined above.

Contact candidate sets for piecewise-smooth surfaces

To generalize $C(x)$ to piecewise smooth surfaces, we take the following approach. Consider the distance $\|f(y) - f(x)\|^2$ as a function of $y$.  Then  for the gradient w.r.t. $y$, we obtain

\[\nabla_y f(y)~(f(y)-f(x)) = 0\]

at the closest point.
The columns of the matrix $\nabla_y f$ are two tangents at $y$, i.e., this condition is equivalent to the condition of $n(y)$ being parallel to $f(y)-f(x)$ in the definition of $C$.

Contact point as a local minimum of distance.

To generalize the minimum condition to arbitrary points of piecewise smooth surfaces,  consider the set of faces $\boundary_i$ containing a point $y$ (one face for face points, two for edge points, and any number for vertex points). For each $\boundary_i$ and each parametric direction at $y$, there is a well-defined tangent $t$.

The condition for a local minimum is that for any $\boundary_i$, and any tangent $t$,
the distance increases along the tangent, i.e.,

\[\forall t,\; \transpose{t} \nabla_y f(y)~(f(y)-f(x)) \geq  0,\]

where the gradient $\nabla_y$ is computed with respect to a parametrization of $\boundary_i$.

If $i > 1$, for any face $\boundary_i$, there are two unit tangent vectors  $t^1_i(y)$ and $t^2_i(y)$ along two edge curves of $\boundary_i$ meeting at $y$ (<ref>). Any tangent vector $t$ at $y$ for $\boundary_i$ is a linear combination of  $t^1_i$ and $t^2_i$ unless these are collinear; however, in general, one cannot infer that $f$ increases along all directions $t$ between $t^1_i$ and $t^2_i$ from increase along these two directions, unless the angle between them is strictly less than $\pi$, i.e., not for points on straight edges, or curved surfaces if the patch $\Omega_i$ boundary has a concave corner.

To treat all cases uniformly, we add a halfway vector $t^3_i$;  corresponding to the direction halfway between $t^1_i$ and $t^2_i$. Note that the issue described above doesn't happen to piecewise linear surfaces, so $t_i^3$ is not used in this case.

Then every unit tangent vector $t$ within $\Omega_i$ can be represented as a linear combination of  $t^1_i$ and $t^3_i$, or $t^3_i$ and $t^2_i$, with non-negative coefficients, and consequently, positivity along these three directions is sufficient for all directions.

Denote $g(y) = \nabla_y f(y)~(f(y)-f(x))$. Then if $t^k_i \cdot g(y) \geq 0$ for $k=1,2,3$, the same holds for any $t = a_\ell t^\ell_i + a_m t^m_i$, $(\ell,m) = 1,2$ or $2,3$.

The only special case is the case of face points, which do not have any edge vectors; in this case, we can pick  $t^k_i$ as arbitrary tangent vectors spaced at equal 120-degree angles in the parametric domain. Note that this doesn't happen to piecewise linear surfaces.

The analog of $\Phi^n(x,y)$ requires a criterion for determining that a unit vector $v$ at a point $x$ of a piecewise surface points outside. Let $P_i$ be the tangent plane to $\Omega_i$,  and let $n_i$ be the normals.  Suppose for at least one face, the projection of the vector $v$ to $P_i$ is inside the face. Then the planar sector defined by the rays along $v$ and the closest projection $v'$, to the plane $P_i$ which has minimal value of $|v \cdot n_i|$, lies entirely on one side of the surface locally at $f(x)$.  Thus if $n_i \cdot v$ which is signed distance from the tip of $v$ to $P_i$ is positive, it is pointing to the outside of the surface.  Define  $d_i^+ = \min_i v\cdot n_i$, over $i$ for which $v \cdot n_i > 0$,  and $v$ projects inside the face $i$, setting it to $\infty$ if this set is empty.
Similarly define $d_i^- = \min_i -v\cdot n_i$ for $v \cdot n_i < 0$,
and similarly $v$ projects inside the face $i$, setting it to $\infty$ if the set is empty (both cannot be empty). Then $v$ points outwards if $d_i^- - d_i^+ > 0$.

This leads to the following definition of a contact candidate set:

2D version of notation for <ref>.

For each point $y$ of $\boundary$, let $I$ be the set of indices of patches $\boundary_i$ containing $y$, and  for each patch, define $t^k_i(y)$, $k=1,2,3$ as above.  Then $y \in C(x)$ for a point $y$ iff

\begin{equation}
\Phi_{ik}(x,y):=\transpose{t^k_i(y)} \nabla^i_y f (y)  (f(y)-f(x))_+  \geq \alpha
\label{eq:pw-contact-alpha}
\end{equation}

for $i \in I$, $k=1,2,3$, with $\nabla^i_y f$ computed with respect to $\boundary_i$, and $\alpha$ is a positive number, and
in addition, $d_i^- - d_i^+ \geq -\alpha$

Contact potential for p.w. smooth surfaces.  Using the new definition of contact candidate sets $C(x)$  we  can define  the contact potential in a similar way to smooth surfaces,  except we use a new angle factor, defined as

\begin{equation}\label{eq:angle-factor}
\gamma^{PS}(x,y) = \prod_{i=1}^{m_y}\prod_{k=1}^{3} H^{\alpha,b}(\Phi_{ik}(x,y))
\end{equation}

where $m_y$ is the number of patches meeting at a vertex, $m_y = 1$ for a face point $y$,
$m_y = 2$  for an edge point, and $m_y > 2$ for a vertex point. The function $H^{\alpha, b}(z)$ is a smoothed Heaviside function (<ref>) with shift, which is 1 for positive arguments $z \ge b$ and 0 below $-\alpha$:

H^{\alpha,b}(z):=H\left(\frac{3}{\alpha + b}(z - b)\right),

H^{\alpha,b}(b)=1,\ H^{\alpha,b}(-\alpha)=0, \text{and}\ H^{\alpha,b}\in C^2(\mathbb{R}),

Additionally, we define

\Phi_{i}^n(x,y)=\transpose{n_i(y)} \nabla^i_y f (y) (f(y)-f(x))_+

where $n_i(y)$ is the face normal of patch $\Omega_i$ at $y$, which is unique per patch on a piece-wise linear surface. To eliminate the contact between opposite sides of the object (Figure <ref>), the following should hold

\Phi_{i}^n(x,y) < 0 \quad \text{for some}\ i\in{1,\dots,m_y}.

Note that not all normals satisfy the constraint above when there is
a collision, so the constraint is on only one normal direction. While this condition for determining the correct side is precise, its mollification is relatively complex. For this reason, we relax this to a simpler, weaker condition, following a similar way of smoothing the constraints on $\Phi_{ik}$ and multiply $\gamma^{PS}(x,y)$ by

$$H^{1,0}\left(-1+\sum_{j=1}^{m_y} H^{\alpha,b}(-\Phi_{i}^n(x,y))\right),$$

which reaches maximum 1 when there exists $i$ such that

\Phi_{i}^n(x,y) \geq b,

and vanishes if for all $i$

\Phi_i^n(x,y) \leq -\alpha.

Since $\Phi_i^n,\ \Phi_{ik}$ are essentially dot products between unit vectors, $\alpha$ acts as a threshold of contact potential on the cosine value of the angle between contact pairs. One can pick $\alpha$ based on the desired angle threshold.

Recall that contacts between vertices and faces may not necessarily lead to non-zero potentials for faces.
As a consequence, the interaction potential in these configurations will vanish outside of a measure-zero set (edge curves and/or vertices), and direct integration of $\psi(x)$ will lead to a zero potential.

For this reason, we treat low-dimensional elements separately. Informally, we can think about the interaction of lower-dimensional elements, as expanding them into areas: e.g., an area of width $L$ along curves, and an area of size $L^2$ assigned to vertices, with the potential constant along the direction in the area orthogonal to the curve, or on the whole area assigned to a vertex.  More precisely, this corresponds to adding line integrals for edge curves weighted by $L$, and points sums for vertices, weighted by $L^2$.  The potentials we integrate for a pair of elements of $G$ and $H$, possibly of different dimensions, can be written in the form

\[
P(x,y; G, H) = \gamma^{PS}(x,y; G)\gamma^{PS}(y,x; H)\pe(x,y)
\]

where $\gamma^{PS}(x,y; G)$ is the angle factor defined in <ref>, for $x$ considered as a point on element $G$.  E.g., if $x$ is a vertex point on a boundary of a face $G$, then it is considered as a face point, with a single face used.  If $G$ however is an edge, both faces are used for $x$.

Then the total potential can be written as

\begin{equation}\label{eq:pwpot}
\sum_{ (g,h) }  L^{2-\dim g - \dim h} \sum_{i \in I_g, j \in I_h, G_i \neq H_j} \int_{G_i}\int_{H_j} P(x,y; G_i, H_j)
\end{equation}

where $g,h$ are one of $V, E, F$, the sum is over 6 unordered pairs  $I_g$ and $I_h$ are sets of indices of elements of types $g$ and $h$,  and the integrals are area, line, or 0-dimensional, i.e., simple evaluations for vertices.

If we define the distance to contact as

\[\min_{y \in C(x)} d(x,y),\]

then the potential <ref> satisfies all required conditions.

Parameters of the potential In principle, whether the potential satisfies contact potential requirements does not depend on the choice of  $L$, $\epsilon_\text{trg}$,

and $\alpha$, as long as these are positive.

* $\alpha$ determines the size of the contact candidate sets, and as a consequence how smoothly the potential depends on $f$:  the closer $\alpha$ is to zero, the less smooth the dependence is.

* $\epsilon_{trg}$ determines the upper bound on how far the potential extends ($\hat{d}$ in IPC) but especially at concave corners, $\epsilon_{trg}$ as well as $f_0$, affect this.

* $L$ determines the strength of potential for low-dimensional contact.

§.§ Discretization

To discretize these integrals, and still satisfy the contact potential conditions, we cannot directly apply an arbitrary quadrature: if we evaluate the potentials at arbitrary points on elements we may miss a contact point, and the potential will not satisfy the barrier property. To address this problem we choose the closest points as our quadrature points,  possibly reducing integration precision, but guaranteeing that the potential remains a barrier.

In this case, the lowest-order discretization of the potential integral is:

\[
\sum_{ (g,h)} \sum L^{2-\dim g - \dim h} \sum_{(i,j)} P(x_i,y_j; G_i, H_j) A(G_i) A(H_j)
\]

where $(x_i,y_j)$ is a pair of closest points on elements $G_i$ and $H_j$, respectively, and $(i,j)$ are pairs of non-adjacent
elements, and $A(G_i)$ and $A(H_j))$ are measures of elements (area for faces, length for edges and 1 for vertices).  Here the adjacent elements are dropped from summation.

For the discrete potential to depend differentiably on the shape parameters (vertex positions in the piecewise-linear case), we need an additional modification.  First, observe that the distance between two elements is already $C^1$ with the only exception being the distance between two parallel edges, which can be mollified as explained in <cit.>, so the term $\pe$ in $P$ does not require modifications.  However, the terms $\gamma(x,y)$ depend on the direction between closest points, which may change non-smoothly.

We use a mollification for the closest point position to obtain a smooth variation. Importantly, it does not affect the barrier property or the contact set properties.

Distance function mollification.

For an edge $[p_0,p_1]$, and a point $q$, we define the mollified closest point as follows. Let $0\leq s < 1$ be the coordinate of $q$ along the edge, and $0 \leq c < 1/2$ be the mollification parameter (we use $10^{-3}$). Our mollification is based on

0 & z \leq 0 \\
z(2-z) & 0 < z < 1 \\
1 & 1\leq z \\
\end{cases}\quad
h_c(s):= h(\frac{s}{c})

Define three basis functions
$m_0(s) = h_c(c-s)$, $m_1(s) = h_c(s-(1-c))$ and $m_2(s) = 1-m_0(s)-m_1(s)$, and mollified closest point as
$m_0(s_c)p_0 + m_1(s_c)p_1 + m_2(s_c)p_c$, where $p_c$ is the closest point, and $s_c$ its coordinate along the edge.   Then the $\gamma$ is computed using the mollified point.

Similarly for a triangle, we use barycentric coordinates $(s_1,s_2,s_3)$ to define a mollified closest point,
with  7 basis functions, one for each vertex and edge and the central one:

\begin{align*}
m^e_i(s_i) &=  m_0(s_i)(1-m_0(s_j))(1-m_0(s_k))\\
m^v_i(s_i) &=  (1-m_0(s_i)) m_0(s_j) m_0(s_k)
\end{align*}

where $i,j,k$ is a permutation of $(1,2,3)$ and $m^e_i$ and $m^v_i$ are the edge and vertex basis functions.
The mollified closest point combines the vertex, edge, and face closest points with these weights.

Note that with the amplification above, the edge-edge mollification in IPC is not needed anymore. The edge-edge mollification was introduced to avoid non-smooth distance with almost parallel edges. In our case, not only the parallel edges case but also the case where edge-vertex distance reduces to vertex-vertex distance is mollified.

Adaptive $\epsilon_{trg}$

We introduce adaptive $\epsilon_{trg}$ for our formulation. We first specify a fixed $\hat{d}_{\max}$ for the simulation, and collect contact pairs in the rest configuration within $\hat{d}_{\max}$ and with nonzero potential values. We pick $\epsilon_{trg}$ for every primitive (vertex/edge/face) so that none of the contact pairs are active. Then the $\epsilon_{trg}$ for every edge (face) is chosen to be the minimum among its neighboring vertices (edges). Benefited from our formulation with angle constraints, for a cube mesh there's no contact for a wide range of $\alpha$, so the adaptive $\epsilon_{trg}$ can be any positive number without causing spurious contact forces between its primitives.

Adjacent elements can be excluded.
Specifically, if adjacent elements are in contact, other, nonadjacent elements will be in contact.
Adjacent elements such that one does not contain another can be vertex-adjacent face and edge, and vertex- and face-adjacent faces.   Suppose an edge and a face are in contact, i.e., an edge is in the plane of the face, and they share a vertex and other points. Then the other edge endpoint is either inside the face, or the edge intersects the boundary of the face, crossing an edge which it does not share a vertex with.  In the first case, we have vertex-face contact, in the second case, we have non-adjacent edge-edge contact.   If it is a face-face contact, the argument reduces to the first, if applied to one of the face edges that is not a common edge with the other face.

Spurious stresses: nonuniform squircle mesh Our method properly handles nonuniform meshes. In this case, we consider a 2D rounded block (of size $\qty{1}{\m}$ x $\qty{1}{\m}$) with maximal edge length $\qty{0.21}{\m}$ in its straight sections and minimal edge length $\qty{0.01}{\m}$ at its corners. Top row: with $\hat{d}=\qty{0.1}{\m}$, IPC introduces spurious contact forces in the refined corners, resulting in a deformation. Bottom row: our method avoids this by calculating $\epsilon(x)$ such that the potential is zero in the rest configuration.

§ EVALUATION

We implement our algorithm in C++ by extending the IPC Toolkit <cit.> (details provided in the supplemental document). We use Eigen <cit.> for linear-algebra operations, Pardiso <cit.> for solving linear systems, auto-differentiation <cit.> for the derivatives of our potential, and PolyFEM <cit.> as the finite element simulation framework. All experiments are run on a system with an AMD Ryzen Threadripper PRO 3995WX 64-Cores (limited to 16 threads) and 440 of memory. For comparisons against IPC <cit.> and “Convergent IPC” <cit.>, we use the open-source implementations provided in the IPC Toolkit (together with PolyFEM). Please see our supplemental video for animations. All simulation parameters and statistics are summarized in <ref>. Our reference implementation, used to generate all results, will be released as an open-source project.

§.§ Resolving Spurious Forces

In this section, we consider a variety of unit tests where the original IPC formulation introduces spurious forces.

Using larger $\epsilon_\text{trg}$  Suppose the input mesh is relatively fine, either globally (a dense mesh) or locally (adaptively refined e.g. to resolve fine features). Having a large $\epsilon_\text{trg}$ allows one to solve the implicit time-stepping problem faster because the barrier is less numerically stiff, and CCD needs to do less work in the line search if the objects are kept further apart by the potential.  An even stronger reason, when for large deformations, the elements may shrink a lot, is given below

We provide an example of this scenario in <ref>, where it is natural to have small edges around the rounded corners of a square and long edges along the sides. However, doing so restricts the range of usable $\dhat$ for IPC. Using a value ($\dhat=\qty{0.1}{\m}$ in this case) larger than the minimum edge length ($h_{\min}=\qty{0.01}{\m}$) results in spurious forces along the corners and artifacts upon simulation. In contrast, our utilization of an adaptive $\epsilon$ allows us to choose a starting $\epsilon_\text{trg}$ that results in zero initial contact force, while still being large for the non-refined regions.

Spurious stresses: rest configuration 3D Top row: two cubes initially separated by less than $\hat{d}=\qty{0.025}{\m}$. IPC artificially repels the two blocks while ours does not. Bottom row: a block with a slit of width less than $\hat{d}=\qty{0.0125}{\m}$. Our method does not introduce spurious forces across the slit.

Spurious stress in the rest configuration Even if $\dhat < h_{\min}$ is satisfied, the IPC formulation may still have spurious stress in the rest configuration. We show two such scenarios in <ref>.

In the first case (<ref> Top), two blocks in the initial configuration sit on the plane with the initial distance between blocks less than
$\dhat$. With IPC, they incorrectly start sliding apart without external force applied where the blocks should stay still at all times.
In the second case (<ref> Bottom), IPC causes the slit to expand at the top without external forces and stress appears at the bottom.

Our approach avoids spurious forces between close objects by using an adaptive $\epsilon_{trg}$. For the contact between primitives on the same object, due to our constraints on the angle of contact pairs, we can use $\epsilon_{trg}$ larger than the edge lengths without activating contact at rest configuration (with a wide range of $\alpha$).

Spurious stresses: compression An initially valid choice for $\hat{d}$ can lead to spurious contact (arrows and color) when using IPC. Left: cube at rest with a minimum edge length of $\qty{0.03}{\m}$. Center: upon compression with $\hat{d} = \qty{0.0125}{\m}$, IPC introduces spurious contacts. Right: our method avoids this by considering the angle when finding candidate contacts. Bottom: IPC introduces artificial contact pressure.

Spurious stresses: expansion Upper left: an initially deflated balloon with rest thickness of $\qty{0.1}{\m}$ is expanded by an outward force of $\qty{5}{\kilo\N}$. Center: with $\hat{d} = 0.02$, the balloon eventually becomes thin enough that spurious contact forces (represented as red arrows) between its inner and outer layers appear. Lower left: a cross-section showing spurious forces on both the inner and outer layer of the balloon. Right: our method can inflate the balloon until it is arbitrarily thin without introducing artificial contact forces. Bottom: artificial contact pressure is introduced by IPC.

Spurious stress under deformation
Even for a well-chosen $\dhat$ parameter and a benign rest configuration, IPC can still add spurious forces upon deformation.

In <ref>, we compress a cube mesh to 33% of its original height. We assign a Poisson ratio of 0 to the cube to show an example where no bulging and/or folding of the surface occurs. The value for $\dhat$ is initially chosen such that no points are in contact, but upon compression distances shrink and IPC introduces spurious contact forces on the sides of the cube (despite them being flat). Our formulation does not have spurious forces because we use the angle of contact to build our collision candidate set.

In <ref>, we inflate a balloon modeled as a volumetric membrane of thickness $\qty{0.1}{\m}$.[While one could model the balloon using codimensional shell elements to avoid this issue, modeling the thickness may be important for analysis or design. For example, with constant outward pressure, the balloon oscillates in thickness, which a shell model would not capture.] As the balloon inflates, its walls get thinner and eventually they become thinner than the initially chosen $\dhat=\qty{0.02}{\m}.$ At this point, IPC treats the inner and outer sides as in contact, introducing forces between the two sides. Our collision candidate set for one side does not include the other as it accounts for the angle of contact, avoiding this issue.

In <ref>, we extrude the extended $\chi$-shaped structure from <cit.> to 3D and simulate its compression. When the shape buckles it forms cusp contacts at the corners. We show that Convergent IPC <cit.> exhibits large contact forces in these regions, IPC <cit.> exhibits less forces, while our method has much less contact forces than both methods. Our method reduces the contact forces around the cusp because of the constraints on the angle between faces – only when edges are close enough to parallel, i.e. the angle of the cusp becomes small enough, the contact is activated. Without this constraint, large $\kappa$ results in spurious stresses as in Convergent IPC, while small $\kappa$ results in convergence issues in IPC – the number of iterations is much higher due to tiny distance of $\qty{e-11}{\m}$ between contact pairs.

Cusp compression Compression of a $\chi$-shaped structure. For comparison, we use the adaptive barrier stiffness for IPC, $\kappa$ equals Young's modulus both for convergent IPC and our method. In the last time step, IPC took 257 iterations, convergent IPC 58, and ours 72.

§.§ Infinite Potential

Setup for corner hitting a plane discussed in <ref>.

Finite Potential for Zero Distance A square impacts a plane at the corner of the square with a velocity of 234. We refine both meshes (only at the corner of the square of efficiency). Here we plot the height of the squares tip over time. Left: for Convergent IPC <cit.>, we see that the minimum distance shrinks with every refinement of the mesh. This is a consequence of the continuous model having a finite potential. Middle: due to our choice of barrier function, our method exhibits the same trajectory (all plots overlap) under refinement.

Consider a corner hitting a plane, and suppose they are both refined (<ref> Right). In this case, the integral of the Convergent IPC <cit.> potential over the plane will be finite, even if the apex of the corner is directly on the plane (violation of <ref>). Although each descritization will be infinite, the implication is that the total potential will decrease as we refine, so the actual distance will decrease to zero.

<Ref> show the configuration: for a gap of $r$ and a distance from the closest point on the plane $\Delta x$, points along the square will be at a distance of $\Delta x + r$. Assume the potential for points on the corner is just $\ln\left(\frac{\Delta x + r}{\dhat}\right)$,[The extra quadratic term does not affect the conclusions.] then integrating over $x$ we get

\[\left( r+L \right) \ln  \left( {\frac {r+L}{\dhat}} \right) -\ln  \left( {\frac {r}{\dhat}} \right) r-L,\]

and for $r = 0$ (i.e. contact), a finite limit $\ln  \left( {\frac {L}{\dhat}} \right) L-L$.

The force (i.e., the potential derivative) at $r = 0$ is $\sim\ln(r)$ though (i.e. is infinite), so for a finite force the static equilibrium problem will always have a valid solution. This is not the case however for a dynamic problem: if the kinetic energy exceeds the finite potential, then the barrier cannot prevent contact. We show this case in <ref>. One can see the Convergent IPC model results in ever decreasing distances as the mesh is refined (without refining $\dhat$). Our method in comparison exhibits the same trajectory for all levels of refinement. This is because we use a barrier function whose integral over the surface is not finite as the distance goes to zero, and because the potential converges under mesh refinement.

Given that the minimum distance doesn't depend on the mesh resolution, one can pick $\epsilon_\text{trg}$,  independent of the edge lengths, based on the deired accuracy on the contact handling.

§.§ Parameter Ablation and Convergence Study

To study the effects of our simulation parameters we conduct ablation studies for $\alpha$ and $\epsilon_\text{trg}$.

Lower values of $\alpha$ (<ref>) make the solution more accurate, preventing spurious forces from appearing but requiring more iterations. This effect is more pronounced in scenes with cusps: for simpler scenes the difference in iteration count is less noticeable (<ref>).

We plot in <ref> the effect of $\epsilon_\text{trg}$ on the number of iterations. As the support of the contact potential increases the problem becomes softer and the number of iteration decreases accordingly.

We also perform a convergence study for three different scenarios show in <ref>. We see convergence under mesh refinement in all of these scenes. Importantly, we use a fixed $\epsilon_\text{trg}$ (unlike Convergent IPC <cit.> which requires co-refinement of $\dhat$).

We note that when our discretization is applied to PL surface, treating it as a piecewise smooth surface, the integrals for edge potentials and summations for vertex potentials are introduced, with the relative scale of potentials determined by the constant $L$.  As the surface is refined, the scale of the constant also needs to be adjusted, for the refined mesh potential to approximate the smooth surface potential.
We leave a rigorous study of convergence of the discrete potential to the potential of the limit smooth surface as future work.

Parameter-study: $\alpha$ Simulation of <ref> with different $\alpha$ including the total number of solver iterations. As $\alpha$ increases, the nonlinear problem becomes softer, hence fewer iterations. However, large $\alpha$ introduces artifacts and spurious stresses in some cases.

Parameter-study: $\epsilon_\text{trg}$ Simulation of two bunnies colliding with varying $\alpha$ and $\epsilon_\text{trg}$. We show the plot (left) of the number of iterations over $\epsilon_\text{trg}$, with a different $\alpha$ for each curve, as well as the initial frame (middle) and the colliding frame (right) of the simulation. For simple scenes, for a large range of $\alpha$, the number of iterations is similar; as $\epsilon_\text{trg}$ increases, it takes fewer iterations to converge.

Convergence of potential under mesh refinement in fixed scenes with fixed $\epsilon_\text{trg}$. We pick 2 configurations in 2D and 1 in 3D, start with the coarse mesh, and compute the potential under mesh refinement. In 2D, the potentials of (A) and (B) converge with an order of 0.82 and 1.57. Due to the nearly singular cusp, (A) converges slower than (B). In 3D, although we cannot further refine the surface mesh due to memory restriction (the number of contact candidates grows quadratically under refinement), it shows a trend of converging.

§.§ 3D Examples

Lastly, we reproduce challenging 3D simulation examples of <cit.>. First, we validate our method on the 3D unit tests proposed by <cit.> in <ref>. We see similar results to those shown in <cit.>, but we highlight one improved result in <ref> where we see reduced spurious tangential movement compared to <cit.>. Second, we reproduce the dolphin funnel (<ref>), trash-compactor (<ref>), and mat twist (<ref>) examples. Each of these examples features large deformations and complex contacts. Just as in <cit.>, we robustly handle these scenarios and prevent intersections and inversions at every step.

Erleben tests We reproduce the test-cases of <cit.>. Top: initial conditions involving challenging exact point-point, point-edge, and edge-edge collisions. Bottom: as in <cit.>, our approach robustly passes all the tests.

Erleben test: cliff edge We reproduce the cliff edge test-case of <cit.>. Middle: <cit.> passes the test but introduces spurious horizontal forces, causing the top block to rotate. Right: our method significantly reduces extra sliding due to the restriction in normal forces.

Dolphin in a funnel We reproduce the funnel test from <cit.> using our method. Top: an elastic dolphin is pulled through a small tube. Middle: this causes extreme deformations. Bottom: the dolphin squeezes through without artifact and recovers its original shape.

Trash-Compactor We reproduce the trash-compactor example from <cit.> using our method. Left: three objects are placed in a compactor. Middle: the objects are compressed. Right: the compactor releases and the shapes return to their original shape without intersections.

Mat-Twist We reproduce the mat-twist example from <cit.>, demonstrating, overall, our method works similarly to the original IPC. Left: our simulation at 10 after 2 rounds of twisting at both ends. Right: at 27.5 after 5.5 rounds of twisting.

§ CONCLUSION

Incremental potential contact <cit.> enabled a qualitative improvement in robustness of accurate simulation of deformable objects with complex contact. In this paper, we revisit the IPC contact formulation to relate it to a family of potentials defined for a broad class of surfaces, that satisfy a set of natural requirements for barrier potentials for contact.  We demonstrate that applying these principles leads to a new formulation that is similar in efficiency to IPC but alleviates some of its shortcomings.

Future work While we have not considered friction in this paper, the extension to friction is straightforward: the IPC formulation applies essentially with no changes. Our derivation assumes surfaces without boundaries at several steps, which needs modifications for co-dimensional objects; a possible direction for future work is handling surfaces with boundaries and co-dimensional surfaces. Additionally, it may be beneficial to use a high-order quadrature to better integrate the continuous formulation (in addition to including the closest-point quadrature).

Simulation statistics For each simulation we report geometry, (minimum, average, and maximum edge length $h$), time step $\Delta t$, material (Young's modulus $E$, Poisson ratio $\nu$, and density $\rho$), $\epsilon_\text{trg}$, maximum memory used, as well as average timing and number of Newton iterations. A value of 0 for iterations indicates the optimization did not run because the initial configuration was at a force equilibrium (i.e., no spurious forces at rest).

#1     #1#1#1     #1     #1     #1     #1#1          #1#1

titleMitsuba Renderer.

Accessed: 2023-10-10.

[Alappat et al(2020)]
authorpersonChristie Alappat, personAchim Basermann, personAlan R. Bishop, personHolger Fehske, personGeorg Hager, personOlaf Schenk, personJonas Thies, and personGerhard Wellein. year2020.
A Recursive Algebraic Coloring Technique for Hardware-Efficient Symmetric Sparse Matrix-Vector Multiplication.
journalACM Transactions on Parallel Computing volume7, number3, Article articleno19 (dateJune year2020), numpages37 pages.

[Armero and Petőcz(1998)]
authorpersonFrancisco Armero and personEva Petőcz. year1998.
Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems.
journalComputer methods in applied mechanics and engineering volume158, number3-4 (year1998), pages269–300.

[Belgacem et al(1998)]
authorpersonF.B. Belgacem, personP. Hild, and personP. Laborde. year1998.
The mortar finite element method for contact problems.
journalMathematical and Computer Modelling volume28, number4 (year1998), pages263–271.
Recent Advances in Contact Mechanics.

[Belytschko et al(2000)]
authorpersonTed Belytschko, personWing Kam Liu, and personBrian Moran. year2000.
booktitleNonlinear Finite Elements for Continua and Structures. Vol. volume26.
publisherJohn Wiley & Sons, Ltd. 38–3926–38–3926 pages.

[Belytschko and Neal(1991)]
authorpersonTed Belytschko and personMark O Neal. year1991.
Contact-impact by the pinball algorithm with penalty and Lagrangian methods.
journalInternat. J. Numer. Methods Engrg. volume31, number3 (year1991), pages547–572.

[Bollhöfer et al(2019)]
authorpersonMatthias Bollhöfer, personAryan Eftekhari, personSimon Scheidegger, and personOlaf Schenk. year2019.
Large-scale Sparse Inverse Covariance Matrix Estimation.
journalSIAM Journal on Scientific Computing volume41, number1 (year2019), pages380–401.

[Bollhöfer et al(2020)]
authorpersonMatthias Bollhöfer, personOlaf Schenk, personRadim Janalik, personSteve Hamm, and personKiran Gullapalli. year2020.
State-of-the-Art Sparse Direct Solvers.
journalParallel Algorithms in Computational Science and Engineering (year2020), pages3–33.

[Buck and Orloff(1995)]
authorpersonGregory Buck and personJeremey Orloff. year1995.
A simple energy function for knots.
journalTopology and its Applications volume61, number3 (year1995), pages205–214.

[Carpenter et al(1991)]
authorpersonNicholas J Carpenter, personRobert L Taylor, and personMichael G Katona. year1991.
Lagrange constraints for transient finite element surface contact.
journalInternational journal for numerical methods in engineering volume32, number1 (year1991), pages103–128.

[Christensen et al(1998)]
authorpersonPW Christensen, personAnders Klarbring, personJong-Shi Pang, and personNiclas Strömberg. year1998.
Formulation and comparison of algorithms for frictional contact problems.
journalInternat. J. Numer. Methods Engrg. volume42, number1 (year1998), pages145–173.

[Deuflhard et al(2008)]
authorpersonPeter Deuflhard, personRolf Krause, and personSusanne Ertel. year2008.
A contact-stabilized Newmark method for dynamical contact problems.
journalInternat. J. Numer. Methods Engrg. volume73, number9 (year2008), pages1274–1290.

[Do Carmo(2016)]
authorpersonManfredo P Do Carmo. year2016.
booktitleDifferential geometry of curves and surfaces: revised and updated second edition.
publisherCourier Dover Publications.

[Du et al(2023)]
authorpersonYinwei Du, personYue Li, personStelian Coros, and personBernhard Thomaszewski. year2023.
titleNo Free Slide: Spurious Contact Forces in Incremental Potential Contact.
[arxiv]2308.01696 [cs.GR]

[Duong and Sauer(2019)]
authorpersonThang Duong and personRoger Sauer. year2019.
A concise frictional contact formulation based on surface potentials and isogeometric discretization.
journalComputational Mechanics volume64, number4 (year2019), pages951–970.

authorpersonKenny Erleben. year2018.
Methodology for Assessing Mesh-Based Contact Point Methods.
journalACM Transactions on Graphics volume37, number3, Article articleno39 (dateJuly year2018), numpages30 pages.

[Ferguson et al(2020)]
authorpersonZachary Ferguson et al year2020.
titleIPC Toolkit.

[Ferguson et al(2023)]
authorpersonZachary Ferguson, personPranav Jain, personDenis Zorin, personTeseo Schneider, and personDaniele Panozzo. year2023.
High-Order Incremental Potential Contact for Elastodynamic Simulation on Curved Meshes. In booktitleACM SIGGRAPH 2023 Conference Proceedings (Los Angeles, CA, USA) (seriesSIGGRAPH '23). publisherAssociation for Computing Machinery, addressNew York, NY, USA, numpages11 pages.

[Ferguson et al(2021)]
authorpersonZachary Ferguson, personMinchen Li, personTeseo Schneider, personFrancisca Gil-Ureta, personTimothy Langlois, personChenfanfu Jiang, personDenis Zorin, personDanny M. Kaufman, and personDaniele Panozzo. year2021.
Intersection-Free Rigid Body Dynamics.
journalACM Transactions on Graphics (Proceedings of SIGGRAPH) volume40, number4, Article articleno183 (dateJuly year2021), numpages16 pages.

[Guennebaud et al(2010)]
authorpersonGaël Guennebaud, personBenoît Jacob, et al year2010.
titleEigen v3.

[Harmon et al(2011)]
authorpersonDavid Harmon, personDaniele Panozzo, personOlga Sorkine, and personDenis Zorin. year2011.
Interference-Aware Geometric Modeling.
journalACM Transactions on Graphics volume30, number6 (dateDec. year2011), pages1–10.

[Harmon et al(2009)]
authorpersonDavid Harmon, personEtienne Vouga, personBreannan Smith, personRasmus Tamstorf, and personEitan Grinspun. year2009.
Asynchronous contact mechanics.
journalACM Transactions on Graphics (Proceedings of SIGGRAPH) volume28, number3 (year2009), pages1–12.

[Harmon et al(2008)]
authorpersonDavid Harmon, personEtienne Vouga, personRasmus Tamstorf, and personEitan Grinspun. year2008.
Robust Treatment of Simultaneous Collisions.
journalACM Transactions on Graphics (Proceedings of SIGGRAPH) volume27, number3 (dateAug. year2008), pages1–4.

[Hauret and Le Tallec(2006)]
authorpersonPatrice Hauret and personPatrick Le Tallec. year2006.
Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact.
journalComputer methods in applied mechanics and engineering volume195, number37-40 (year2006), pages4890–4916.

[Hiermeier et al(2018)]
authorpersonMichael Hiermeier, personWolfgang A Wall, and personAlexander Popp. year2018.
A truly variationally consistent and symmetric mortar-based contact formulation for finite deformation solid mechanics.
journalComputer Methods in Applied Mechanics and Engineering  volume342 (year2018), pages532–560.

[Hüeber and Wohlmuth(2006)]
authorpersonS. Hüeber and personB.I. Wohlmuth. year2006.
booktitleMortar methods for contact problems.
publisherSpringer Berlin Heidelberg, addressBerlin, Heidelberg, pages39–47.

authorpersonAmin Joodaky. year2020.
titleMechanics and Design of Polymeric Metamaterial Structures for Shock Absorption Applications.
thesistypePh. D. Dissertation. schoolPurdue University Graduate School.

[Kaldor et al(2008)]
authorpersonJonathan M. Kaldor, personDoug L. James, and personSteve Marschner. year2008.
Simulating knitted cloth at the yarn level. In booktitleACM SIGGRAPH 2008 Papers (Los Angeles, California) (seriesSIGGRAPH '08). publisherAssociation for Computing Machinery, addressNew York, NY, USA, Article articleno65, numpages9 pages.

[Kane et al(1999)]
authorpersonCouro Kane, personEduardo A Repetto, personMichael Ortiz, and personJerrold E Marsden. year1999.
Finite element analysis of nonsmooth contact.
journalComputer Methods in Applied Mechanics and Engineering volume180, number1-2 (year1999), pages1–26.

[Kloosterman et al(2001)]
authorpersonG Kloosterman, personRudolf MJ van Damme, personAntonius H van den Boogaard, and personJ Huetink. year2001.
A geometrical-based contact algorithm using a barrier method.
journalInternat. J. Numer. Methods Engrg. volume51, number7 (year2001), pages865–882.

[Krause and Zulian(2016)]
authorpersonRolf. Krause and personPatrick. Zulian. year2016.
A Parallel Approach to the Variational Transfer of Discrete Fields between Arbitrarily Distributed Unstructured Finite Element Meshes.
journalSIAM Journal on Scientific Computing volume38, number3 (year2016).

[Lan et al(2022)]
authorpersonLei Lan, personDanny M. Kaufman, personMinchen Li, personChenfanfu Jiang, and personYin Yang. year2022.
Affine Body Dynamics: Fast, Stable and Intersection-Free Simulation of Stiff Materials.
journalACM Transactions on Graphics (Proceedings of SIGGRAPH) volume41, number4, Article articleno67 (dateJuly year2022), numpages14 pages.

[Lan et al(2021)]
authorpersonLei Lan, personYin Yang, personDanny Kaufman, personJunfeng Yao, personMinchen Li, and personChenfanfu Jiang. year2021.
Medial IPC: Accelerated Incremental Potential Contact with Medial Elastics.
journalACM Transactions on Graphics (Proceedings of SIGGRAPH) volume40, number4, Article articleno158 (dateJuly year2021), numpages16 pages.

[Larionov et al(2021)]
authorpersonEgor Larionov, personYe Fan, and personDinesh K. Pai. year2021.
Frictional Contact on Smooth Elastic Solids.
journalACM Trans. Graph. volume40, number2, Article articleno15 (dateapr year2021), numpages17 pages.

[Laursen and Love(2002)]
authorpersonTA Laursen and personGR Love. year2002.
Improved implicit integrators for transient impact problems—geometric admissibility within the conserving framework.
journalInternat. J. Numer. Methods Engrg. volume53, number2 (year2002), pages245–274.

[Li et al(2020)]
authorpersonMinchen Li, personZachary Ferguson, personTeseo Schneider, personTimothy Langlois, personDenis Zorin, personDaniele Panozzo, personChenfanfu Jiang, and personDanny M. Kaufman. year2020.
Incremental Potential Contact: Intersection- and Inversion-free Large Deformation Dynamics.
journalACM Transactions on Graphics volume39, number4, Article articleno49 (dateJuly year2020), numpages20 pages.

[Li et al(2023)]
authorpersonMinchen Li, personZachary Ferguson, personTeseo Schneider, personTimothy Langlois, personDenis Zorin, personDaniele Panozzo, personChenfanfu Jiang, and personDanny M. Kaufman. year2023.
titleConvergent Incremental Potential Contact.
[arxiv]2307.15908 [math.NA]

[Li et al(2021)]
authorpersonMinchen Li, personDanny M. Kaufman, and personChenfanfu Jiang. year2021.
Codimensional Incremental Potential Contact.
journalACM Transactions on Graphics (Proceedings of SIGGRAPH) volume40, number4, Article articleno170 (year2021).

[Maas et al(2012)]
authorpersonSteve A. Maas, personBenjamin J. Ellis, personGerard A. Ateshian, and personJeffrey A. Weiss. year2012.
FEBio: Finite Elements for Biomechanics.
journalJournal of Biomechanical Engineering volume134, number1 (dateFeb. year2012).

[Otaduy et al(2009)]
authorpersonMiguel Otaduy, personRasmus Tamstorf, personDenis Steinemann, and personMarkus Gross. year2009.
Implicit Contact Handling for Deformable Objects.
journalComputer Graphics Forum  volume28 (dateApril year2009), pages559–568.

[Paoli and Schatzman(2002a)]
authorpersonLaetitia Paoli and personMichelle Schatzman. year2002a.
A numerical scheme for impact problems I: The one-dimensional case.
journalSIAM J. Numer. Anal. volume40, number2 (year2002), pages702–733.

[Paoli and Schatzman(2002b)]
authorpersonLaetitia Paoli and personMichelle Schatzman. year2002b.
A numerical scheme for impact problems II: The multidimensional case.
journalSIAM journal on numerical analysis volume40, number2 (year2002), pages734–768.

[Puso and Laursen(2004)]
authorpersonMichael A. Puso and personTod A. Laursen. year2004.
A mortar segment-to-segment contact method for large deformation solid mechanics.
journalComputer Methods in Applied Mechanics and Engineering volume193, number6 (year2004), pages601–629.

[Sauer and De Lorenzis(2013)]
authorpersonRoger A Sauer and personLaura De Lorenzis. year2013.
A computational contact formulation based on surface potentials.
journalComputer Methods in Applied Mechanics and Engineering  volume253 (year2013), pages369–395.

[Schneider et al(2019)]
authorpersonTeseo Schneider, personJérémie Dumas, personXifeng Gao, personDenis Zorin, and personDaniele Panozzo. year2019.

[Strzelecki and von der Mosel(2013)]
authorpersonPaweł Strzelecki and personHeiko von der Mosel. year2013.
Tangent-Point Repulsive Potentials for a Class of Non-smooth m-dimensional Sets in $\mathbb{R}^n$. Part I: Smoothing and Self-avoidance Effects.
journalJournal of Geometric Analysis volume23, number3 (year2013), pages1085–1139.

[Taylor and Papadopoulos(1993)]
authorpersonRobert L Taylor and personPanayiotis Papadopoulos. year1993.
On a finite element method for dynamic contact/impact problems.
journalInternational journal for numerical methods in engineering volume36, number12 (year1993), pages2123–2140.

[Verschoor and Jalba(2019)]
authorpersonMickeal Verschoor and personAndrei C Jalba. year2019.
Efficient and accurate collision response for elastically deformable models.
journalACM Transactions on Graphics volume38, number2, Article articleno17 (dateMarch year2019), numpages20 pages.

[Vola et al(1998)]
authorpersonDidier Vola, personElaine Pratt, personMichel Jean, and personMichel Raous. year1998.
Consistent time discretization for a dynamical frictional contact problem and complementarity techniques.
journalRevue Européenne des éléments finis volume7, number1-3 (year1998), pages149–162.

[Vouga et al(2011)]
authorpersonEtienne Vouga, personDavid Harmon, personRasmus Tamstorf, and personEitan Grinspun. year2011.
Asynchronous Variational Contact Mechanics.
journalComputer Methods in Applied Mechanics and Engineering volume200, number25 (year2011), pages2181–2194.

authorpersonPeter Wriggers. year1995.
Finite Element Algorithms for Contact Problems.
journalArchives of Computational Methods in Engineering  volume2 (dateDec. year1995), pages1–49.

[Wriggers and Laursen(2006)]
authorpersonPeter Wriggers and personTod A Laursen. year2006.
booktitleComputational contact mechanics. Vol. volume2.

[Xie et al(2023)]
authorpersonTianyi Xie, personMinchen Li, personYin Yang, and personChenfanfu Jiang. year2023.
A Contact Proxy Splitting Method for Lagrangian Solid-Fluid Coupling.
journalACM Transactions on Graphics (Proceedings of SIGGRAPH) volume42, number4 (year2023), pages1–14.

[Youett et al(2019)]
authorpersonJonathan Youett, personOliver Sander, and personRalf Kornhuber. year2019.
A globally convergent filter-trust-region method for large deformation contact problems.
journalSIAM Journal on Scientific Computing volume41, number1 (year2019), pagesB114–B138.

[Yu et al(2021a)]
authorpersonChris Yu, personCaleb Brakensiek, personHenrik Schumacher, and personKeenan Crane. year2021a.
Repulsive Surfaces.
journalACM Trans. Graph. volume40, number6 (year2021).

[Yu et al(2021b)]
authorpersonChris Yu, personHenrik Schumacher, and personKeenan Crane. year2021b.
Repulsive Curves.
journalACM Trans. Graph. volume40, number2 (year2021).

§ VERIFICATION OF PROPERTIES

The properties and are accomplished by construction.

is also attained by construction, as all functions used in the definition of the potential are $C^1$.  One subtle point is that the distance functions used in the discrete version of the potential are also $C^1$, except the edge-to-edge distance that require mollification. However, the closest point on an element to a given point $f(x)$ is not a $C^1$ function of $f(x)$, and for this reason we introduce mollification whenever these points are used.

Our arguments are informal and we leave a complete proof with all the necessary technical assumptions as future work.

Generality In the case of smooth surfaces, the potential is clearly well-defined as it involves integrating a function bounded from above, as there is a uniform bound on the distance between a point $x$ and a point in its contact set $y$.

For piecewise smooth surfaces, including PL surfaces, the situation is more subtle, as points may be arbitrarily close to their contact set, if points are located on faces sharing an edge (Figure <ref>). However, the contact area for such points becomes increasingly small, so if the point-point potential does not grow very quickly, the integrated potential is well-defined.

The contact set for two point $x_1$ and $x_2$ on a concave corner may be arbitrarily close, but the contact area for such points becomes increasingly

We assume that the angle $\gamma$ between the normals of two faces sharing an edge, is never zero.
Then for a point at distance $r$ from the edge the distance to the contact set is  $d \sim r \sin \gamma$,

and the area of the contact set, in 3D, $\sim r^2\sin\alpha$. Thus, if potential power does not exceed $2$, then the integral of the potential over the contact set for a point at distance $r$ to the edge, remains bounded, as $d^{-n} r^2$ is bounded. Similarly, in $2D$, the potential needs to grow no faster than $1/r$.

Barrier The interplay between the discretization and the continuum formulation is quite subtle in this case.

In the continuum case, clearly, it is not sufficient to choose a point-to-point potential function that becomes infinite for coinciding points: the integral of such function may be finite, which would mean that the barrier is not preventing the contact.

Consider first the integral of the potential over one surface $S_y$ with respect to the variable $y$, with $x$ fixed. Further, assume an isolated contact point with the distance between surfaces growing as $r^n$ with the distance to the contact point. For a smooth contact $n=2$, for a point contact (e.g. with a cone), $n = 1$.

Consider a model problem with two surfaces in planar contact on a disk of unit radius, at distance $h$, and assume $\hat{d}$ and $\alpha$ to be large. Then the total contact potential is approximated by

\[
2\pi \int_{r=0}^R\int_{q=0}^R\int_{\theta=0}^{2\pi} ( (r\cos\theta  - q)^2 + r^2 \sin^2\theta + h^2)^{p} d\theta dr dq
\]

where $R$ is the contact zone radius,  can be computed explicitly, and  for the point-to-point potential with $p=2$,  can be shown to be  $4\pi^2 ( \ln(2) -\ln(h))^2$, i.e., proportional to $\ln(h)^2$.  If we further assume that one of the surfaces has a cone point, i.e., the distance between the surface and the plane increases linearly w.r.t. the distance along the plane to the contact point, with $h$ replaced with $(h+r)$ above, then the rate of growth reduces to $\ln(h)$, but still results in a repulsive force.

For a potential with a lower rate of growth, the limit interaction potential can be finite, hence it does not result in a barrier in the continuum case.

We note however, that for our construction of the discrete potential this does not necessarily cause a problem for any fixed resolution: as the quadrature for the potential is defined using the minimal distance between contact pairs, for any barrier potential (e.g., $\ln(d)$ used in <cit.>)
one still has a barrier property for any discrete approximation but the strength of the barrier decreases with refinement, and the minimal distance between objects decreases.