PBI: Position-Based Dynamics Handles Updated Lagrangian Inelasticity

Abstract.

Position-based Dynamics (PBD) and its extension, eXtended Position-based Dynamics (XPBD), have been predominantly applied to compliant constrained elastodynamics, with their potential in finite strain (visco-) elastoplasticity remaining underexplored. XPBD is often perceived to stand in contrast to other meshless methods, such as the Material Point Method (MPM). MPM is based on discretizing the weak form of governing partial differential equations within a continuum domain, coupled with a hybrid Lagrangian-Eulerian method for tracking deformation gradients. In contrast, XPBD formulates specific constraints, whether hard or compliant, to positional degrees of freedom. We revisit this perception by investigating the potential of XPBD in handling inelastic materials that are described with classical continuum mechanics-based yield surfaces and elastoplastic flow rules. Our inspiration is that a robust estimation of the velocity gradient is a sufficiently useful key to effectively tracking deformation gradients in XPBD simulations. By further incorporating implicit inelastic constitutive relationships, we introduce a plasticity in-the-loop updated Lagrangian augmentation to XPBD. This enhancement enables the simulation of elastoplastic, viscoplastic, and granular substances following their standard constitutive laws. We demonstrate the effectiveness of our method through high-resolution and real-time simulations of diverse materials such as snow, sand, and plasticine, and its integration with standard XPBD simulations of cloth and water.

Figure 1. PBI supports simulating a wide range of classical continuum elastoplastic material models such as Von-Mises plasticine, Drucker-Prager sand, and Cam Clay snow, as well as their interactions with traditional PBD materials such as Position-based Fluid.

1. Introduction

Position-based Dynamics (PBD) [Müller et al., 2007] and its extension eXtended Position-based Dynamics (XPBD) [Macklin et al., 2016] are widely adopted in compliant constrained dynamics, particularly favored for their performance and simplicity for graphics applications such as rigid bodies [Müller et al., 2020], soft bodies [Bender et al., 2014], cloth [Müller et al., 2007], rods [Umetani et al., 2015] and hair [Müller et al., 2012]. When simulating mesh-based elasticity, it is straightforward to model XPBD constraints with FEM hyperelastic energies defined over explicit mesh topology. This allows for effective simulations of elasticity while maintaining the stability and speed of XPBD. For example, quardratic energy potentials [Chen et al., 2023] can be formulated using XPBD-style constraints. Macklin and Muller [2021] and Ton-That et al. [2023] reformulated stable Neo-Hookean [Smith et al., 2018] to demonstrate XPBD’s capability in handling nonlinear elasticity.

For inelasticity, on the other hand, two significant challenges emerge. Firstly, topology changes during material splitting and merging introduce great complexity in maintaining a high quality mesh, often necessitating remeshing. Secondly, while there have been explorations that enhances Position-based Fluids (PBF) [Macklin and Müller, 2013] with the conformation tensor [Barreiro et al., 2017] for viscoelastic fluids, it remains underexplored for XPBD to model physically-grounded finite strain (visco-) elastoplastic constitutive laws from classical continuum mechanics, such as von-Mises [Mises, 1913], Drucker-Prager [Drucker and Prager, 1952] and Herschel-Bulkley [Herschel and Bulkley, 1926] flow rules. Being able to simulate them would greatly improve XPBD’s versatility and intuitive controllability of material parameters.

Figure 2. Noodles. We simulate noodles modeled using Von Mises plasticity as it is pressed through a cylindrical mold.

In comparison with XPBD’s success in mesh-based materials, Material Point Method (MPM)’s development in graphics over the past decade has majorly focused on inelastic phenomena with topology change. MPM is based on the weak form of governing PDEs and employs a hybrid Lagrangian-Eulerian approach for spatial discretization. It handles topology changes, self collision, and finite strain deformation without overhead, and thus has been used for many continuum inelastic phenomena such as snow [Stomakhin et al., 2013], lava [Stomakhin et al., 2014], sand [Klár et al., 2016], mud [Tampubolon et al., 2017], metal [Wang et al., 2020], foam [Ram et al., 2015], and fracture [Wolper et al., 2019]. While substantial progress has been made in these areas, MPM exhibits several notable drawbacks, including excessive numerical dissipation due to particle-grid transfers [Jiang et al., 2015], artificial stickiness that hampers material separation [Fei et al., 2021], and resolution-dependent gaps between colliding materials [Jiang et al., 2017]. None of these issues are present in XPBD. This raises a natural question: Can we simulate MPM-style phenomena using XPBD instead?

Table 1. PBI keeps PBD’s pure particle nature of the Degrees of Freedom while allowing MPM-style plasticity and granular material modeling through an updated Lagrangian treatment of the deformation and classical continuum mechanics-based elastoplastic flow rules.

	Lagrangian	DOF	Plasticity	Granular
MPM	Updated	Grid	Flow Rule	Continuum
PBD	Total	Particle	N/A	Sphere Approx.
Ours	Total & Updated	Particle	Flow Rule	Continuum

Towards addressing this question, we make an important observation: the primary factor that facilitates the modeling of inelasticity lies not in the hybrid Lagrangian-Eulerian nature of MPM, but rather in its use of an updated Lagrangian formulation for the deformation gradient tensor. In particular, one considers the time $n+1$ velocity $\mathbf{v}^{n+1}$ to be defined over the previous time $n$ domain $\Omega^{n}$ through the Lagrangian velocity $\mathbf{V}$ of a particle traced back using the inverse deformation map $\bm{\phi}^{-1}(\mathbf{x},t)$: $\mathbf{v}^{n+1}(\mathbf{x})=\mathbf{V}(\bm{\phi}^{-1}(\mathbf{x},t^{n}),t^{n+1})$ [Jiang et al., 2016]. This enables the derivation of the rate form of the deformation gradient $\mathbf{F}$ given by $\dot{\mathbf{F}}=(\nabla\mathbf{v})\mathbf{F}$, which can be further discretized into $\mathbf{F}^{n+1}=(\mathbf{I}+\Delta t\nabla\mathbf{v}^{n+1})\mathbf{F}^{n}$, allowing one to track the deformation gradient without referring to a material space configuration.

Inspired by this observation, if we can track the deformation gradient tensor $\mathbf{F}$ using an updated Lagrangian view in XPBD, then by treating $\mathbf{F}$ as a function of XPBD degrees of freedom, we can modify XPBD to resemble a “material point” approach. As detailed in later sections, this task reduces to robustly computing and differentiating the velocity gradient tensor. We present developments surrounding these ideas. By further incorporating an implicit plasticity treatment and additional stability-enhancing components, our method leverages the efficiency and simplicity of PBD while capturing the complex inelastic material responses typically associated with MPM; see Table. 1. In summary, our contributions include:

• An updated Lagrangian augmentation for XPBD that tracks meshless deformation gradients and per-particle constraints;

• PBI, a fully implicit plasticity-aware, velocity-based augmented XPBD algorithm capable of handling continuum mechanics-based elastoplastic/viscoplastic laws;

• An investigation for practical stability enhancements, such as XSPH and position correction, and validations of our method with various practical examples.

Figure 3. PBI fits into traditional XPBD pipeline and naturally couples updated Lagrangian materials (viscoplastic paint) and mesh-based geometry (cloth).

Figure 4. PBI simulates Hershel-Bulkley shear thinning ($h=0.3$), viscoplastic ($h=1.0$), and shear thickening materials ($h=3.0$), where $h$ controls a power law flow rate detailed in [Yue et al., 2015].

2. Related Work

Inelasticity with PBD

Müller et al. [2007] introduced PBD, which replaces internal forces with positional constraints and produces appealing, stable and real-time simulations. Its first-order convergence was studied by Plunder and Merino-Aceituno [2023]. XPBD [Macklin et al., 2016], an extension of PBD, utilizes the compliant-constraint framework [Tournier et al., 2015] to uniformly handle soft and hard constraints to simulate elasticity. Our work follows the latest XPBD paradigm [Macklin et al., 2019] with substeps. Other PBD materials include rigid body [Müller et al., 2020], soft body [Bender et al., 2014], cloth [Müller et al., 2007], hair [Müller et al., 2012], elastic rod [Umetani et al., 2015], sand [Macklin et al., 2014], fluid [Macklin and Müller, 2013] with surface tension [Xing et al., 2022] and their unified couplings [Macklin et al., 2014; Frâncu and Moldoveanu, 2017; Abu Rumman et al., 2020; Robinson-Mosher et al., 2008]. We refer to Bender et al. [2017] for a comprehensive survey.

For continuum materials, Bender et al. [2014] defined a constraint for the elastic strain energy. Müller et al. [2015] constrained the strain tensor directly instead. Macklin and Muller [2021] reformulated stable Neo-Hookean using XPBD. Plastic deformation and fracture can be modeled by shape matching [Chentanez et al., 2016; Jones et al., 2016; Falkenstein et al., 2017], prioritizing efficiency over accuracy. Macklin et al. [2014] simulated sand as colliding spheres with friction. SideFX Houdini’s Vellum PBD solver further added spring-like cohesion for snow. Without further utilizing continuum mechanics-based inelastic models, these approaches have limited mechanical intuition and physical parameter controllability. A step forward was proposed by Barreiro et al. [2017], which enhances PBF with the conformation tensor for viscoelastic fluids. Nevertheless, it does not incorporate finite strain continuum mechanics, limiting its suitability in modeling general elastoplastic and viscoplastic laws.

Inelasticity with MPM

MPM was introduced by Sulsky et al. [1994] as a hybrid Lagrangian/Eulerian approach for solids. Since its adoption in graphics [Hegemann et al., 2013; Stomakhin et al., 2013], MPM has gained significant attentions for its automatic topology change and material versatility. Snow plasticity was first done by projecting principal stretches [Stomakhin et al., 2013]. Phase change was modeled by [Stomakhin et al., 2014] through a dilational/deviatoric splitting of the constitutive model. Yue et al. [2015] adopted Herschel-Bulkley viscoplasticity for foam. Ram et al. [2015] used the Oldroyd-B model for viscoelastic fluids. Fei et al. [2019] developed an analytical plastic flow approach for shear-dependent liquids. Following Drucker-Prager yield criterion, Klár et al. [2016]; Daviet and Bertails-Descoubes [2016] modeled sand as continuum granular materials and Tampubolon et al. [2017] further added wetting. Wolper et al. [2019, 2020] captured dynamic fracture using Non-Associated Cam-Clay (NACC) plasticity and damage mechanics. Advocating implicit integrators, stiff plastic materials like metal was simulated with Newton-Krylov MPM [Wang et al., 2020], while Fang et al. [2019] used Alternating Direction Methods of Multipliers (ADMM) for viscoelasticity and elastoplasticity and Li et al. [2022] proposed a variational implicit inelasticity formulation.

Inelasticity with Other Discretizations

Smoothed Particle Hydrodynamics (SPH) was originally for incompressible flow. Clavet et al. [2005] added dynamic-length springs for viscoelasticity. Jones et al. [2014] and Müller et al. [2004] solved Moving Least Squares (MLS) for elastoplasticity. Gerszewski et al. [2009] first used deformation gradient tensor with multiplicative elastoplastic decomposition in SPH. Alduán and Otaduy [2011] and Yang et al. [2017] modeled granular materials based on Drucker Prager yielding. Gissler et al. [2020] developed an implicit SPH snow solver similarly to Stomakhin et al. [2013]’s MPM treatment. Using power diagram-based particle-in-cell [Qu et al., 2022] and MLS-MPM [Hu et al., 2018], power plastics [Qu et al., 2023] simulated inelastic flow with an XPBD-style Gauss-Seidel solver. Peridynamics [Silling, 2000] defines pairwise forces and integrates particle interactions. He et al. [2017] combined peridynamics with projective dynamics [Bouaziz et al., 2023] and modeled Drucker-Prager plasticity. Chen et al. [2018] used isotropic linear elasticity with plasticity and simulated fracture.

3. Method

We start with briefly reviewing XPBD [Macklin et al., 2016], which lets hyperelasticity be governed by Newton’s equations of motion through a potential $U(\bm{x})$: $\bm{M}\ddot{\bm{x}}=-\nabla U^{T}(\bm{x})$, where $\bm{M}$ is the mass matrix and $\bm{x}=[x_{1},x_{2},...,x_{n}]^{T}$ is the unknown position states. XPBD assumes $U(\bm{x})$ can be further expressed as $U(\bm{x})=\frac{1}{2}\mathbf{C}(\bm{x})^{T}\bm{a}^{-1}\mathbf{C}(\bm{x})$, where $\bm{C}=[C_{1}(\bm{x}),C_{2}(\bm{x}),...,C_{m}(\bm{x})]^{T}$ contains $m$ constraints, and $\bm{a}$ is a block diagonal compliance matrix. The elastic internal force $\bm{f}$ and Lagrange multiplier $\bm{\lambda}$ are then shown to be

(1) $$\bm{f} =-\nabla\mathbf{C}(\bm{x})^{T}\bm{a}^{-1}\mathbf{C}(\bm{x}),$$

(2) $$\bm{\lambda} =-\tilde{\bm{a}}^{-1}\mathbf{C}(\bm{x}).$$

where $\tilde{\bm{a}}=\frac{\bm{a}}{h^{2}}$ and $h=\Delta t$ is the time step size.

Figure 5. Candy Camponotus. We simulate the brittle fracture of a candy shaped like a camponotus falling onto the ground.

3.1. Rewriting StVK Elasticity as Constraints

For modeling elasticity, we adopt the St. Venant-Kirchhoff (StVK) model with Hencky strains. The elastoplastic behavior of isotropic materials is characterized in the principal stretch space $\bm{\Sigma}$ via singular value decomposition (SVD) [Stomakhin et al., 2012] of the deformation gradient $\bm{F}$. For StVK we have energy density $\Psi$

(3) $$\Psi=\mu\text{tr}\left(\log\left(\bm{\Sigma}\right)^{2}\right)+\frac{\lambda}{2}\left(\text{tr}\left(\log\left(\bm{\Sigma}\right)\right)\right)^{2}=\Phi/V_{0},$$

where $V_{0}$ is an element’s rest volume and $\Phi$ is the element’s total potential energy, assuming piecewise constant element deformations.

To convert $\Phi$ into constraints, our first option is to separately handle each term. As done by Macklin and Muller [2021] for Neo-Hookean, by utilizing $\Phi=\frac{1}{2}\frac{1}{\alpha}C^{2}$, we could define two constraints for the $\mu$ and $\lambda$ term respectively as

(4) $$\alpha_{\mu}=1/(2\mu V_{0}), \text{ and}\ C_{\mu}=\sqrt{\text{tr}\left(\log\left(\bm{\Sigma}\right)^{2}\right)};$$

(5) $$\alpha_{\lambda}=1/(\lambda V_{0}), \text{ and}\ C_{\lambda}=\text{tr}\left(\log\left(\bm{\Sigma}\right)\right).$$

Alternatively, as done by Qu et al. [2023] on power diagrams, we can absorb Lamé parameters from $\alpha$ to $C$ and define a single $\bm{F}$-dependent constraint through $\Phi=V_{0}\Psi(\bm{F})=\frac{1}{2\alpha}C(\bm{F})^{2}$ to achieve

(6) $$\alpha={1}/{V_{0}},\text{ and}\ C(\bm{F})=\sqrt{2\Psi(\bm{F})},$$

which clearly halves the number of constraints by allowing more nonlinearity in each $C(\bm{x})$. In practice we find both options effective, and adopt the single-constraint version for efficiency.

3.2. Gradient is All You Need

Unlike mesh-based representations [Macklin and Muller, 2021] which use linear FEM to model the deformation map $\bm{x}=\bm{\phi}(\bm{X},t)$ and compute deformation gradient $\bm{F}=\partial\bm{\phi}(\bm{X},t)/\partial{\bm{X}}$ using an undeformed reference state $\bm{X}\in\Omega^{0}$, meshless inelastic materials cannot utilize simplex elements due to extreme deformation. We follow Jiang et al. [2016] to take an updated Lagrangian view for deformation rate:

(7) $$\frac{\partial}{\partial t}\bm{F}(\bm{X},t) =\frac{\partial}{\partial t}\frac{\partial\bm{\phi}}{\partial\bm{X}}(\bm{X},t)=\frac{\partial\bm{v}}{\partial\bm{x}}(\bm{\phi}(\bm{X},t),t){\frac{\partial\bm{\phi}}{\partial\bm{X}}(\bm{X},t)},$$

with $\bm{V}(\bm{X},t)=\partial\bm{\phi}(\bm{X},t)/\partial t$ being the Lagrangian velocity whose Eulerian counterpart is $\bm{v}(\bm{x},t)=\bm{V}(\bm{\phi}^{-1}(\bm{x},t),t)$. With time discretization from $t^{n}$ to $t^{n+1}$ and the assumption that the velocity $\bm{v}$ at time $t^{n+1}$ being $\bm{v}^{n+1}(x)$ for $\bm{x}\in\Omega^{n}$, we have

(8) $$\frac{\partial}{\partial t}\bm{F}(\bm{X},t^{n+1})=\frac{\partial\bm{v}^{n+1}}{\partial\bm{x}}(\bm{\phi}(\bm{X},t^{n}))\bm{F}(\bm{X},t^{n}).$$

Taking $\frac{\partial}{\partial t}\bm{F}_{p}(\bm{X}_{p},t^{n+1})\approx({\bm{F}^{n+1}_{p}-\bm{F}^{n}_{p}})/{h}$ for a particle $\bm{X}_{p}$ we get

(9) $$\bm{F}_{p}^{n+1}=\bm{F}_{p}^{n}+h\frac{\partial\bm{v}^{n+1}}{\partial\bm{x}}(\bm{x}_{p}^{n})\bm{F}_{p}^{n}=\left(\bm{I}+h\frac{\partial\bm{v}^{n+1}}{\partial\bm{x}}(\bm{x}_{p}^{n})\right)\bm{F}_{p}^{n}$$

as the evolution of $\bm{F}_{p}^{n+1}$ given $\bm{v}^{n+1}$ and $\bm{F}_{p}^{n}$. With updated Lagrangian, the reference space is thus always $\Omega^{n}$ and there is no need to store $\Omega^{0}$; see Fig. 6. Therefore, to express constraints as functions of positions $C(\bm{F}(\bm{v}^{n+1}(\bm{x}^{n+1})))$, robustly estimating $\partial\bm{v}/\partial\bm{x}$ ($\forall x\in\Omega^{n}$) and its derivative is all one needs.

Figure 6. Deformation Gradient Evolution. The dotted line (bottom) illustrates the evolution of the deformation gradient in the updated Lagrangian view, transitioning from $\bm{F}_{p}^{0}$ (initial configuration) to $\bm{F}_{p}^{n+1}$ (updated configuration), with $\bm{F}_{p}^{n}$ (current configuration) serving as the reference state. This facilitates tracking large deformations. The solid line loop (top) depicts iterations of our PBI algorithm to simulate $t^{n}\rightarrow t^{n+1}$, alternating between an XPBD iteration and a fixed point iteration. During iteration $k$, $\bm{F}_{p}^{E,tr}$ is first estimated based on the current gradient of $\bm{v}^{(k)}$ ($\S$ 3.2). Then, plasticity is applied through projection to obtain $\bm{F}_{p}^{(k+1)}$ ($\S$ 3.3), and finally, $\bm{v}^{(k+1)}$ is updated by solving constraints ($\S$ 4.1).

Accurately estimating meshless velocity gradients is generally challenging. Most meshless shape functions require a dense neighborhood to fulfill the kernel’s normalization condition. Known as neighborhood deficiency, significant accuracy degradation would occur especially for first-order derivatives (such as velocity gradient) in sparse regions. Kernel gradient correction [Bonet and Lok, 1999] and the reproducing kernel particle method [Liu et al., 1995] are examples of strategies in mitigating this problem. Here we adopt reweighting-based kernel gradient correction [Bonet and Lok, 1999]

(10) $$\bm{L}_{p}^{-1}=\sum_{b}V_{b}^{n}\nabla W_{b}(\bm{x}_{p})\otimes(\bm{x}_{b}-\bm{x}_{p}),$$

with SVD based pseudo inverse and Wendland kernels for $W$ and $\nabla W$, where $V_{b}^{n}=V_{b}^{0}\text{det}(\bm{F}_{p}^{n})$ is the time $n$ volume of particle $b$. As in many SPH methods, the correction is indispensable for maintaining simulation stability. See Westhofen et al. [2023] for discussions of similarly viable gradient estimation choices. Our discrete velocity gradient at particle $p$ in $\Omega^{n}$ is then

(11) $$\frac{\partial\bm{v}}{\partial\bm{x}}(\bm{x}_{p}^{n}) =\sum_{b\neq p}V_{b}^{n}(\bm{v}_{b}-\bm{v}_{p})\left(\bm{L}_{p}\nabla W_{b}(\bm{x}_{p}^{n})\right)^{T}.$$

Combining the gradient estimation with Eq. 9 and differentiating per-particle constraint $C_{p}(\bm{F}_{p})$ (Eq. 6) reveals

(12) $$\nabla_{\bm{x}_{b}}{C_{p}}|_{b\neq p} =V_{b}^{n}\frac{\partial C_{p}}{\partial\bm{F}_{p}}{\bm{F}_{p}^{n}}^{T}(\bm{L}_{p}\nabla W_{b}(\bm{x}_{p}^{n})),$$

(13) $$\nabla_{\bm{x}_{p}}C_{p} =-\sum_{b\neq p}\nabla_{\bm{x}_{b}}C_{p},$$

which provides us all necessary constraint derivatives in XPBD.

Figure 7. Dam Breach. Our method can be seamlessly coupled with PBF [Macklin and Müller, 2013] to simulate sand and water mixture.

Figure 8. Hourglass. Sand in hourglass accumulates at the bottom. The material is modeled with Drucker-Prager plasticity [Klár et al., 2016].

Figure 9. Comparison on notched sand block fall. Inital (left), our fully implicit treatment (middle), and semi-implicit plasticity (right).

3.3. Implicit Plasticity

Plasticity in continuum mechanics is typically solved with return mapping, denoted as $\mathcal{Z}(\cdot)$, which adjusts strains according to a plastic flow rule. It projects an elastic predictor $\bm{F}^{E,tr}$ onto the yield surface to ensure an inequality constraint on the stress.

To make plasticity implicit, we propose to alternate between (1) an XPBD iteration with a projected stress and (2) a stress projection. This is essentially a fixed point iteration similarly to Li et al. [2022]:

(14) $$\bm{F}_{p}^{(k+1)}\leftarrow\mathcal{Z}\left(\bm{F}_{p}^{E,tr}\left(\bm{v}^{(k)}\left(\bm{F}_{p}^{(k)}\right)\right)\right),$$

where $\bm{v}^{(k)}(\bm{F}_{p}^{(k)})$ is the updated velocity by previous $k$ XPBD iterations based on $\bm{F}_{p}^{(k)}$ in the previous iteration. In contrast to Li et al. [2022]’s fixed point iteration on $\mathcal{Z}$ which functions as an independent outer loop of a full Newton optimization, our design establishes a fixed point on $\bm{F}$, updating variables directly impacted by the fixed point iteration within XPBD iterations; see Fig. 6 top and Alg. 1. Resultingly, our implicit plasticity treatment introduces negligible extra cost on top of what implicit elasticity already necessitated.

Nonetheless, convergence in fixed-point iterations depends on an initial guess sufficiently close to the solution, among other conditions of the implicit function. In this paper, we do not monitor quantitative plasticity convergence since few XPBD iterations are needed for visually plausible results. We do emphasize the importance of implicit plasticity and compare it with a semi-implicit treatment which only applies plasticity at the end of a time step; see $\S$ 5.1.

4. Algorithm

Here we detail the PBI pipeline and its seamlessly integration into existing XPBD. Our pseudocode for advancing a time step using velocity-based XPBD is summarized in Alg. 1.

4.1. Algorithm Overview

Similarly to MPM, we use material particles to discretize the continuum. Each particle $p$ is governed by a constitutive model-induced constraint $C_{p}$ (Eq. 6) and a plastic return mapping operator $\mathcal{Z}$. We use $C_{p}$ to denote particle-wise inelasticity constraints ($|\{p\}|=N=\#\text{ particles}$) and $C_{i}$ to denote traditional PBD constraints ($|\{i\}|=M$).

Due to our dependency on velocity gradients, it is more natural to reparametrize XPBD with velocities rather than positions as primary unknown variables. Closely resembling position-based XPBD, we solve for velocities $\bm{v}^{n+1}$ and Lagrange multipliers $\bm{\lambda}^{n+1}$ that satisfies

(15) $$\bm{M}(\bm{v}^{n+1}-\tilde{\bm{v}}^{n})-\nabla\bm{C}(\bm{v}^{n+1})^{T}\bm{\lambda}^{n+1} =\bm{0},$$

(16) $$\bm{C}(\bm{v}^{n+1})+\tilde{\bm{a}}\bm{\lambda}^{n+1} =\bm{0},$$

by updating per-particle inelastic constraint $C_{p}$’s corresponding Lagrangian multiplier $\lambda_{p}$ with

(17) $$\Delta\lambda_{p}=\frac{-C_{p}-\tilde{\alpha}_{p}\lambda_{p}}{\sum_{b=1}^{N}\frac{1}{m_{b}}|\nabla_{\bm{x}_{b}}C_{p}(\bm{x})|^{2}+\tilde{\alpha_{p}}},$$

where $\alpha_{p}=1/V_{p}^{0}$ and $\tilde{\alpha}_{p}=\alpha_{p}/h^{2}$. The velocity update is given by:

(18) $$\Delta\bm{v}=\left(\bm{M}^{-1}\nabla\mathbf{C}(\bm{x})^{T}\Delta\bm{\lambda}\right)/h.$$

The velocities and multipliers are jointly updated by a colored Gauss-Seidel iteration (see $\S$ 4.3 and Alg. 1):

(19) $$\lambda_{p}^{(k+1)}\leftarrow\lambda_{p}^{(k)}+\Delta\lambda_{p},\quad\bm{v}^{(k+1)}\leftarrow\bm{v}^{(k)}+\Delta\bm{v}.$$

Algorithm 1 Simulating $t^{n}\rightarrow t^{n+1}$ with PBI

Neighbor search using $\bm{x}^{n}$ $\triangleright$ $\S$ 4.2

Evaluate kernel gradient correction $\bm{L}_{p}$ $\triangleright$ $\S$ 3.2

$\bm{v}\leftarrow\bm{v}^{n}+h\bm{M}^{-1}\bm{f}_{\text{ext}}$

$\bm{\lambda}=\bm{0}$

for number of XPBD iterations do

     for all $p\in\{C_{p}\}$ looping with colored Gauss-Seidel do $\triangleright$ $\S$ 4.3

         $(\nabla\bm{v})_{p}\leftarrow$ evaluateVelocityGradient$(\bm{x}_{p}^{n},\bm{v}_{p})$

         $\bm{F}_{p}\leftarrow(I+h(\nabla\bm{v})_{p})\bm{F}_{p}^{n}$

         $\bm{F}_{p}\leftarrow\mathcal{Z}(\bm{F}_{p})$ $\triangleright$ $\S$ 3.3

         if $C_{p}(\bm{F}_{p})\neq 0$ then

              $\Delta\lambda_{p}=\frac{-C_{p}-\tilde{\alpha}\lambda}{\sum_{i=1}^{N}\frac{1}{m_{i}}|\nabla_{\bm{x}_{i}}C_{p}(\bm{x})|^{2}+\tilde{\alpha}}$

              $\lambda_{p}\leftarrow\lambda_{p}+\Delta\lambda_{p}$

              $\Delta\bm{v}=\frac{1}{h}\bm{M}^{-1}\nabla C_{p}(\bm{x})^{T}\Delta\lambda_{p}$

              $\bm{v}\leftarrow\bm{v}+\Delta\bm{v}$

         end if

     end for

     for all $i\in\{C_{i}\}$ looping with colored Gauss-Seidel do

         $\Delta\bm{v}=\frac{1}{h}\bm{M}^{-1}\nabla C_{i}(\bm{x})^{T}\Delta\lambda_{i}$ (e.g., position correction)$\triangleright$ $\S$ 4.4

         $\bm{v}\leftarrow\bm{v}+\Delta\bm{v}$

     end for

end for

$\bm{v}^{n+1}\leftarrow\bm{v}$

Perform XSPH smoothing of $\bm{v}^{n+1}(\bm{x}^{n})$ $\triangleright$ $\S$ 4.4

Update $\bm{F}_{p}^{n+1}$ and apply constitutive models $\triangleright$ $\S$ 4.5

$\bm{x}^{n+1}\leftarrow\bm{x}^{n}+h\bm{v}^{n+1}$

Note: $\alpha=1/V_{p}^{0}$ and $\tilde{\alpha}=\alpha/h^{2}$ for each constraint.

4.2. Particle Neighbor Search

We reconstruct the neighbor information for each particle at the beginning of each timestep, similarly to Macklin and Müller [2013]. A comprehensive overview of CPU- and GPU-based neighborhood search methods is surveyed by Ihmsen et al. [2014]. Each material particle is assigned the same kernel radius in our discretization scheme. We adopt a uniform spatial-grid-based method for neighborhood searches following Hoetzlein [2014]. Particles are spatially stored in cells while neighbor lists $\mathcal{N}_{p}=\{b\mid\|x_{p}-x_{b}\|_{2}\leq 2r\}$ are determined by querying adjacent grid cells using the Wendland kernel’s support radius of $2r$.

4.3. Colored Gauss-Seidel

Original PBD and XPBD frameworks solve constraints iteratively using nonlinear Gauss-Seidel [Müller et al., 2007; Macklin et al., 2014, 2016]. In contrast, Macklin and Müller [2013] adopted a Jacobi-style iteration for fluids, solving each constraint independently to enhance parallelism. We found that for high resolution and often high stiffness simulations considered in this paper, Jacobi iterations too slowly propagate information and often suffer from non-convergence (also noted by Macklin et al. [2014]). Thus we implement colored Gauss-Seidel to maximize convergence, parallelism, and GPU throughput. We assign particles into cells with $\Delta x=2r$ ($\S$ 4.2). $3^{d}$ colors are specified for eliminating dependencies between constraints in a $d$-dimensional simulation. We process all cells of the same color in parallel while constraints $C_{p}\in c_{i}$ corresponding to all particles in the same cell $c_{i}$ are computed serially.

Table 2. Parameters and Statistics. We summarize the parameters and timing statistics, including maximum particle numbers, the Wendland kernel radius, the average time per frame, the XPBD iterations per substep, and the time step size $\Delta t$ for various demos described in $\S$ 5.2. Material-related parameters are detailed in the last two columns. The abbreviations are NACC for Non-Associated Cam-Clay [Wolper et al., 2019], DP for Drucker-Prager [Klár et al., 2016], VM for Von Mises [Li et al., 2022], and HB for Herschel-Bulkley [Yue et al., 2015]. In addition to the basic settings of the material (density $\rho$, Youngs Modulus $E$, and Poisson Ratio $\nu$), we include model-specific parameters arranged as follows: 1) NACC: $(\rho,E,\nu,\alpha_{0},\beta,\xi,M)$; 2) DP: $(\rho,E,\nu,\phi_{f},c_{0})$; 3) VM: $(\rho,E,\nu,\sigma_{\gamma})$; and 4) HB: $(\rho,E,\nu,\sigma_{\gamma},h,\eta)$. See references for detailed explanations of these parameters.

demo	particle #	radius	ave sec/frame	iter #	$\Delta t_{\text{frame}}$	$\Delta t_{\text{step}}$	material	material parameters
(Fig. 13) Hitman	1.05M	$1/512$	41.2	10	1/100	$4\times 10^{-5}$	NACC	$(4,2\times 10^{4},0.3,-0.005,0.05,30,1.85)$
(Fig. 13) Snow Dive	2.48M	$1/512$	84.7	5	1/100	$4\times 10^{-5}$	NACC	$(4,1\times 10^{4},0.3,-0.0005,0.05,30,1.85)$
(Fig. 2) Noodles	1.18M	$1/256$	67.1	10	1/40	$1\times 10^{-4}$	VM	$(1,2\times 10^{4},0.3,76.9)$
(Fig. 8) Hourglass	1.01M	$1/1024$	47.9	5	1/24	$1\times 10^{-4}$	DP	$(1,3.537\times 10^{5},0.3,35,0)$
(Fig. 7) Dam Breach	4.00M	$1/384$	198.5	7	1/24	$2.5\times 10^{-4}$	DP	$(1,400,0.4,30,0.0007)$
(Fig. 5) Camponotus	1.12M	$1/1024$	52.0	10	1/100	$4\times 10^{-5}$	NACC	$(2,2\times 10^{4},0.35,-0.02,0.5,1,2.36)$
(Fig. 3) Cloth	1.10M	$1/512$	43.7	10	1/100	$5\times 10^{-5}$	HB	$(100,14754,0.475,50,1,10)$
(Fig. 16) Wrist	20K	$1/256$	0.015	5	1/100	$2\times 10^{-4}$	HB	$(100,2250,0.125,10,1,10)$

The efficiency of this implementation heavily depends on the average Particle Per Cell (PPC), as particles within the same cell are traversed sequentially. To optimize PPC and ensure an even particle distribution, we utilize Poisson disk sampling [Bridson, 2007] during the initial particle placement.

Figure 10. Hydraulic test on a stiff aluminum wheel: initial (left), w/ (middle) and w/o (right) position correction. The simulation suffers from artificial fracture and instability without our correction.

Figure 11. PBI effectively captures viscoplastic coiling; as the shear modulus $\mu$ increases from left to right, coils become more elastic.

4.4. XSPH and Position Correction

The goal of eXtended Smoothed Particle Hydrodynamics (XSPH) is to incorporate artificial viscosity for mitigating nonphysical oscillations in dynamics observed in SPH-based simulations, which is more observable when the material is stiff and the constraints become hard to solve. We follow Schechter and Bridson [2012]’s simpler XSPH-style noise damping after velocity update by blending in surrounding particle velocities in $\Omega^{n}$:

(20) $$\bm{v}^{n+1,new}_{p}\leftarrow\bm{v}_{p}^{n+1}+c\sum_{b}V_{b}^{n}(\bm{v}_{b}^{n+1}-\bm{v}_{p}^{n+1})\cdot W_{b}(\bm{x}_{p}^{n}).$$

XSPH encourages smooth and coherent motion for inelastic materials in this paper. We use $c=0.01$ in all examples.

In particle-based simulations, including but not limited to SPH, FLIP [Brackbill and Ruppel, 1986; Zhu and Bridson, 2005], and APIC [Jiang et al., 2015], non-physical particle clumping and uneven particle distributions are common issues due to accumulation of advection errors. While MPM is insensitive to the problem due to its structural grid nature, it is crucial for pure particle-based methods like SPH and our approach to maintain reasonably even particle distributions. Without doing so can strongly impair simulation quality and convergence, particularly when material stiffness is high.

Various techniques in graphics have been proposed to address this by shifting particle positions [Ando et al., 2012; Ando and Tsuruno, 2011; Kugelstadt et al., 2019]. However, we point out that while these methods work well for fluids, they are problematic for updated Lagrangian simulations because such positional shifting is transparent to the evolution of deformation gradients, leading to discrepancy between the positions of points and their perceived deformations. Fortunately within XPBD we can directly adopt a point-point distance constraint [Macklin et al., 2014]

(21) $$C(\bm{x}_{p},\bm{x}_{b})=\left\|\bm{x}_{p}-\bm{x}_{b}\right\|_{2}-r+\epsilon\geq 0$$

inside nonlinear iterations. Here $r$ represents the particle kernel radius and $\epsilon$ is a small gap threshold used to determine when corrective action is needed for nearby particle pairs. We set $\epsilon$ to $0.25r$.

4.5. Deformation Gradient Update

Strain variables within XPBD iterations are temporary. After arriving at the post-XSPH velocity $\bm{v}^{n+1}$, we update both the deformation gradient state and position state using the same velocity – an essential subtlety to maintain their consistency:

(22) $${\bm{F}}_{p}^{n+1}=\mathcal{Z}\left(\left(\bm{I}+h\frac{\partial\bm{v}^{n+1}}{\partial\bm{x}}(\bm{x}_{p}^{n})\right)\bm{F}_{p}^{n}\right),$$

where inelastic return mapping is also applied to ensure the stored elastic deformation gradient is within the yield region.

Figure 12. A typical breakdown of the total computational cost of our framework. We take the Hourglass example (Fig. 8) for demonstration.

Figure 13. Hitman and Snow Dive. We successfully reproduce realistic and complex snow behaviors, such as a snowball hitting a person (top) and a person falling into a snow ground (bottom), using NACC [Wolper et al., 2019] constitutive model.

5. Results

Here we evaluate and benchmark our Position-based Inelasticity (PBI) framework in terms of visual results against traditional XPBD and MPM methods, as detailed in $\S$ 5.1. Additionally, we present various demonstrations in $\S$ 5.2 that illustrate PBI’s effective handling of diverse phenomena. We use Intel Core i9-14900KF CPU with 32GB memory and NVIDIA GeForce RTX 4090. We model common inelastic materials including Cam-Clay (NACC) [Wolper et al., 2019] snow and fracture, Drucker-Prager [Klár et al., 2016; Tampubolon et al., 2017] sand, Von Mises [Li et al., 2022] plasticine and metal, and Herschel-Bulkley [Yue et al., 2015] foam.

5.1. Evaluation

Intuitive Parameters

We simulate sand collapsing with varying friction angles $\phi_{f}$. Our method reproduces characteristic piling shapes. Similarly we can easily control the fluidity of viscoplastic goo (Fig. 11). A smaller $\mu$ gives a more fluid-like appearance, while a larger $\mu$ leads to more elastic behavior. In Fig. 4, we compare viscoplastic, shear thinning, and shear thickening materials by only altering the Herschel-Bulkley power parameter $h$. Upon impact with a ground plane, the shear thickening material exhibits low flow rates under high stress, behaving elastically and bouncing off. Conversely, the shear thinning material flows immediately due to its higher flow rate.

Figure 14. Vanilla XPBD v.s. MPM v.s. PBI. We simulate two sand blocks collide (top) and sand column collapse (bottom) with vanilla XPBD [Macklin et al., 2014, 2016] (left), MPM [Klár et al., 2016] (middle), and our method (right). Vanilla XPBD fails to accurately replicate the correct friction angle. While both MPM and our method successfully model the continuum behavior of sand, our method results in a more uniform particle distribution. In the rightmost image, we quantitatively compare these results by plotting the average distance between each particle and its nearest neighbor for each simulation frame, with the shaded areas indicating the respective standard deviations.

Comparisons to Vanilla XPBD and MPM

We compare PBI sand with both vanilla XPBD [Macklin et al., 2014, 2016], which employs a point-wise friction model, and explicit MPM with Drucker-Prager plasticity [Klár et al., 2016]; see Fig. 14. This comparison includes simulations of two sand blocks collide (top) and sand column collapse (bottom). All methods apply a timestep of $\Delta t=0.1$ ms, with identical initial sampling positions for all particles across the methods. The vanilla XPBD (left) approach fails to accurately replicate the correct friction angle upon sand settling. While both MPM (middle) and PBI (right) successfully model the continuum behavior of sand, our method achieves a more uniform particle distribution, avoiding the sparsity, clumping, and artificial grid $\Delta x$-gap phenomena typically caused by MPM solvers. For a quantitative analysis, we also plotted the average distance between each particle and its nearest neighbor per frame, displayed on the far right of Fig. 14. Each particle was initially positioned at intervals equal to the particle radius. The average distance relative to the initial state in the MPM decreases rapidly post-collision, accompanied by an increase in the variance of the distance, whereas our method maintains both the relative distance and variance stably throughout the simulation, indicating a more consistent particle distribution. It is also noteworthy that our method aligns more closely with the approach of [Yue et al., 2018] compared to MPM, opting for Discrete Element Method (DEM) to capture more discrete behaviors near the free surface.

Implicit Plasticity

We evaluate our fully implicit plasticity treatment by replacing it with a semi-implicit method, which only performs return mappings at the end of the time steps as in Stomakhin et al. [2013] $\S$ 4.5), while XPBD iterations only address elasticity. As shown in Fig. 9, the semi-implicit approach (right) can lead to severe artifacts. This occurs because the forces generated by stresses outside the yield surface cause the continuum to behave more like a purely elastic body. This artifact results from the semi-implicit method’s failure to account for plasticity during the XPBD solve, leading to an overestimation of the material’s resistance to tensile deformation. In contrast, our method (middle) fully incorporates plasticity in the XPBD iterations and avoids such artifacts.

Figure 15. Cantilever beams modeled with StVK constitutive model using our method with different stiffness (top). Relative residual errors in $\bm{h(\bm{v},\bm{\lambda})}=\mathbf{C}(\bm{v})+\tilde{\bm{a}}\bm{\lambda}$ for a single frame taken from the above examples (bottom).

Figure 16. Real-time Vision Pro^TM Interaction. Interactive manipulation of a viscoelastic fluid, consisting of up to 20K particles, is simulated at 30fps. Hand movements are tracked using an Apple Vision Pro^TM device and VisionProTeleop [Park and Agrawal, 2024] open-sourced protocol, with these motions applied to a virtual hand interacting with the coiling fluid.

Figure 17. We simulate viscoplastic monsters falling to the ground using varying numbers of particles. We plot the computation time for each frame for 8K, 56K, 400K, and 3M (left) particles, demonstrating consistent behaviors across different particle counts. We also plot the average cost per particle (right). The running overhead of our algorithm decreases significantly as the number of particles increases, showing strong superlinear scalability.

Convergence

We study the convergence of our method using cantilever beams of varying stiffness, $E=10^{4},10^{5},10^{6}$, respectively (Fig. 15). We set the density at $100$ and the timestep at $\Delta t=5$ ms. We also compare the relative residual errors of the Gauss-Seidel and Jacobi solvers in our method. PBI with Gauss-Seidel can converge stably with a large timestep. In contrast, the Jacobi solver only converges with softer materials and struggles with high stiffness. This is consistent with observations about XPBD in prior work. Given that the materials discussed in our paper are predominantly very stiff, we opted for grid-colored Gauss-Seidel as our solver. Although using a large timestep is feasible with sufficient iterations, it becomes cost-inefficient if too many iterations are required. Thus, following Macklin et al. [2019], we employ a small timestep.

Position Correction

To validate the importance of position correction, we conducted a hydraulic test on a highly stiff aluminum wheel, as shown in Fig. 10. Without position correction, areas with significant deformation and stress suffer from gradient estimation errors due to uneven particle distribution, resulting in artificial fractures and eventually simulation instability. With position correction, however, we can reliably simulate high-stiffness materials under extensive deformation. This example also shows our capability in animating metal ductility using Von Mises plasticity.

Scalability

To demonstrate the scalability of our method, we simulate viscoplastic monsters hitting the ground using 8K, 56K, 400K, and 3M particles, respectively. We maintain a constant simulation time step of $\Delta t=0.1$ ms, perform 10 iterations of XPBD per substep, and apply consistent material parameters across all simulations. Our method consistently replicates material behavior at varying resolutions, as depicted in Fig. 17. We measure and plot the computation time for each individual frame (left) alongside the average computation time per particle (right). The average computation times per particle for 8K, 56K, 400K, and 3M are 0.75 ms, 0.178 ms, 0.075 ms, and 0.038 ms, respectively. These results highlight strong scalability; as the number of particles increases, the colored Gauss-Seidel solver can more effectively exploit GPU resources, significantly reducing the overall computational overhead per particle.

Timing Breakdown

Fig. 12 illustrates the GPU computational cost breakdown for the Hourglass simulation example. In the breakdown, Collision Detection refers to the time spent on constructing the LBVH [Karras, 2012] and querying for collision between point-triangle pairs. Neighbor Search covers the time taken to build the background grid and neighborhood list by querying adjacent cells (see $\S$ 4.2 for details). Solve Inelasticity Constraints involves our colored Gauss-Seidel solver for inelasticity constraints, including fixed-point implicit plasticity treatment. Solve Other Constraints accounts for the time spent on resolving all other constraints, such as point-triangle distance constraints for boundary conditions, position correction, as well as stretching, bending, and density constraints in other examples. XSPH and Update Particle State are detailed in $\S$ 4.4 and $\S$ 4.5, respectively. The majority of our framework’s computational time is spent on solving inelasticity constraints, while additional stability enhancements like XSPH and position correction contribute a relatively small overhead.

5.2. Demos

Complex materials with up to millions of particles, such as snow (Fig. 13), sand (Fig. 8), mud (Fig. 2, Fig. 7), viscoplastic paint (and its coupling with traditional XPBD cloth) (Fig. 3), and brittle fracture (Fig. 5) can be simulated with PBI. The timing and parameters are summarized in Table. 2.

Real-time Interaction

We further showcase our method in interactive applications where very small time steps are impractical. The convergence and stability of our approach enable interactive performance in moderately complex scenarios involving 20K particles. For this application, we employ the Jacobi solver due to its significant parallelism capabilities on GPUs for small-scale simulation. In this example, we use an Apple Vision Pro VR device and VisionProTeleop [Park and Agrawal, 2024] to track hand motions and enable a virtual hand to interact with viscous fluids.

6. Discussion

In summary, PBI is a novel updated Lagrangian enhancement for XPBD, enhancing its capability for simulating complex inelastic behaviors governed by continuum mechanics-based constitutive laws. Further incorporating an implicit plasticity treatment and stability enhancements, PBI can be easily integrated into standard XPBD to open up its new simulation possibilities.

While PBI supports a broad range of material behaviors, Maxwell viscoelastic materials [Fang et al., 2019] necessitate more specialized treatment. Also, interactions between sand and water mixtures occur primarily at the material boundary. Simulating actual porous media [Tampubolon et al., 2017] and fluid sediment mixture [Gao et al., 2018] with proper momentum transfer are interesting future work.

Another direction is to optimize parallelism for solving inelastic per-particle constraints. Although grid-colored Gauss-Seidel significantly improves performance over sequential iterations in large-scale simulations, it underperforms in small-scale cases on modern GPUs due to low utilization, obstructing many interactive rate experiments. A tailored real-time solver would be interesting.