Frictional Contact Solving for Material Point Method

Accurately handling contact with friction remains a core bottleneck for Material Point Method (MPM), from reliable contact point detection to enforcing frictional contact laws (non-penetration, Coulomb friction, and maximum dissipation principle). In this paper, we introduce a frictional-contact pipeline for implicit MPM that is both precise and robust. During the collision detection phase, contact points are localized with particle-centric geometric primitives; during the contact resolution phase, we cast frictional contact as a Nonlinear Complementarity Problem (NCP) over contact impulses and solve it with an Alternating Direction Method of Multipliers (ADMM) scheme. Crucially, the formulation reuses the same implicit MPM linearization, yielding efficiency and numerical stability. The method integrates seamlessly into the implicit MPM loop and is agnostic to modeling choices, including material laws, interpolation functions, and transfer schemes. We evaluate it across seven representative scenes that span elastic and elasto-plastic responses, simple and complex deformable geometries, and a wide range of contact conditions. Overall, the proposed method enables accurate contact localization, reliable frictional handling, and broad generality, making it a practical solution for MPM-based simulations in robotics and related domains.

💡 Research Summary

This paper addresses one of the most persistent challenges in Material Point Method (MPM) simulations: robust and accurate frictional contact handling. The authors propose a complete pipeline that integrates contact detection and contact resolution directly into the implicit MPM time‑stepping loop, thereby enabling large time steps, reusing existing linearizations, and remaining agnostic to material models, interpolation kernels, and particle‑to‑grid transfer schemes.

Contact detection is performed at a sub‑cell level using particle‑centric geometric primitives. Each particle carries a deformed tetrahedron defined by its deformation gradient. By querying tetrahedron‑tetrahedron proximity with a variant of the GJK algorithm, the method identifies candidate contacts and computes for each a contact point, a unit normal, and an orthonormal contact frame (two tangential directions and the normal). This particle‑centric approach eliminates the diffusion inherent in grid‑node based detection and guarantees that the detected geometry matches the actual deformed shape of the material.

Contact resolution casts the problem as a Nonlinear Complementarity Problem (NCP) over contact impulses λ. The NCP encodes the Signorini non‑penetration condition, the Coulomb friction law, and the maximum dissipation principle. Mathematically, the complementarity reads
 K_μ ∋ λ ⟂ σ + Γ_μ(σ) ∈ K*_μ,
with σ = G λ + g, where G = H A⁻¹ Hᵀ is the Delassus operator, H is the contact Jacobian that maps grid velocities to contact‑frame velocities, A⁻¹ is the admittance matrix obtained from the implicit MPM linearization, and g is the free contact velocity when no impulses are applied.

To solve the NCP, the authors adopt an Alternating Direction Method of Multipliers (ADMM) scheme, a strategy previously successful in rigid‑body and FEM friction solvers. Each ADMM iteration consists of a linear step that solves a global system involving (A + ρ Hᵀ H) (the same matrix factorization used for the implicit MPM step) and a proximal projection step that projects each contact impulse onto the friction cone K_μ. Convergence is monitored via primal and dual residuals, and the algorithm typically converges within a few hundred iterations, with most time steps achieving residuals below 10⁻⁶.

A key implementation detail is the barycentric lift used to deposit contact information onto the grid. The contact point inside a deformed tetrahedron is expressed as a weighted combination of the tetrahedron’s four vertices. Four “proxy” particles are created at these vertices, each carrying the corresponding barycentric weight. These proxies participate in the standard particle‑to‑grid (P2G) transfer, ensuring that contact forces are applied only to grid nodes that already contain mass, thereby avoiding empty‑cell artifacts and preserving locality.

The method is naturally compatible with multi‑object simulations. Each object retains its own background grid; the global contact Jacobian is assembled by concatenating per‑object blocks, while the global admittance matrix becomes block‑diagonal. The ADMM solve then operates on the combined system, preserving the independence of internal dynamics while enforcing global contact constraints.

Experimental validation spans seven representative scenarios: (1) a cube falling onto a plane, (2) multi‑object interaction, (3) a cube sliding down an inclined plane, (4) a plastic cube sliding between vertical walls, (5) a cube entering a fluid‑like medium, (6) a soft robotic finger rubbing a plane, and (7) a sphere impacting a soft quadruped robot. For each scenario the authors test multiple interpolation kernels (linear, quadratic, CK‑MPM), transfer schemes (basic, APIC, MLS), and material models (linear elasticity, St. Venant‑Kirchhoff, co‑rotational, Neo‑Hookean). Results show accurate contact localization, correct sticking‑to‑sliding transitions, and stable convergence across all combinations. Table I reports the average number of active contacts, mean ADMM iterations per step, and residual statistics, demonstrating that more sophisticated kernels (quadratic, CK‑MPM) reduce the number of contacts and iterations, while the method consistently achieves residuals well below the prescribed tolerance.

Performance analysis indicates that the additional cost of the contact solve is modest: the admittance matrix is factorized once per time step per object, and the ADMM iterations reuse this factorization, leading to an overhead of roughly 0.5–1.2 ms per step on a CPU implementation. This is substantially lower than explicit or penalty‑based approaches that require much smaller time steps for stability.

In conclusion, the paper delivers a general, efficient, and mathematically rigorous solution for frictional contact in implicit MPM. By leveraging particle‑centric detection, a global NCP formulation, and ADMM with reuse of the implicit MPM linearization, the method achieves high accuracy without sacrificing computational efficiency. The authors also outline future directions, including GPU acceleration, handling of multiple simultaneous contacts per particle, incorporation of hysteresis or rate‑dependent friction models, and integration with differentiable simulation pipelines for robotics learning and control. This work therefore represents a significant step toward making MPM a practical tool for high‑fidelity, friction‑aware simulations in robotics, graphics, and computational mechanics.