⚙️ Physics Simulation · Continuum Mechanics
📅 April 2026⏱ 15 min🔴 Advanced

Material Point Method (MPM): Snow, Sand & Fracture Simulation

Snowflakes that crunch underfoot, sand dunes that collapse realistically, and glass that shatters with physically accurate crack patterns — all three were made possible in film and games by the Material Point Method. MPM is a hybrid continuum simulation technique that represents material as Lagrangian particles carrying physical state, while solving the equations of motion on an Eulerian background grid. The result is a method that handles large deformations, fragmentation, and phase transitions that would tear apart an FEM mesh and flood through an SPH particle simulation.

1. Why MPM? Limitations of FEM and SPH

The two traditional approaches to continuum simulation both have fundamental limitations that MPM was designed to overcome:

Finite Element Method (FEM)

FEM discretises the material into a mesh of elements and solves the PDE on that mesh. The physical accuracy is excellent. The problem: the mesh is tied to the material's topology. When a material fractures, tears, or undergoes extreme deformation (stretching by 10× or more), the mesh elements become inverted (negative Jacobian), and the simulation fails catastrophically. Remeshing at every fracture event is expensive and requires complex topological surgery.

Smoothed Particle Hydrodynamics (SPH)

SPH is Lagrangian throughout — particles carry all state and interact through kernel-weighted sums with their neighbours. No mesh needed. Large deformations are handled naturally. The problem: accurately capturing solid mechanics (stress, strain, wave propagation) requires much higher particle densities than fluid simulation, and tensile instability causes SPH particles to clump together under tension, producing non-physical fracture patterns.

MPM's key insight: Use particles for transport (they carry mass, velocity, stress — no mesh distortion possible) and use a background grid for solving the equations of motion (clean spatial derivatives, natural collision handling, easy boundary conditions). The grid is reset to zero every timestep, so it never accumulates distortion. The particles remember the material's history.

MPM was originally developed in 1994 by Sulsky, Chen, and Schreyer as a generalisation of the Fluid-Implicit-Particle (FLIP) method to solid mechanics. It remained a research method until Disney Research developed the elastoplastic snow model for the film Frozen (2013), and published the algorithm in detail — triggering enormous interest in the computer graphics community.

2. The MPM Algorithm: Particle-Grid-Particle

Each MPM timestep consists of three phases — transferring state from particles to the grid, solving physics on the grid, and transferring updated state back to the particles:

// MPM timestep: Particle → Grid → Particle // === PHASE 1: Particle to Grid (P2G) === // Reset grid for all grid nodes i: m_i = 0; v_i = 0 // Splat particle mass and momentum onto grid via interpolation weights for each particle p: for each grid node i in p's support (typically 3×3×3 nodes): w = weight(x_p, x_i) // interpolation weight m_i += w * m_p v_i += w * m_p * v_p // weighted momentum f_i -= w * V_p * σ_p ∇w // stress divergence → grid force // Compute grid velocity from accumulated momentum for all grid nodes i with m_i > 0: v_i = v_i / m_i // momentum → velocity // === PHASE 2: Grid Solve === for all grid nodes i: v_i_new = v_i + Δt * (f_i / m_i + g) // solve Newton (g = gravity) apply boundary conditions on v_i_new // === PHASE 3: Grid to Particle (G2P) === for each particle p: v_p_new = 0; ∇v_p = 0 for each grid node i in p's support: w = weight(x_p, x_i) v_p_new += w * v_i_new // gather velocity ∇v_p += v_i_new ⊗ ∇w // velocity gradient (for deformation) // Update deformation gradient F F_p = (I + Δt * ∇v_p) · F_p // F tracks cumulative deformation // Apply plasticity / constitutive update to F_p // Advect particle position x_p += Δt * v_p_new

The grid serves purely as a computational scratch-pad: it is initialised to zero at the start of each step and discarded at the end. All material history (deformation gradient F, plastic strain, density) lives on the particles, which persist across frames. This separation is what makes topology change trivial — a fracturing solid is just particles that stop interacting, with no connectivity table to update.

3. Interpolation: Linking Particles to Grid

The interpolation weights w(x_p, x_i) are the core of the MPM algorithm — they determine how much each particle contributes to each grid node and vice versa. Standard MPM uses B-spline kernels defined on the background grid's cell spacing h:

// 1D quadratic B-spline weight (most common for MPM) // For distance |d| = |x_p - x_i| / h (grid units): w(d) = 0.75 − d² for |d| < 0.5 w(d) = 0.5 * (1.5 − |d|)² for 0.5 ≤ |d| < 1.5 w(d) = 0 for |d| ≥ 1.5 // 3D: w(x_p, x_i) = w(dₓ) * w(d_y) * w(d_z) (separable) // Support radius: ±1.5 cells → 3×3 nodes in 2D, 3×3×3 in 3D (27 nodes) // Weight gradient (needed for force computation): ∇w = (dw/dₓ * w_y * w_z, w_x * dw/d_y * w_z, w_x * w_y * dw/d_z) / h

Higher-order cubic B-splines give smoother results (used in some film productions) but expand the support to 4×4×4 nodes (64 nodes per particle in 3D), roughly tripling the P2G/G2P cost. The quadratic B-spline is the standard choice for most simulations.

Grid resolution tradeoff: Finer grid cells give higher spatial accuracy but require proportionally more particles (target: 4–8 particles per grid cell) to avoid holes in the representation. A 128³ grid with 8M particles is a typical mid-resolution MPM simulation. Film-quality snow in Frozen used particle counts in the hundreds of millions.

4. Constitutive Models: Snow, Sand, Elastic

The constitutive model determines what material the simulation represents by defining how stress σ relates to deformation F. Each particle carries its own deformation gradient F (a 3×3 matrix representing the local stretch and rotation relative to the rest configuration), and the constitutive model maps F to a Cauchy stress tensor σ.

Neo-Hookean Elasticity (Rubber, Soft Tissue)

// Neo-Hookean hyperelastic model // F = deformation gradient (updated each step as F ← (I + ΔtL)·F) J = det(F) // volume change ratio // Piola-Kirchhoff stress (first PK): P = μ(F − F⁻ᵀ) + λ log(J) F⁻ᵀ Cauchy stress: σ = (1/J) P Fᵀ Material parameters: μ = shear modulus (stiffness to shape change) λ = bulk modulus (stiffness to volume change) Both relate to Young's modulus E and Poisson's ratio ν: μ = E / (2(1+ν)), λ = Eν / ((1+ν)(1−2ν))

Elastoplastic Snow Model (Stomakhin et al. 2013, Disney)

// Snow: elastic deformation + plastic flow via singular value decomposition // Polar decomposition: F = R·S (rotation × stretch) // SVD: F = U·Σ·Vᵀ where Σ = diag(σ₁, σ₂, σ₃) are principal stretches // Elastic trial step: treat all deformation as elastic F_E_trial = (I + ΔtL)·F_E_prev // Plasticity: clamp singular values to [1-θ_c, 1+θ_s] // θ_c = critical compression (snow compacts: ~2.5%) // θ_s = critical stretch (snow separates: ~7.5%) SVD(F_E_trial) → U, Σ, Vᵀ Σ_clamped = clamp(Σ, 1−θ_c, 1+θ_s) F_E = U·Σ_clamped·Vᵀ // elastic part (clamped) F_P = Vᵀ·Σ⁻¹·U^ · F_E_trial // plastic part absorbs the residual J_E = det(F_E) // Modified Neo-Hookean stress with plastic hardening: μ_eff = μ · exp(ξ(1 − J_P)) // ξ = hardening coefficient P = 2μ_eff(F_E − F_E⁻ᵀ) + λ(J_E−1)·J_E·F_E⁻ᵀ

Drucker-Prager Sand Model

// Sand: granular, cohesionless, pressure-dependent yield // Yield condition (Drucker-Prager): f(σ) = ‖dev(σ)‖ + sin(φ) · tr(σ) ≤ 0 dev(σ) = σ − (1/3)tr(σ)I (deviatoric stress) φ = friction angle of the sand (~30° for dry sand) // If f > 0 (yielded): project stress back to the yield surface // — see Klar et al. 2016 for the return mapping algorithm // In compression: allow compaction (particles can overlap in the grid) // In tension: particles separate freely (no tensile cohesion)

5. Fracture and Topology Change

Fracture is where MPM truly shines compared to mesh-based methods. Since there is no connectivity structure — only particles and a background grid — topology change is automatic: regions of high tension simply cease to transmit stress when material separates.

// Fracture via damage mechanics (simplest approach) // Each particle carries a damage variable d ∈ [0, 1] // Damage growth when tensile strain exceeds threshold: if (ε_max > ε_crit): d += Δd // accumulate damage if d > 1.0: d = 1.0 // fully broken // Apply damage to stress: σ_effective = (1 − d) · σ_elastic // Fully damaged particles (d→1): no longer transmit tensile stress // They still contribute mass and can collide as fragments // Crack path emerges automatically from the damage field topology

More sophisticated models use cohesive zone formulations, where particles in a fracture zone experience reduced tensile strength following a specific traction-separation law (linear, exponential, or bilinear). The crack front automatically follows the path of maximum stored elastic energy, producing realistic fracture patterns without any explicit crack tracking.

Multi-material MPM simulations (fluid + solid + gas simultaneously) handle phase interfaces naturally: each material has its own particles and constitutive model, and the grid naturally mediates interactions — a solid object impacting a fluid just adds its momentum splat to the same grid, and the grid velocity update reflects the combined interaction.

6. APIC and MLS-MPM: Modern Variants

APIC: Affinely Particle-In-Cell (Jiang et al. 2015)

Standard MPM (and FLIP) transfer schemes suffer from numerical dissipation (energy loss each timestep) or noise amplification. APIC solves this by having each particle carry not just a velocity but a local affine velocity field — a 3×3 matrix C_p that represents the linear velocity gradient in the particle's neighbourhood:

// APIC: each particle stores (v_p, C_p) — velocity + affine matrix // P2G transfer (with affine correction term): momentum_i += w_ip * m_p * (v_p + C_p · (x_i − x_p)) // ↑ affine correction: more accurate momentum splat // G2P transfer: v_p = Σᵢ w_ip * v_i_grid C_p = (4/h²) Σᵢ w_ip * v_i_grid ⊗ (x_i − x_p) // C_p carries the velocity gradient → no need for separate ∇v computation // Benefits: exactly conserves angular momentum (fixes the FLIP rotation leak) // zero numerical dissipation without FLIP-style noise

MLS-MPM: Moving Least Squares MPM (Hu et al. 2018)

MLS-MPM further simplifies APIC by using the Moving Least Squares framework to unify the P2G and G2P transfers. The result is a formulation where the entire MPM step can be written in ≈60 lines of code with no separate stress divergence computation — the Cauchy stress enters naturally through the affine matrix update. MLS-MPM is the basis of the popular taichi simulation framework and is what most modern MPM tutorials implement.

// MLS-MPM P2G (full): stress enters as part of affine momentum // D_p = (1/4) h² I (for quadratic B-splines) momentum_i += w_ip * (m_p * v_p + (m_p * C_p − Δt * V_p * σ_p D_p⁻¹)(x_i−x_p)) // ↑ stress term folded into affine momentum // No separate force loop! The stress divergence is built into the affine splat. // Entire 2D MPM engine: ~50 lines of code

7. Implementation Sketch in JavaScript

A minimal 2D MLS-MPM elastic solid can be implemented in under 100 lines. The key data structures are flat typed arrays for performance:

// 2D MLS-MPM core data structures const N_GRID = 64; // grid resolution (64×64) const DX = 1.0 / N_GRID; const INV_DX = N_GRID; const N_P = 8000; // particle count // Particle state (SoA layout for cache efficiency) const px = new Float32Array(N_P); // position x const py = new Float32Array(N_P); // position y const vx = new Float32Array(N_P); // velocity x const vy = new Float32Array(N_P); // velocity y const Cxx = new Float32Array(N_P); // affine matrix C (4 components) const Cxy = new Float32Array(N_P); const Cyx = new Float32Array(N_P); const Cyy = new Float32Array(N_P); const Jdet = new Float32Array(N_P).fill(1); // volume ratio (J = det F) // Grid (reset each step) const gm = new Float32Array((N_GRID+1)**2); // grid mass const gvx = new Float32Array((N_GRID+1)**2); // grid momentum/velocity x const gvy = new Float32Array((N_GRID+1)**2); // grid momentum/velocity y function step(dt) { gm.fill(0); gvx.fill(0); gvy.fill(0); // P2G reset for (let p = 0; p < N_P; p++) { // P2G const bx = Math.floor(px[p] * INV_DX - 0.5); const by = Math.floor(py[p] * INV_DX - 0.5); // ... quadratic B-spline weights ... // ... splat mass + momentum + stress ... } for (let i = 0; i < (N_GRID+1)**2; i++) { // Grid update if (gm[i] > 0) { gvx[i] = gvx[i]/gm[i]; gvy[i] = gvy[i]/gm[i] + dt * -9.8; // gravity } // boundary conditions: zero velocity at grid edges } for (let p = 0; p < N_P; p++) { // G2P // gather velocity and C matrix from grid // update J: J *= 1 + dt*(dv_x/dx + dv_y/dy) // advect: px[p] += dt*vx[p]; py[p] += dt*vy[p] } }
Performance: A 2D MLS-MPM simulation with 8 000 particles runs at ~120 fps in JavaScript. A 3D simulation with 100 000 particles requires a WebGL or WebGPU compute shader implementation to run in real-time — the P2G/G2P loops are embarrassingly parallel over particles, and the grid update is parallel over grid nodes.

8. Applications: Disney, Games, Engineering

Film VFX

Games (Real-Time MPM)

Recent advances in GPU-accelerated MPM allow real-time use in games. The key optimisation is running the grid at coarser resolution (e.g. 32³ cells) with smaller particle counts (10 000–50 000), using WebGPU or Vulkan compute shaders for the particle-grid transfers.

Scientific and Engineering Applications

MPM continues to be one of the most active research areas in physics simulation — new constitutional models appear regularly, and the convergence of MPM with machine learning (learned constitutive models, neural correction terms) is an active frontier in both graphics and engineering simulation.