Geometric Contact Potential¶

In addition to the original implementation of Li et al. [2020], we also implement the Geometric Contact Potential (GCP) from Huang et al. [2025].

Motivation¶

A representative frictionless contact problem can be viewed as a geometrically constrained optimization:

where \(u\) is a deformation and \(g\) measures the distance to contact. Barrier-based solvers replace this with an unconstrained problem:

where \(b(u)\) is a barrier potential that increases to infinity as the configuration approaches contact.

All existing contact barrier potentials are defined by aggregating repulsion terms between pairs of points or elements. The key question is how to choose the strength of repulsion so that it is high when two points are close to contact, vanishes when they are far, and satisfies a number of desirable properties. In particular, one important difficulty is distinguishing true contact — arising from different objects or parts of the same surface moving toward each other — from points that remain close simply because they are neighbors on the undeformed rest shape.

The IPC barrier potential [Li et al., 2020] uses a purely distance-based criterion, which can produce spurious repulsive forces between nearby surface points that are not actually in contact. This is especially problematic when the barrier activation distance \(\hat{d}\) is relatively large compared to mesh features. GCP addresses this by introducing geometric criteria — the local minimum constraint and the exterior direction constraint — that identify true contact candidates, eliminating spurious forces and decoupling \(\hat{d}\) from the discretization.

Requirements¶

GCP is designed to satisfy the following requirements simultaneously. No prior potential in the literature satisfies all of them:

1. Finiteness: The total contact potential integral is finite for any piecewise-smooth surface not in contact.

2. Barrier: The potential grows to infinity as the distance to contact approaches zero, guaranteeing that all configurations remain contact-free (when combined with CCD).

3. No spurious forces: In the undeformed (rest) configuration, both the potential and its gradient are zero, so no artificial forces arise.

4. Localization: The potential has a locality parameter \(\hat{d}\) and vanishes if the distance to contact exceeds \(\hat{d}\).

5. Differentiability: The potential depends differentiably (and piecewise twice differentiably) on mesh vertex positions, enabling second-order solvers.

6. Discretization: The discrete version of the potential satisfies all of the above requirements exactly, not just in the limit of refinement.

Formulation¶

The GCP pointwise potential takes the general form:

where \(f\) is the deformation map, \(p_{\epsilon}\) is a barrier function that diverges as its argument approaches zero and vanishes for distances exceeding \(\epsilon(x)\), and \(\gamma(x, y)\) is a directional factor that restricts the potential to a geometrically meaningful interaction set \(C(x, f)\).

The total potential is obtained by integrating (or in the discrete case, summing) over all pairs of surface points.

The formulation involves four main components:

1. Interaction sets \(C(x, f)\): the set of points that are close to being in true contact with \(x\).

2. Adaptive locality \(\epsilon(x)\): a per-point barrier extent.

3. Barrier function \(p_\epsilon\): the distance-based barrier.

4. Directional factor \(\gamma(x, y)\): a smooth factor supported on \(C(x, f)\).

Barrier Function¶

The barrier function is:

where \(r = n - 1\) (\(n \in \{2, 3\}\) is the dimension of the scene), and \(h_\epsilon(z) := \tfrac{3}{2} B^3(2z/\epsilon)\) is a cubic \(C^2\) spline with support \(|z| \leq \epsilon\). Here \(B^3\) is the cubic B-spline basis function:

Because \(h_\epsilon\) vanishes for \(z \geq \epsilon\) and \(z^{-r}\) diverges as \(z \to 0\), the barrier satisfies both the localization and barrier requirements.

Interaction Sets¶

The key insight of GCP is defining interaction sets based on a geometric characterization of contact: points that are close to being local minima of the distance function and for which the displacement vector points toward the surface exterior.

For smooth surfaces, two points \(x\) and \(y\) on the boundary are in contact if (a) they coincide in space, \(f(x) = f(y)\), and (b) their normals have opposite orientation, \(n(x) = -n(y)\). The interaction set captures points close to satisfying both conditions, via two constraints:

Local minimum constraint. A point \(y\) is near a local minimum of the distance from \(f(x)\) if:

where \((\cdot)_+\) denotes normalization to a unit vector and \(0 < \alpha < 1\). This condition checks whether the distance direction is nearly aligned with the surface normal at \(y\) — equivalently, whether \(y\) is near a closest point on the surface to \(f(x)\).

Exterior direction constraint. The displacement vector \(f(y) - f(x)\) should point toward the surface exterior at \(y\) (i.e., \(f(x)\) is on the outward-normal side of \(f(y)\)):

This eliminates interactions between opposite sides of a thin volumetric shell.

The interaction set for smooth surfaces is then:

Note that both constraints are applied symmetrically: \(y\) must be near a local minimum as seen from \(x\), and \(x\) must be near a local minimum as seen from \(y\).

Piecewise Smooth Surfaces¶

Smooth surfaces are typically approximated by piecewise-linear meshes in practice. On piecewise-smooth surfaces (including triangle meshes), normals are not uniquely defined at edges and vertices, so the constraints must be generalized.

Local minimum constraint (generalized). For a point \(y\) on a piecewise-smooth surface, let \(I\) be the set of indices of patches containing \(y\). For each patch \(\partial\Omega_i\), define tangent vectors \(t^k_i(y)\) for \(k = 1, 2, 3\) (two edge tangents and their angle bisector, see figure below). The local minimum condition becomes:

where \(\alpha_t\) is the tangent-direction smoothing parameter.

Exterior direction constraint (generalized). At an edge point, the constraint is evaluated by projecting onto the plane perpendicular to the edge and applying a 2D inside-cone test. At a vertex, where multiple face normals \(n_i(y)\) are available, a conservative smooth formulation is used. The per-face exterior constraint is:

and the mollified constraint ensures that \(\Phi^e_i(x, y) \geq 0\) for at least one incident face \(i\).

Contact between different element types. On piecewise-smooth surfaces, contact can occur between any pair of element types: Face–Face, Face–Edge, Face–Vertex, Edge–Edge, Edge–Vertex, and Vertex–Vertex. Since edge and vertex contact sets have zero measure, directly integrating the potential over the surface would yield zero. GCP addresses this by treating low-dimensional elements separately with an element-measure weighting parameter \(L\):

where \(g, h \in \{\text{Face}, \text{Edge}, \text{Vertex}\}\), and the integrals are area, line, or point evaluations depending on the element dimension.

Directional Factor \(\gamma\)¶

To keep the potential differentiable while restricting it to the interaction set, GCP uses smooth mollification. The directional factor is built from two components:

Local minimum factor (using the smoothed Heaviside function \(H^\alpha\)):

This equals 1 when all tangent-direction constraints are satisfied, and vanishes smoothly as any constraint is violated.

Exterior direction factor (for vertices with multiple normals):

which equals 1 when at least one face normal satisfies the exterior constraint, and vanishes when none do. For edge points, a simpler form using the 2D cone test is used.

The combined factor is:

The symmetry ensures that both points in a pair must satisfy the geometric criteria.

Smooth Heaviside Function¶

The smooth Heaviside function \(H(z; \alpha, \beta) \in C^1(\mathbb{R})\) used throughout the construction satisfies:

The parameters \(\alpha\) and \(\beta\) control the transition region. Since the inputs to \(H\) are dot and cross products of unit vectors, we have \(-1 \leq z \leq 1\), and the basic requirement is:

Adaptive Barrier Localization¶

To satisfy the no spurious forces requirement, the per-point barrier extent \(\epsilon(x)\) is set adaptively:

where \(d_c(x, f_0)\) is the distance to the interaction set in the rest configuration and \(\hat{d}\) is the global maximum activation distance. This ensures that the potential is exactly zero in the rest pose, eliminating spurious forces by construction.

Discretization¶

For the barrier property to hold exactly in the discrete setting, GCP uses closest points on element pairs as quadrature points (rather than standard quadrature). This may reduce integration precision, but guarantees the barrier property at the discrete level. The discrete potential is:

where \((x_i, y_j)\) is a pair of closest points on elements \(G_i\) and \(H_j\), and \(A(\cdot)\) denotes the element measure (area for faces, length for edges, 1 for vertices).

Distance Mollification¶

Since closest points may not depend smoothly on vertex positions (e.g., when a closest point moves from the interior of a face to its boundary), GCP introduces mollification factors \(M(x, y)\) that smoothly vanish near these non-smooth transitions.

For Edge–Vertex contact with \(Q\) on edge \(AB\):

which vanishes as the closest point approaches either endpoint of the edge.

For Face–Vertex contact with face \(ABC\):

which vanishes as the closest point approaches the boundary of the triangle.

For Edge–Edge contact between edges \(AB\) and \(CD\):

Here \(h_c(s) = h((s-1)/c)\) with \(h(z) = z(2-z)\) for \(0 \leq z < 1\) and \(h(z) = 1\) for \(z \geq 1\). These mollifications are \(C^1\) smooth and ensure that when one contact type vanishes due to mollification, another contact type becomes active, preserving the barrier property.

Usage¶

GCP is implemented as separate collision and potential classes. A basic example of computing the GCP potential is as follows.

#include <ipc/ipc.hpp>

// Set up parameters
double dhat = 1e-3;     // barrier activation distance
double alpha_t = 0.5;   // local minimum constraint support
double beta_t = 0.0;    // local minimum constraint offset
double alpha_n = 0.1;   // exterior direction constraint support
double beta_n = 0.0;    // exterior direction constraint offset
int r = 2;              // barrier exponent (dimension - 1)

ipc::SmoothContactParameters params(dhat, alpha_t, beta_t, alpha_n, beta_n, r);

// Build collision set
bool use_adaptive_dhat = true;
ipc::SmoothCollisions collisions;
if (use_adaptive_dhat)
    collisions.compute_adaptive_dhat(collision_mesh, vertices, params);
collisions.build(collision_mesh, vertices, params, use_adaptive_dhat);

// Compute potential
ipc::SmoothContactPotential barrier_potential(params);
double b = barrier_potential(collisions, collision_mesh, vertices);

// Compute gradient
Eigen::VectorXd grad = barrier_potential.gradient(collisions, collision_mesh, vertices);

// Compute Hessian
Eigen::SparseMatrix<double> hess = barrier_potential.hessian(
    collisions, collision_mesh, vertices);

import ipctk

# Set up parameters
dhat = 1e-3       # barrier activation distance
alpha_t = 0.5     # local minimum constraint support
beta_t = 0.0      # local minimum constraint offset
alpha_n = 0.1     # exterior direction constraint support
beta_n = 0.0      # exterior direction constraint offset
r = 2             # barrier exponent (dimension - 1)

params = ipctk.SmoothContactParameters(dhat, alpha_t, beta_t, alpha_n, beta_n, r)

# Build collision set
use_adaptive_dhat = True
collisions = ipctk.SmoothCollisions()
if use_adaptive_dhat:
    collisions.compute_adaptive_dhat(collision_mesh, vertices, params)
collisions.build(collision_mesh, vertices, params, use_adaptive_dhat)

# Compute potential
barrier_potential = ipctk.SmoothContactPotential(params)
b = barrier_potential(collisions, collision_mesh, vertices)

# Compute gradient
grad = barrier_potential.gradient(collisions, collision_mesh, vertices)

# Compute Hessian
hess = barrier_potential.hessian(collisions, collision_mesh, vertices)

Parameter Choices¶

The SmoothContactParameters structure contains the following parameters:

The \(\alpha\) and \(\beta\) parameters for each constraint must satisfy:

As \(\alpha + \beta\) decreases, the support of the Heaviside function shrinks and the function becomes less smooth, making the potential harder to optimize.

Additional internal parameters that may affect behavior:

• Adaptive dhat ratio (default 0.5): Controls the ratio \(\epsilon(x) / d_c(x, f_0)\) in the adaptive barrier localization. Can be set via SmoothContactParameters::set_adaptive_dhat_ratio().

• Element measure \(L\): For vertices, this is set to the average edge length around the vertex; for edges, to the edge length; for faces, \(L\) is not needed. This determines the strength of the potential for low-dimensional contact (edge–edge, edge–vertex, vertex–vertex).

Friction¶

We implement the original friction formulation of Li et al. [2020] following the same style as GCP. Details are covered in the paper Huang et al. [2025].

Comparison with IPC¶

The following summarizes the key differences between GCP and the original IPC barrier: