Prerequisites — The Math & Physics You Need First

The Material Point Method (MPM) sits at the intersection of continuum mechanics, numerical PDE discretization, and the particle-in-cell family of fluid solvers. You do not need a graduate degree to follow it, but a handful of concepts are genuinely essential — skip them and the algorithm reads as a pile of unmotivated index gymnastics. This document builds those concepts from the ground up.

The canonical reference these notes track is the SIGGRAPH 2016 course “The Material Point Method for Simulating Continuum Materials” by Jiang, Schroeder, Teran, Stomakhin, and Selle (PDF, mpm.graphics). Its own stated prerequisites are modest: “minimal concepts of continuum mechanics. Familiarity with multivariable calculus, linear algebra and common numerical algorithms is assumed.” This document fills in exactly those.

Concept dependency map

How the prerequisite ideas build on one another, and where each feeds into the MPM algorithm:

1. Lagrangian vs. Eulerian frames

This is the single most important conceptual split in all of computational physics, and MPM’s defining trick is that it uses both at once.

• Lagrangian frame — you ride along with the material. Your coordinate is the material (undeformed) position X, fixed to a given chunk of stuff for all time. A mesh-based solid solver — the Finite Element Method (FEM), the standard technique for simulating deforming solids — is Lagrangian: each node is glued to a material point and moves with it. Quantities like the velocity V(X, t) are functions of which particle you are watching.
• Eulerian frame — you stand at a fixed point in space x and watch whatever material flows past. A grid-based fluid solver is Eulerian. The velocity v(x, t) is a function of where in space you are looking, regardless of which particle currently occupies that spot.

The two views are linked by the flow map (a.k.a. deformation map) φ(X, t), which sends a material point X to its current world position x:

$$\mathbf{x} = \phi(\mathbf{X}, t).$$

Why MPM cares. Lagrangian methods track material naturally (mass is exactly conserved, history-dependent material state rides along on each particle) but their mesh tangles and degrades under large deformation. Eulerian methods handle arbitrary deformation and topology change effortlessly (there is no mesh to tangle) but smear material identity and history across a fixed grid. MPM keeps material state on Lagrangian particles and borrows a disposable Eulerian grid only to compute spatial derivatives (forces) each step — getting the conservation and history-tracking of the former with the large-deformation/contact robustness of the latter.

Material derivative. Because a fixed spatial point sees different particles over time, the time-rate-of-change of a field following the material is not just ∂/∂t. The material derivative D/Dt adds an advection term: $$\frac{D}{Dt} f = \frac{\partial f}{\partial t} + (\nabla f)\cdot \mathbf{v}.$$ This is why Eulerian fluid solvers need an advection step and Lagrangian ones don’t. (Course notes §5.4.)

2. Continuum mechanics essentials

MPM treats matter under the continuum assumption: material is modeled as a continuous medium with smooth fields (density, velocity, stress) rather than discrete molecules. This holds for essentially everything graphics simulates — elastic solids, flesh, cloth, sand, snow, mud, water, smoke. Crucially, an MPM “particle” is not a molecule — it is a small chunk of continuum, a quadrature (sampling) point for the material domain, analogous to a Gauss point in FEM.

2.1 The deformation gradient F — the star of the show

The deformation gradient is the Jacobian of the flow map. It is the local linear map from a tiny material fiber to its deformed world-space counterpart:

$$\mathbf{F}(\mathbf{X}, t) = \frac{\partial \phi}{\partial \mathbf{X}}(\mathbf{X}, t) = \frac{\partial \mathbf{x}}{\partial \mathbf{X}}.$$

Intuitively, take two nearby material points; F maps the vector between them in the rest state to the vector between them now. It captures all the local stretching, shearing, and rotation at a material point. In 2D it is a 2×2 matrix, in 3D a 3×3 matrix.

Key facts (course notes §5.2):

• Rigid motion ⇒ $\mathbf{F} = \mathbf{R}$ (a pure rotation matrix). No stretch, no internal force. The identity I is the special case of no motion at all.
• $J = \det(\mathbf{F})$ is the local volume ratio (deformed volume ÷ rest volume).
  • $J > 1$ → locally expanded, $J < 1$ → compressed.
  • $J = 1$ → volume-preserving (all rigid motions, and incompressible flow).
  • $J \to 0$ → material crushed to zero volume; $J < 0$ → inverted (a triangle turned inside-out). Robust solvers must tolerate these degeneracies.
• F evolves by the velocity gradient (course notes Eq. 24, the relation that drives the per-particle F update in MPM): $$\dot{\mathbf{F}} = (\nabla \mathbf{v})\,\mathbf{F}.$$

Everything elastic in MPM is a function of F: the strain energy, the stress, and therefore the force.

2.2 Stress: Cauchy and Piola-Kirchhoff

Stress is force per unit area transmitted through the material. Two measures matter for MPM:

• Cauchy stress σ — force per current (deformed) area. This is the “true,” physical stress in world space. Symmetric for non-polar materials (conservation of angular momentum).

• First Piola-Kirchhoff stress P — force per rest area. It is the natural object when your energy is written in terms of F, because

  $$\mathbf{P} = \frac{\partial \Psi}{\partial \mathbf{F}},$$

  where Ψ(F) is the strain-energy density. P and σ are related by $\sigma = \tfrac{1}{J}\mathbf{P}\mathbf{F}^{T}$. MPM’s force computation is cleanest in terms of P (course notes §6.1), which is why hyperelastic models below are given as Ψ and its derivative P.

Strain is the dimensionless measure of how much the material has deformed (distinct from stress, the force response). For small deformations the linear (infinitesimal) strain tensor suffices; for the large deformations MPM targets, the full nonlinear kinematics of F are required, which is why MPM works directly with F rather than a linearized strain.

3. Constitutive models — turning deformation into force

A constitutive model is the material’s “personality”: the rule mapping deformation (F) to stress (P). MPM’s flexibility comes from the fact that this is just a per-particle function — change it (or store different parameters per particle) and you get a different material, with no change to the solver core. This is how one MPM framework simulates elastic, snow, sand, and fluid in the same scene.

3.1 Hyperelasticity (energy-based elasticity)

A hyperelastic material derives its stress from a stored strain-energy density Ψ(F), so $\mathbf{P} = \partial\Psi/\partial\mathbf{F}$. Two workhorse models (course notes §6.2–6.3):

Neo-Hookean. A simple nonlinear elastic model. With Lamé parameters μ (shear stiffness) and λ (resistance to volume change), and $J = \det(\mathbf{F})$:

$$\mathbf{P}(\mathbf{F}) = \mu\left(\mathbf{F} - \mathbf{F}^{-T}\right) + \lambda \log(J)\,\mathbf{F}^{-T}.$$

(This is exactly the form nialltl’s tutorial codes for elastic bodies — see Implementation.) It is cheap but blows up as $J \to 0$ because of the $\log(J)$ and $\mathbf{F}^{-T}$ terms.

Fixed Corotated. The graphics-favored model (introduced by Stomakhin et al. 2012 and used in the Disney snow paper). It penalizes deviation from the nearest rotation R, so it stays well-behaved under large rotations and inversion:

$$\mathbf{P}(\mathbf{F}) = 2\mu(\mathbf{F} - \mathbf{R}) + \lambda(J - 1)\,J\,\mathbf{F}^{-T},$$

where R is the rotation factor of F. Getting R requires the polar decomposition / SVD (§4 below).

3.2 Plasticity — permanent deformation (snow, sand, mud)

Elastic materials spring back. Plastic materials deform permanently past a yield threshold — and that is precisely what makes snow pack, sand pile, and mud slump. The standard trick is a multiplicative split of the deformation gradient into an elastic and a plastic part:

$$\mathbf{F} = \mathbf{F}_E\,\mathbf{F}_P.$$

Only the elastic part F_E generates stress; the plastic part absorbs permanent deformation. A return mapping projects F_E back onto an admissible (yield) region whenever stress exceeds the yield criterion.

• Snow plasticity (Stomakhin et al. 2013, course notes §6.5): the singular values of F_E are clamped to a critical compression/stretch range; hardening makes packed snow stiffer. This single mechanism produces fracture, packing, and the wet-vs-powdery spectrum.
• Drucker-Prager (Klar et al. 2016; Daviet & Bertails-Descoubes 2016): the pressure-dependent yield criterion from soil mechanics, the engineering-favored model for sand and granular flow. The yield surface is a cone in principal-stress space — material can shear freely under low confining pressure but resists under compression, exactly how a sand pile behaves. The return mapping is cleanest in logarithmic (Hencky) strain space: take the SVD of F_E, work with the log of the singular values, split into volumetric (trace) and deviatoric parts, and project — if the deviatoric magnitude exceeds the friction cone for the current pressure, scale it back to the cone; pure tension (positive trace) collapses to the cone apex (the material can’t sustain it). The cone apex sits at the stress origin, so cohesionless granular has zero shear strength at a free surface (no confining pressure → it flows like a frictional liquid there).
  • Cohesion (Klar et al. 2016, §cohesion): translate the cone along the hydrostatic axis into tension by an amount $c \ge 0$. The material then sustains tension/shear up to $c$ before yielding — grains stick together and to surfaces — while still flowing plastically under load. $c = 0$ recovers cohesionless sand; larger $c$ models damp/sticky granular (wet sand, packed snow). It’s a one-parameter shift of the same return mapping, not a new model.
  • The apex volume-gain artifact + volume correction (Tampubolon et al. 2017, §4.3.4). The apex projection is a silent trap: it zeroes net-tension strain (correct — cohesionless granular carries no tension, which is what lets an impact disperse it), but in doing so it forgets the expansion — the particle ends up stress-free at its new, spread-out spacing, and any return toward its original density is penalized elastically. Dispersed material therefore re-settles permanently fluffed: under-dense beds and gap-riddled craters after any energetic event, even though initial settling (compression-only — the apex never fires) packs perfectly. Two signature reads: displaced grains that appear to repel each other on landing, and thin stress-free “roof” lines over excavated voids that masquerade as overhangs able to carry weight — which cohesionless sand cannot physically hold. The standard cure: track per particle the log-volume discarded at apex projections (one scalar, $v_c$, init 0) and add $(v_c/d)\,\mathbf{I}$ to the trial Hencky strain before each projection. While debt is outstanding the particle still reads at/beyond the apex, so recompression is stress-free until the pre-dispersal density is restored (the debt drains as the bed compacts); the elastic/shear branches then clear the residual into $\mathbf{F}_E$. Klar et al. 2016 §7.2 covers the complementary shear-branch rule (the cone projection must preserve the trace — non-associative flow); the apex debt is the Tampubolon extension.

4. The linear algebra: SVD & polar decomposition

The corotated and plasticity models above both need to extract a rotation from F and to manipulate its singular values. Two closely related factorizations do this (course notes §6.4):

• Singular Value Decomposition (SVD): any matrix factors as $\mathbf{F} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{T}$, where U and V are rotations (orthogonal) and Σ is diagonal with the non-negative singular values σ₁, σ₂(, σ₃) — the principal stretches along principal axes. Clamping these singular values is how snow/sand plasticity is applied.
• Polar decomposition: $\mathbf{F} = \mathbf{R}\mathbf{S}$, where R is the nearest rotation (the R in the corotated model) and S is a symmetric stretch. It is one short step from the SVD: $\mathbf{R} = \mathbf{U}\mathbf{V}^{T}$ and $\mathbf{S} = \mathbf{V}\boldsymbol{\Sigma}\mathbf{V}^{T}$.

Because the matrices are tiny (2×2 or 3×3), specialized fast closed-form/iterative SVD routines are used in practice rather than general LAPACK calls (see McAdams/Selle’s fast 3×3 SVD, and Higham’s note on SVD via polar decomposition). This per-particle SVD is often the most expensive scalar kernel in an MPM step.

5. Interpolation & B-spline shape functions

MPM constantly moves data between scattered particles and grid nodes. The weights for that transfer are interpolation (shape) functions N(x) centered on each grid node — the same idea as FEM basis functions, but on a regular Cartesian grid.

MPM almost always uses B-spline kernels rather than the hat (linear) functions of classic PIC, because B-splines are smooth (continuous gradients), which kills a notorious MPM artifact called cell-crossing instability — the force jolt a particle feels when it crosses a grid-cell boundary under non-smooth weights.

• Quadratic B-spline: support of 3 nodes per axis (a 3×3 stencil in 2D, 3×3×3 in 3D). This is the most common choice; nialltl’s and most real-time implementations use it.
• Cubic B-spline: support of 4 nodes per axis (4×4 / 4×4×4), smoother still, more expensive.

Each particle therefore “talks to” 9 (2D) or 27 (3D) grid nodes per step. The weights N and, for classic MPM, their gradients ∇N are evaluated from the particle’s position relative to the grid. (One of the contributions of MLS-MPM — Moving Least Squares MPM, covered in Theory — is eliminating the explicit ∇N evaluation.)

6. The PIC / FLIP / APIC family — MPM’s direct ancestors

MPM is, historically and mechanically, the generalization of the particle-in-cell fluid methods to solid continuum mechanics. Understanding this lineage demystifies the transfer steps. All of them share the same scatter-to-grid / solve-on-grid / gather-to-particle rhythm; they differ only in what gets transferred and how much numerical dissipation results.

APIC is the modern default and the one MPM tutorials use. Each particle carries, in addition to velocity, a small affine velocity matrix C that records the local velocity gradient. Transferring this extra information preserves the rotational (angular-momentum) content that plain PIC throws away and plain FLIP corrupts. APIC is the bridge to MLS-MPM (see Theory), which shows the affine matrix and the MPM velocity gradient are the same object, fusing the two and simplifying the math.

Where to go next

You now have the vocabulary. Continue to:

• Theory — the history, the full particle-to-grid → grid → grid-to-particle algorithm with the math, use cases, limitations, and how MPM compares to other simulation methods (smoothed-particle hydrodynamics, finite elements, position-based dynamics).
• Implementation — GPU optimization and the reference codebases (taichi_mpm, nialltl, EA’s position-based MPM).
• References — the primary papers behind every claim above.

Sources for this document

• Jiang, Schroeder, Teran, Stomakhin, Selle. The Material Point Method for Simulating Continuum Materials, SIGGRAPH 2016 Course Notes. PDF · mpm.graphics — §5 Kinematics, §6 Hyperelasticity, §8 Material Particles. (Primary; the source of the deformation-gradient, stress, constitutive-model, and SVD material here.)
• CB-Geo. LearnMPM educational notes. geoelements.org/LearnMPM
• Jiang, Schroeder, Selle, Teran, Stomakhin. The Affine Particle-In-Cell Method, SIGGRAPH 2015 — origin of APIC. (See References.)
• Klar et al. Drucker-Prager Elastoplasticity for Sand Animation, SIGGRAPH 2016 — the sand constitutive model. (See References.)
• N. Higham. Faster SVD via Polar Decomposition (2015). nhigham.com
• nialltl. MLS-MPM guide — practical Neo-Hookean and B-spline weighting. nialltl.neocities.org/articles/mpm_guide