Material Point Method (MPM) — How to Simulate Snow, Sand, and Soft Bodies

1. Why MPM? Lagrangian vs Eulerian

Continuum mechanics simulations split into two families:

• Lagrangian (particle/mesh): elements follow the material. Great for tracking large deformation and fracture. Used in FEM for solid mechanics. Breaks down when mesh becomes too distorted (tearing, extreme flow).
• Eulerian (grid): grid is fixed; material flows through it. Great for fluid — no mesh tangling. Hard to track sharp material interfaces or history-dependent plasticity.

MPM sits in between: particles carry the material history (deformation gradient, stress, velocity, density) while a fresh background grid is recreated every time step to handle contact, pressure projection, and boundary conditions. The grid is discarded after each step — no tangling possible.

2. Algorithm Overview: Four Steps per Frame

Each simulation step performs four key operations:

1. P2G (Particle → Grid): splat particle mass and momentum onto nearby grid nodes using interpolation weights
2. Grid forces: compute and apply forces (elastic stress, gravity, collision impulses) to grid node velocities
3. Boundary conditions: enforce walls, sticky/slippery surfaces, objects on the grid
4. G2P (Grid → Particle): update particle velocities and positions from the grid; update deformation gradient

3. P2G — Particle-to-Grid Transfer

Each particle affects the nearest 3×3 (2D) or 3×3×3 (3D) grid nodes. The interpolation kernel is a quadratic B-spline — chosen for its compact support (2 cells), C¹ continuity, and partition of unity property:

The elastic stress contribution to the grid force comes from the First Piola-Kirchhoff stress P = ∂ψ/∂F, where ψ is the strain energy density and F is the deformation gradient:

4. Grid Update — Forces and Boundary Conditions

Once mass and momentum are rasterized onto the grid, forces are integrated explicitly:

Imposing BCs on the grid is trivial compared to Lagrangian FEM where contact requires complex gap-function formulations. This is one of MPM's key practical advantages.

5. G2P — Grid-to-Particle Transfer

After the grid is updated, information flows back to particles. The critical operation is updating the deformation gradient F, which encodes how the material has stretched and sheared since the start of the simulation:

Particle positions are updated last, after the plastic return mapping has clamped F to the yield surface (for snow/sand models):

6. Constitutive Models

The deformation gradient F drives the stress computation. Different materials use different strain energy densities:

Neo-Hookean Elastic Solid

Snow (Disney 2013)

Snow fractures under tension and compacts under compression. This is modeled via multiplicative plasticity:

Drucker-Prager Sand

Sand has a friction-angle yield criterion — it can sustain compression but not tension, and shear strength scales with pressure:

Weakly Compressible Fluid (SPH-like pressure)

7. APIC and FLIP Variants

The basic PIC scheme (reconstruct velocity purely from grid) is numerically diffusive — angular momentum is lost. Two improvements exist:

• FLIP (1988): only transfer the velocity increment from the grid: v_p += (v_grid_new − v_grid_old). Preserves more energy, but accumulates noise over time.
• APIC (Jiang et al. 2015): Affine Particle-In-Cell. Each particle stores a local velocity field approximation matrix C_p. P2G uses v_p + C_p·(x_i − x_p); G2P reconstructs C_p from the grid. Conserves angular momentum exactly with zero numerical dissipation — the current state of the art for most MPM applications.

8. MLS-MPM Core Loop in JavaScript

The following shows the essential 2-D MLS-MPM loop — streamlined for clarity. A full runnable implementation with canvas rendering needs ~200 lines:

// MLS-MPM core: 64×64 grid, quadratic B-spline weights
const n = 64; // grid resolution
const dt = 1e-4;
const dx = 1 / n;
const inv_dx = n;

// Particle state
let pos, vel, C, J; // position, velocity, affine matrix, volume ratio

function mpmStep() {
  // 1. Reset grid
  const gm  = new Float32Array(n * n);      // mass
  const gmv = new Float32Array(n * n * 2); // momentum xy

  // 2. P2G
  for (let p = 0; p < pos.length; p++) {
    const [px, py] = pos[p];
    const base = [Math.floor(px * inv_dx - .5), Math.floor(py * inv_dx - .5)];

    // Quadratic B-spline weights  (3 nodes per axis)
    const fx = [px * inv_dx - base[0], px * inv_dx - base[0] - 1, px * inv_dx - base[0] - 2];
    const fy = [py * inv_dx - base[1], py * inv_dx - base[1] - 1, py * inv_dx - base[1] - 2];
    const wx = [.5 * (1.5-fx[0])**2, .75-fx[1]**2, .5*(.5+fx[2])**2];
    const wy = [.5 * (1.5-fy[0])**2, .75-fy[1]**2, .5*(.5+fy[2])**2];

    // MLS-MPM stress term (Neo-Hookean fluid: P = k₀(1 - 1/J) I)
    const stress = -dt * 4 * inv_dx * inv_dx * 400 * (1 - 1 / J[p]);

    for (let di = 0; di < 3; di++)
      for (let dj = 0; dj < 3; dj++) {
        const w = wx[di] * wy[dj];
        const ix = base[0] + di, iy = base[1] + dj;
        if (ix < 0 || ix >= n || iy < 0 || iy >= n) continue;
        const idx = ix * n + iy;
        const dpos = [(di - fx[di]) * dx, (dj - fy[dj]) * dx];
        gm[idx]        += w;
        gmv[idx*2]   += w * (vel[p][0] + C[p][0]*dpos[0] + stress * dpos[0]);
        gmv[idx*2+1] += w * (vel[p][1] + C[p][3]*dpos[1] + stress * dpos[1]);
      }
  }
  // 3. Grid update: normalize by mass, add gravity, enforce boundary
  // 4. G2P: reconstruct velocity and update F, x  (loop over particles)
}