Theory, History & Comparisons

This document covers what MPM is, where it came from, how the algorithm works step by step, what it’s good for, where it breaks down, and how it stacks up against the neighboring simulation methods.

1. What MPM is

The Material Point Method (MPM) is a hybrid Eulerian–Lagrangian continuum simulation method. Material is represented by Lagrangian material points (particles) that carry mass, velocity, and deformation state, moving freely through space; each step, their data is transferred to a fixed Eulerian Cartesian background grid, where forces are computed and momentum is updated, and the result is transferred back to the particles. The grid is a disposable scratchpad — reset to zero at the start of every step and discarded at the end. There is no Lagrangian mesh connectivity between particles, which is exactly what lets MPM handle large deformation, self-collision, fracture, and topology change automatically while still supporting rich, history-dependent elasto-plastic materials. (Verified against the SIGGRAPH 2016 course notes and the Disney 2013 snow paper.)

2. History

MPM did not appear from nowhere — it is the end of a decades-long lineage of particle-in-cell methods (covered in Prerequisites §6).

• 1960s — PIC. Harlow and colleagues at Los Alamos develop the Particle-In-Cell method for fluids: particles carry mass through a fixed grid that does the math. Stable but heavily dissipative.
• 1986 — FLIP. Brackbill & Ruppel introduce the Fluid-Implicit-Particle method, transferring velocity changes rather than velocities to slash dissipation.
• 1994–1995 — MPM is born. Deborah Sulsky, Zhen Chen, and Howard Schreyer at the University of New Mexico generalize FLIP from fluids to solid / history-dependent continuum mechanics. The founding papers are Sulsky, Chen & Schreyer, “A particle method for history-dependent materials” (Comput. Methods Appl. Mech. Engrg., 1994) and Sulsky, Zhou & Schreyer (1995). (Verified, 2-0.) For two decades MPM lived mostly in computational engineering and geomechanics (soil, snow, landslides, impact).
• 2013 — MPM enters computer graphics. Stomakhin, Schroeder, Chai, Teran & Selle, “A Material Point Method for Snow Simulation” (SIGGRAPH 2013), done at Walt Disney Animation Studios, is the landmark. It introduced a user-controllable elasto-plastic snow model and powered the snow in Frozen (2013). Disney’s in-house MPM engine is called Matterhorn (Stomakhin was its lead developer).
• 2015 — APIC. Jiang et al., The Affine Particle-In-Cell Method — the stable, angular-momentum-conserving transfer scheme that modern MPM uses.
• 2016 — the canonical course & sand. The SIGGRAPH 2016 course notes consolidate the method into its standard pedagogical form; Klar et al. bring Drucker-Prager plasticity for sand.
• 2018 — MLS-MPM. Hu et al. fuse Moving Least Squares (MLS) shape functions with APIC for a ~2× speedup and a simpler kernel (see §5).
• 2024 — PB-MPM. EA’s research lab SEED ships Position-Based MPM (PB-MPM), targeting real-time games and making MPM stable at large timesteps (see Implementation).

MPM has been used at Disney in Frozen, Big Hero 6, and Zootopia (stated verbatim in the SIGGRAPH 2016 course abstract), and broadly across graphics research and geotechnical engineering.

3. Why a hybrid? The core idea

Every simulation method makes a trade between tracking material and handling deformation:

• Pure Lagrangian (the Finite Element Method, FEM): mesh glued to material. Excellent accuracy and conservation, but the mesh tangles under large deformation, and fracture, self-collision, and material coupling each need bolt-on machinery. The SIGGRAPH 2016 abstract names exactly this: FEM “usually requires additional computational components in the case of large deformation, mesh distortion, fracture, self-collision and coupling between materials.”
• Pure Eulerian (grid fluids): no mesh to tangle, topology change is free, but material identity and history smear across the grid, and advection is lossy.

MPM splits the responsibilities:

The particles never need neighbor connectivity, so the material can deform, split, and merge arbitrarily; the grid is regular, so derivatives and contact are cheap and robust. From the course notes, the resulting advantages:

• Like FEM, MPM derives from the weak form of conservation of momentum — physically accurate discretization.
• Boundary conditions, wall collisions, and external forces apply easily on the grid and particles.
• Automatic self-collision/contact, because particle motion is interpolated from a single-valued grid velocity field — two pieces of material at the same grid node simply move together.
• Automatic splitting and merging from the particle representation — ideal for fluids and granular materials.
• Multi-material / multiphase coupling for free: give particles different constitutive models; they interact through the shared grid.

The cost (course notes are explicit about this): there is “a sacrifice of some accuracy” — pure hyperelastic solids are not simulated as faithfully as in FEM, and grid contact is famously sticky (no-slip). You trade some accuracy for self-collision and topology change for free.

4. The algorithm — one simulation step

MPM repeats a four-phase cycle each timestep. The grid is reset every step; all persistent state lives on particles.

Phase 1 — Reset grid

Set mass and momentum to zero at every grid node. The grid carries nothing between steps.

Phase 2 — Particle-to-Grid (P2G)

For each particle p, over its 3×3 (2D) or 3×3×3 (3D) B-spline stencil of nearby nodes i:

1. Evaluate interpolation weights $w_{ip}$ from the particle’s position.
2. Compute the particle’s stress from its deformation gradient F via the constitutive model (e.g. fixed corotated; see Prerequisites §3).
3. Scatter mass and momentum (and the affine contribution, with APIC) to the nodes, weighted by $w_{ip}$:
   • node mass $m_i \mathrel{+}= w_{ip}\, m_p$
   • node momentum $(m\mathbf{v})_i \mathrel{+}= w_{ip}\, m_p (\mathbf{v}_p + \mathbf{C}_p(\mathbf{x}_i - \mathbf{x}_p))$
4. Add internal forces (the stress divergence) and external forces (gravity, etc.) to the nodes.

This scatter is the algorithm’s performance heart and its parallelization nightmare: many particles write to the same node simultaneously (the atomic write contention problem — see Implementation).

Phase 3 — Grid update

At each node with nonzero mass:

1. Recover velocity from momentum: $\mathbf{v}_i = (m\mathbf{v})_i / m_i$.
2. Integrate forces forward (symplectic/explicit Euler shown): $\mathbf{v}_i^{\,n+1} = \mathbf{v}_i^{\,n} + \Delta t\, \mathbf{f}_i / m_i$. (Implicit/backward-Euler integration — course notes §11 — is used when stability at larger steps is needed; it solves a linear system via Newton’s method.)
3. Enforce boundary conditions and collision objects directly on the grid velocity field.

Phase 4 — Grid-to-Particle (G2P)

For each particle, over the same stencil (weights are identical since the particle hasn’t moved yet):

1. Gather the updated grid velocities into the particle’s new velocity (and, with APIC, reconstruct its affine matrix C = velocity gradient).
2. Update the deformation gradient using the velocity gradient (discretizing $\dot{\mathbf{F}} = (\nabla\mathbf{v})\mathbf{F}$, course notes §9.4): $\mathbf{F}_p^{\,n+1} = (\mathbf{I} + \Delta t\, \nabla\mathbf{v}_p)\,\mathbf{F}_p^{\,n}$. For plastic materials, the return mapping is applied here.
3. Advect the particle: $\mathbf{x}_p^{\,n+1} = \mathbf{x}_p^{\,n} + \Delta t\, \mathbf{v}_p^{\,n+1}$.

The governing equations underneath

MPM discretizes conservation of mass and conservation of momentum ($\rho \frac{D\mathbf{v}}{Dt} = \nabla\cdot\boldsymbol\sigma + \rho \mathbf{g}$). Mass conservation is automatic — particles carry fixed mass. Momentum balance is discretized in its weak (variational) form: multiply by a test function, integrate by parts, and replace the integral over the domain with a sum over material points (the particles act as quadrature points, analogous to FEM’s Gauss points). The P2G force scatter is precisely the assembled weak-form momentum equation. (Verified against course notes §7, 3-0.)

5. Important variants

APIC (Affine Particle-In-Cell, 2015)

Augments each particle with an affine velocity matrix C, transferring local rotational information that plain PIC dissipates and plain FLIP corrupts. The result is stable and low-dissipation, with exact conservation of angular momentum across transfers. APIC is the de-facto transfer scheme for modern MPM. (See Prerequisites §6.)

MLS-MPM (Moving Least Squares MPM, 2018)

Hu, Fang, Ge, Qu, Zhu, Pradhana & Jiang, A Moving Least Squares Material Point Method with Displacement Discontinuity and Two-Way Rigid Body Coupling (SIGGRAPH 2018). Its contributions:

• A new stress-divergence discretization via Moving Least Squares that lets “all MPM simulations run two times faster than before.”
• It eliminates computing the shape-function gradients ∇N — standard MPM needs both the weights N and their gradients per particle; MLS-MPM reuses the APIC affine matrix as the velocity gradient, so only N is needed. This is the single biggest reason modern real-time MPM work starts from MLS-MPM.
• It naturally derives APIC (and Polynomial PIC) from a Galerkin-style weak form — unifying the transfer schemes with the discretization.
• The companion CPIC (Compatible Particle-In-Cell) algorithm adds material cutting, dynamic open boundaries, and two-way rigid-body coupling via colored distance fields — phenomena standard MPM couldn’t do.
• Famous for its 88-line reference implementation (see Implementation).

Disney snow MPM (2013)

The graphics debut. A user-controllable elasto-plastic constitutive model (fixed-corotated elasticity + singular-value-clamping plasticity with hardening) on the MPM framework, targeting wet/dense snow that behaves as both solid and fluid — a regime prior methods handled poorly. Its semi-implicit integration has conditioning independent of the particle count.

PB-MPM (Position-Based MPM, 2024)

EA SEED’s real-time-oriented reformulation; stable at any timestep. The transfer cycle is wrapped in a constraint-projection loop, so each material becomes “what shape does F want to be, and how fast do we relax toward it” rather than a stress law. What that buys is unconditional stability, not accuracy: stiffness becomes a function of iteration count and Δt, so materials are relaxation rates rather than Young’s moduli and cannot be calibrated to physical units. The close reading of the reference implementation — the displacement form, the per-material projections, the fused tile kernel — is Implementation §3.3.

Transfer filters — XPIC(m) and FMPM(k) (2017, 2021)

A different axis from APIC. Write one G2P→P2G round trip as the operator $SS^{+}$ ($S^{+}$ is the mass-weighted transfer); its complement $I - SS^{+}$ isolates the null space — particle velocity content the grid cannot see, and the source of the noise FLIP is famous for. XPIC(m) (Hammerquist & Nairn 2017) removes it recursively, $P = I - (I - SS^{+})^{m}$, converging to exact null-space removal as $m \to \infty$; $m = 1$ is plain PIC. It is computed by $m-1$ grid round trips, never by forming matrices.

FMPM(k) (Nairn & Hammerquist 2021) reframes the same thing: XPIC(k) is exactly a $k$-term truncated Neumann series for the consistent-mass-matrix inverse, $\tilde{m}^{-1} = (I + A + A^{2} + \dots)\,m^{-1}$ with $A = I - S^{+}S$ and residual $A^{k}$. Seen that way the filtered velocity can replace the lumped one everywhere, strain updates included. Four properties matter for anyone implementing it:

• The incremental form is the one to build. $\Delta v(1) = p^{+}/m$; then $k-1$ times, one gather plus one mass-weighted scatter forms $S^{+}S\,\Delta v$, subtract for the next increment, accumulate. Each order costs one extra P2G+G2P-shaped pair over the same grid, and the increment magnitude doubles as a convergence monitor.
• Constraints must live inside the loop. Zero each increment on BC-controlled nodes and apply contact corrections per increment. Constraining only at the end degrades with order and goes unstable at interfaces; the incremental form reproduces the no-interface result exactly.
• The particle update gets cheaper, not dearer: $V^{+} = S\,v^{+}(k)$ and $X \mathrel{+}= \tfrac{1}{2}(V^{+} + V)\,\Delta t$, with no Lagrangian inconsistency between the position and velocity updates — filtering is what makes using the particle velocity in the position update safe.
• Filtering and APIC are complementary, not substitutes. A filter targets grid-invisible velocity noise; APIC’s affine state is the published fix for the angular-momentum deficit of plain transfers. Neither replaces the other.

Two cautions. Dissipation is not monotone in order — in the source’s impact benchmark FMPM dissipates more than FLIP below order 4 and less at order ≥ 4 — the practical point being that at finite resolution all MPM dissipates; a scheme that gains energy is the red flag, since exact-conservation claims usually hide injection. And multi-field filtering is unresolved: filtering per material field does not reduce to single-field XPIC even across a perfectly bonded interface, and the authors treat it as an open concern in their own implementation.

GIMP and CPDI — shape functions from particle domains (2004, 2011)

GIMP derives shape weights as particle-domain averages of the grid functions, $S_{ip} = (1/V_p)\int \chi_p N_i \,\mathrm{d}V$. uGIMP is the analytic hat⊗box convolution; B2-GIMP convolves the quadratic spline with the particle box; CPDI (Sadeghirad et al. 2011) tracks the deformed domain through F, with corner positions from semi-side vectors $\mathbf{r}_i = \mathbf{F}\,(\text{half-size} \cdot \hat{\mathbf{e}}_i)$ recomputed each step, so no corner state is stored.

The reason this matters for a B-spline solver is mostly confirmation rather than homework: a quadratic B-spline is the same convolution with a fixed cell-width box, so B-spline transfers already carry the C¹ continuity that kills cell-crossing noise — GIMP’s whole reason to exist. Gao et al. 2017 make it exact, proving uGIMP with $L_p = \Delta x$ is identical to the quadratic B-spline. What does not come for free is CPDI’s tracking of a deformed domain, which is a genuine extension — at the cost of corner-sampled transfers that are hostile to a divergence-sensitive GPU stencil.

Multi-field (multimaterial) contact MPM (2000–)

Plain MPM stores one velocity per grid node, so two distinct bodies whose particles reach the same nodes get mass-averaged into a single motion — they weld instead of colliding. The classic cure (Bardenhagen, Brackbill & Sulsky 2000; Bardenhagen et al. 2001) gives each body/material its own velocity field on the shared grid and resolves a contact law at nodes where fields overlap: detect approach + proximity, then apply an inelastic normal match plus Coulomb friction (separable contact) — or deliberately merge the fields (“stick,” which reproduces single-field behavior for that pair). Design details that matter in practice: contact normals from the difference of the fields’ mass/volume gradients (Nairn 2013’s average-gradient method; Nairn, Hammerquist & Smith 2020 supersede gradient-based normals with logistic regression and add null-space filtering); detection must combine approach and a position-separation threshold (approach-only fires ~a stencil-width early); nodes with 3+ fields are handled pairwise (explicit 3+-body treatments perform worse — Nairn’s docs note the nodal fields simply don’t carry enough information). Related lines: multi-mesh MPM (Hu & Chen 2003; Pan et al. 2008) gives each body a whole background grid; field-gradient partitioning (Homel & Herbold 2017) creates the two fields dynamically from a damage gradient for fracture self-contact. Full citations: References Contact & multi-body coupling. Distinct from CPIC (§ MLS-MPM above), which handles rigid/thin boundaries by particle-node coloring rather than by per-body velocity fields.

The calibration half. The architecture above is widely reimplemented; the constants usually are not, and a shipping engineering code has answers for most of them:

• The separation constant is derivable, not tunable. Extrapolated surface positions systematically lag true edges by a shape-function-dependent bias. Inverting a 1D model of an edge approaching a node gives a subtraction of 0.8Δx at 2 particles per cell (1.07Δx at 1 ppc) for GIMP functions, or a power-law inverse with exponent ≈0.58. The method — redo that inversion for whichever shape function and ppc you actually run — is what turns a hand-tuned threshold into a computed one. A second-order refinement advances the separation by the step’s own relative motion, $\delta_n \mathrel{+}= (\Delta \mathbf{p} \cdot \hat{\mathbf{n}})\,\Delta t / m_\mathrm{red}$. For a quadratic B-spline at the standard two particles per axis, a planar edge on a node sees particle rows at 1/4, 3/4, and 5/4 cell with weights 11/16, 9/32, and 1/32. Each extrapolated position therefore lags the true edge by 27/64 cell, so the two-field contact-position offset is 27/32 = 0.84375 cell. The corresponding ideal offsets are 1.0 at one particle per axis and 0.827160 at three; they are packing-dependent predictions, not alternative tuning values for the $q = 2$ default.
• Take normals from volume gradients, not mass gradients, so density-mismatched pairs (an elastic body on sand, anything against a fluid) do not skew the normal toward the denser material. Use one shared pair normal; per-material own normals break momentum conservation. With uniform particle mass the two are proportional and the distinction is moot — it starts to matter the day mixed-density or variable-$J$ pairs ship.
• Prescreen on mass ratio. Skip contact when one field holds less than ~1e-5 of the node’s mass: a nearly-empty field’s momentum is one fringe particle’s noise, and the impulse scales with reduced mass so nothing physical is lost. Cheap, and it fixed a stalled sliding case in the source.
• Cap contact-imposed grid strain, not absolute mass. A light but legitimate field can require a large nodal velocity correction. For the full contact correction $\Delta \mathbf{v}$, define $\dot{\varepsilon}_{\max} = \max_{\mathrm{body},\,\mathrm{component}} |\Delta \mathbf{v}| / \Delta x$ and scale both bodies by $\min\!\left(\gamma / (\dot{\varepsilon}_{\max}\,\Delta t),\; 1\right)$. Bardenhagen et al. 2001 derive this collapsed-cell criterion and recommend $0.5 \le \gamma \le 1$; the shared pair scale preserves equal-and-opposite momentum. This is distinct from the relative fringe screen above.
• Contact area from extrapolated volumes ($A \approx \sqrt{2\,\min(V_a, V_b)\,V_\mathrm{tot}} \,/\, h_\perp$) converts impulse laws into traction laws, which is what makes adhesion expressible at all: hold while $(S/S_a)^{2} + (N/N_a)^{2} < 1$, an elliptical failure surface, optionally with velocity-dependent μ.
• Compliant contact without the penalty trap. Model the interface as bilinear (zero tension, finite compression stiffness) and integrate the resulting 1-DOF oscillator analytically within the step; when $\varphi = \kappa\,\Delta t > \pi/2$ the spring is too stiff for the step, so fall back to ordinary stick. Soft contact that cannot explode. The same machinery with tension stiffness plus a debond threshold is breakable glue.
• Regression normals are the current quality ceiling — fit a separating plane through the near-node particles (linear warm start, then logistic Gauss–Newton), normal from the plane, separation from SVM-style margins with a particle-radius correction. Reported to supersede every gradient option, but the per-node iteration over particle lists is unattractive on a real-time GPU budget.
• The consistency ledger: re-impose BCs after every contact pass (BCs win), and make any post-force momentum change also adjust the nodal force by $\Delta \mathbf{p}/\Delta t$. Validate by demanding that two identical materials with a welded interface reproduce single-field results exactly — a transparency check worth building before trusting any of the above.

6. Use cases — what MPM is good at

MPM’s sweet spot is large-deformation continua and the transitions between material states — exactly the phenomena that defeat mesh-based methods. The SIGGRAPH 2016 course lists demonstrated materials: elastic objects, snow, lava, sand, and viscoelastic fluids.

Beyond graphics and games, MPM is heavily used in geotechnical / civil engineering (slope stability, landslides, soil-structure interaction) and high-velocity impact mechanics — see the open-source CB-Geo MPM and Anura3D engines, also listed in References.

The material axis — constitutive models

A material’s behavior in MPM is set by its constitutive model — the stress↔︎strain law evaluated per particle — which is orthogonal to the transfer scheme (§5, APIC/MLS-MPM/PB-MPM). A solver can carry an arbitrary mix of constitutive models on one grid; adding a material is “add a stress function,” not “add a solver.” Beyond the snow / sand / fluid models that dominate graphics MPM, the models below round out the practical palette — several are the standard MPM demo materials (§6 table) and the rest are well-trodden in the engineering literature.

Implementation shapes that decide cost

The model is orthogonal to the transfer scheme, but the form a model is written in decides whether it is affordable per particle per substep. The shapes below are the ones that matter at real-time budgets, drawn from the engineering-MPM literature (Nairn’s material-model catalog, following de Souza Neto/Perić/Owen) rather than the graphics line.

• Split pressure from deviatoric stress. Track $P$ and $\mathbf{s}$ separately, with the pressure update as an overridable EOS plug-in point and a J2 radial return whose direction is constant during the return — which collapses the whole projection to one scalar equation for $\lambda$. Linear hardening then has a closed-form $\lambda$ and needs no iteration, which is the shape a GPU wants. Nonlinear hardening laws need Newton with bracketing (start at rate 1/s, ×10 until the sign flips) — plain Newton diverges for some laws.
• The volumetric energy $U(J)$ is a real choice, not a detail. Three common forms, with a sharp verdict: $U = (K/2)\left(\tfrac{1}{2}(J^{2}-1) - \ln J\right)$ is the only one that behaves at both extremes ($P \to \infty$ as $J \to 0$, $P \to -\infty$ as $J \to \infty$). The common graphics choice $U = (K/2)(J-1)^{2}$ — the $\lambda (J-1) J \mathbf{I}$ term in fixed-corotated — has a tangent bulk modulus that vanishes under crush ($K(J) = KJ \to 0$ as $J \to 0$): particles go soft exactly when compacted hardest. The forms only diverge far from $J = 1$, i.e. impacts and crushes.
• Von Mises without an SVD. The cheapest credible “holds its dent” material is Neo-Hookean elasticity plus a J2 return taken directly on the deviatoric Kirchhoff stress $\mathbf{s} = \mu\,\mathrm{dev}(\bar{\mathbf{B}})$ with $\bar{\mu} = \mu\,\bar{I}_1/3$, updating the elastic left Cauchy–Green tensor as $\bar{\mathbf{B}} \leftarrow \bar{\mathbf{B}}_\mathrm{trial} - \tfrac{2}{3}\lambda \bar{I}_1 \mathbf{n}$. State: 8 scalar components: 6 for symmetric algorithmic $\bar{\mathbf{B}}$, 1 for authoritative total $J$, and 1 for $\alpha$. No eigendecomposition anywhere, closed form for linear hardening.
• Viscoelasticity as a Prony series with recursive internal variables — no strain history stored. The Maxwell-element form is the simplest to implement: $\mathbf{s} = \mathbf{s}_0 + \sum_k \mathbf{s}_k$ with a per-element ODE $\mathrm{d}\mathbf{s}_k = 2 G_k \,\mathrm{d}\mathbf{e} - (\mathbf{s}_k/\tau_k)\,\mathrm{d}t$ — n deviatoric tensors, each an exponential decay plus a drive term. Dissipation $\sum_k \mathbf{s}_k^{2} / (2 G_k \tau_k)$ is the energy sink that settles a gel. Objectivity constraint: stored internal tensors must rotate with the particle, the same treatment stored stress needs.
• Damage with resolution-independent dissipation. 1D theory (strength envelope $F(\delta)$, damage $D = \delta/(\delta + F(\delta)/E)$), evolution costing one divide per step, and crack-band scaling: solve $\delta_{\max}$ from fracture toughness via $G_{Ic} = (V_p/A_c)\int F(\delta)\,\mathrm{d}\delta$, which makes dissipated energy per particle $G_c A_c$ regardless of resolution. Two pieces of explicit guidance worth carrying: chasing the “true” softening-law shape is a fool’s errand — linear softening is the stable recommendation — and there is a mesh limit $\Delta x < \eta E G_c / \sigma_0^{2}$ ($\eta = 2$ for linear) below which softening snaps back and explodes. Scalar isotropic damage (one $D$, regula-falsi only on damaging particles) is the budget tier; phase-field fracture is explicitly not real-time (an extra grid PDE, cells $\le \ell/2$).
• Every law reports its dissipated energy $\mathrm{d}\Phi \ge 0$ through one interface — plastic work, viscous flow, damage. That is both the physically correct settle path and, if routed, the source term a “friction melts snow” feature needs.

Two liquid-pressure cautions belong with this list. True Tait ($P = C K_0 \left(e^{(1-J^{*})/C} - 1\right)$, $C \approx 0.0894$) and the Murnaghan/Monaghan power law ($P = (K/\gamma)\left(J^{-\gamma} - 1\right)$) agree for $J^{*} > 0.7$ and then diverge badly, with true Tait far stiffer under deep compression; its tangent modulus $K(P) = (P + C K_0) J^{*} / C$ grows with pressure, so a CFL estimate should follow the tangent modulus rather than the nominal one. And general EOS forms carry a shared trap — Mie-Grüneisen’s pressure denominator vanishes at a finite compression — so clamp EOS stiffness and compression, and warn when the clamp engages.

References for these constitutive models: the SIGGRAPH 2016 MPM course (Jiang et al., in References — derives elastic, plastic, and viscoelastic models), nialltl’s Neo-Hookean walkthrough, and cebas thinkingParticles’ MPM operator docs as a concrete per-material parameter reference: Visco, VonMises, Elastic, Linear.

7. Limitations & pitfalls

MPM is powerful but not free. The well-documented weaknesses, especially acute for real-time use:

• Small timestep / CFL restriction. Explicit MPM needs a small Δt for stability — the Courant–Friedrichs–Lewy (CFL) condition: a particle shouldn’t cross more than a fraction of a cell per step. This is the barrier to real-time use and the problem PB-MPM and implicit integration attack.
• High memory cost — bandwidth, not capacity, is what bites. MPM stores a lot (many particles each carrying F, C, mass/volume and material state, plus the grid), but raw capacity is rarely the binding constraint on a modern GPU — it bites mainly on memory-tight targets (mobile / console shared RAM). The real-time weakness is memory bandwidth / traffic: every substep streams all of that per-particle and per-grid state through the P2G and G2P transfers, and at the CFL-limited step count that traffic — not the particle count sitting in memory — is the cost. (The scatter half of that traffic is also a contention problem; see the next bullet.) Practical implication: optimizing MPM means cutting data movement per step (compact state, fused passes, cache/tile-friendly layouts), not just fitting more particles in RAM.
• P2G scatter is the parallelization bottleneck. Concurrent writes to shared grid nodes cause atomic contention — the dominant GPU optimization target (see Implementation).
• Numerical fracture. Material can split apart spuriously when particles in a cell drop below the count needed to represent the continuum — sometimes a feature (free fracture), often an unwanted artifact.
• Cell-crossing instability. Force discontinuities as a particle crosses a cell boundary; mitigated (not eliminated) by smooth B-spline kernels.
• Sticky / no-slip grid contact. The automatic contact is famously no-slip; true frictional/separating contact needs extra modeling — multi-velocity-field contact (§5) for body-vs-body, or CPIC for rigid/thin boundaries. Corollary: distinct same-material bodies sharing one velocity field don’t just stick, they weld into one connected body under sustained contact (a resting stack), even when a transient collision bounces cleanly.
• Volume/quadrature error & particle-count sensitivity. Particles double as quadrature points; too few per cell degrades accuracy (ringing, “self-collision” artifacts).
• Energy dissipation in transfers. Mitigated by APIC/MLS-MPM, but the PIC↔︎FLIP dissipation–noise trade-off is intrinsic to the family.

8. MPM vs. neighboring methods

The methods MPM is most often weighed against: the Finite Element Method (FEM), Smoothed Particle Hydrodynamics (SPH), the FLIP/PIC particle-in-cell family, Position-Based Dynamics (PBD), and — for water specifically — height-field / shallow-water fluids, the real-time incumbent. What each method is on its own terms — how it works, where it ships — lives in Related Methods; this section is the MPM-relative comparison. The chart places them by the two axes that matter most.

The height field sits in the far corner away from MPM, and that is the honest picture: it gives up material identity and the ability to represent large-deformation geometry (no overhang, no overturning wave) — which is exactly why it is cheap enough to have won real-time water.

vs. FEM (Finite Element Method)

Both derive from the weak form of momentum balance and both excel at elastic solids. FEM is more accurate for moderate-deformation elasticity (no transfer dissipation, no quadrature/particle-count issues). MPM wins on large deformation, fracture, self-collision, and multi-material coupling, which FEM needs bolt-on remeshing/contact machinery to handle. Rule of thumb: FEM for a deforming-but-intact solid; MPM when it tears, flows, or mixes.

vs. SPH (Smoothed Particle Hydrodynamics)

Both are particle methods with no mesh, so both handle topology change easily. But SPH computes derivatives via a kernel sum over neighbors, requiring an expensive per-step neighbor search; MPM computes derivatives on the fixed background grid and needs no neighbor search. SPH is also prone to tensile instability (particles clump under tension); MPM’s grid-node field approximation sidesteps this. MPM generally gives better-behaved stress and contact; SPH is simpler to implement for pure fluids.

vs. FLIP / PIC

MPM is the generalization of FLIP/PIC to solid continuum mechanics (it adds the deformation gradient and constitutive stress). Plain FLIP/PIC handle only (mostly inviscid) fluids; MPM handles arbitrary elasto-plastic continua. The transfer-dissipation lessons (PIC dissipative, FLIP noisy, APIC balanced) carry straight over.

vs. height-field / shallow-water fluids

The comparison that matters if the material in question is water, because this is what real-time water actually uses. A height field collapses the liquid to one height per column and solves the depth-averaged shallow water equations on a 2D grid, so cost scales with surface area rather than volume — complete scenes in single-digit milliseconds, versus a 3D MPM domain that is orders of magnitude larger for the same lake. Overturning waves, waterfalls, and splashes — which the representation cannot express at all — are patched by spawning particles where the field is detected to be failing and reabsorbing them on landing. For a river, lake, or shoreline, MPM does not beat this on cost and should not claim to. MPM’s case is the four things the height field structurally forgoes: material variety (one solver for snow, sand, mud, and mixtures rather than water plus bespoke extensions), phase change and mixing, genuine 3D flow (caves, overhangs, arbitrary topology), and material identity/history. Choosing MPM for water alone is choosing the expensive tool; choosing it because the water must also be snow, or must mix with sand, is where it pays for itself.

vs. PBD (Position-Based Dynamics)

PBD enforces geometric constraints directly on positions — fast, stable, and ubiquitous in real-time games, but not physically grounded (stiffness is iteration- and timestep-dependent; it approximates rather than solves the governing equations). MPM is continuum-accurate but heavier. The course notes note MPM “provides a unified particle simulation framework similar to Position Based Dynamics (PBD).” Notably, EA’s PB-MPM bridges the two — borrowing PBD’s compliant-constraint stability to make MPM real-time-stable (see Implementation).

* The real-time barrier is what MLS-MPM (speed) and PB-MPM (stability at any timestep) directly target.

† Water-only, and included precisely because of that: the N/A cells are what buys the cost. See Related Methods §1.5.

Where to go next

• Implementation — GPU optimization of the P2G bottleneck, and the reference codebases (taichi_mpm, nialltl, EA PB-MPM).
• Games — whether any commercial game uses MPM, how novel an MPM-native game would be, and the commercial landscape.
• References — every paper cited above, with URLs.
• Prerequisites — the math foundations, if any step above was opaque.

Sources for this document

• Jiang, Schroeder, Teran, Stomakhin, Selle. The Material Point Method for Simulating Continuum Materials, SIGGRAPH 2016 Course Notes — PDF · mpm.graphics. (Primary source for the definition, algorithm, governing equations, and method comparisons.)
• Stomakhin, Schroeder, Chai, Teran, Selle. A Material Point Method for Snow Simulation, SIGGRAPH 2013 — Disney Research.
• Hu et al. A Moving Least Squares Material Point Method…, SIGGRAPH 2018 — project page.
• Sulsky, Chen, Schreyer (1994); Jiang et al. APIC (2015); Klar et al. Sand (2016). Full citations in References.
• de Vaucorbeil et al. Material Point Method after 25 Years survey, and CB-Geo LearnMPM — see References.