Most MPM material is written in 3D, and 2D is usually introduced as the easy case you build first. That framing is misleading in one specific way, and the mistake it invites is expensive: a 2D MPM solver is not a 3D solver with a thin third axis. It is a different continuum assumption, it obeys different arithmetic, and it costs an order of magnitude less. This document sets out why you would want the 2D path, what separates it from 3D, what a 3D solver can and cannot be made to reproduce, and a recipe for building 2D if you decide you want it.
It assumes the deformation gradient, the transfer cycle, and B-spline weighting from Prerequisites and Theory.
1. Why build a 2D solver at all
Five reasons hold up. One common one does not.
Cost, and not by a small factor. The grid scales with the square of resolution rather than the cube, and the transfer stencil touches nine nodes per particle rather than twenty-seven. At 256 cells across the play area a 2D grid is roughly 66 thousand nodes; a 3D grid at the same in-plane fidelity is 17 million (§5). That is not an optimization, it is a change in which hardware can run the simulation at all — a 2D continuum solver at good resolution fits comfortably on integrated graphics and mobile parts, where the 3D equivalent does not fit at any frame rate.
Games whose play space really is a plane. In a side-on or top-down game the simulation plane is the play plane, so nothing is spent simulating depth the player can neither see nor act in. It is also where MPM’s advantage over the alternatives is widest. The cheap incumbents for 2D destructible matter — cellular automata, height fields, sprite particles — stay convincing until the material has to hold a slope, shear, pack under its own weight, or move between solid and flowing behavior. That is exactly what a continuum solver gives you and what they structurally cannot; Games covers the incumbents and where they run out. A 2D solver can ship that behavior on hardware where 3D MPM is still out of reach.
Developing a constitutive law. A material law is the same law in either dimension, and the plane is far cheaper to explore. A parameter sweep that takes a night in 3D takes minutes in 2D, and a wrong yield surface or a sign error in the stress is usually visible in the first frame rather than after a settle. The compact reference implementations most people learn the method from are 2D for this reason.
Debuggability. In 2D the entire grid fits on screen at once, every particle can be drawn carrying its own deformation state, and the state is small enough to hash and diff step by step against a reference run. The same inspection in 3D needs slices, cutaways, and sampling — you are working through a keyhole, and the bug is often in the part you did not slice.
Comparison against published results. Much of the MPM benchmark literature — granular column collapse, dam break, the classic elastic and plastic bar tests — is reported in plane strain. Reproducing it is the cheapest real evidence that a solver is correct at all, and it requires plane strain specifically, not a thin slab (§3).
What it is not: practice for 3D. The transfer cycle carries over and the constitutive laws carry over. Rotation extraction does not (§6.4), the grid accumulator layout does not, and no calibrated material constant does, because the volume ratio driving the law means something different in the plane (§3). Build 2D because you want a 2D result, not as a warm-up.
2. 2D is a confinement assumption
When the third axis is not simulated, its absence is not neutral. It encodes a boundary condition: strain along that axis is identically zero. This is plane strain, and the material behaves as though it were infinitely thick out of the plane — a prism extending without limit, seen end-on.
Given identical material parameters and identical loading, the two configurations do not merely differ in cost — they settle to different shapes.
3. Why a thin slab does not converge to it
The intuition that a slab approaches the 2D case as it gets thinner is exactly backward. Thinning it introduces free surfaces on both faces and moves the material toward plane stress, the opposite limit, in which out-of-plane stress rather than strain is zero. The two limits bracket the 3D case; neither is reachable from the other by refinement.
Three further consequences are worth knowing before choosing a slab:
- The only way to recover plane strain is to impose it. Zero the out-of-plane velocity component and the corresponding row and column of both the affine velocity matrix and the deformation gradient, every step. That works — and at that point the determinant driving the constitutive law is a 2×2 determinant and you have written a 2D solver inside a 3D one, paying 3D cost for 2D physics.
- The volume ratio changes meaning. The determinant of the deformation gradient is an area ratio in the plane and a volume ratio in space. Every constant calibrated against one — bulk modulus and the exponent in a Tait-style equation of state, hardening coefficients, yield surfaces — needs recalibration for the other.
- A slab has a floor thickness. The quadratic B-spline spans three cells per axis, so a slab thinner than three cells has no interior at all: every particle sits in the boundary layer, and grid boundary conditions switch on and off against the stencil rather than bounding a well-posed interior.
4. What a 3D solver can and cannot replicate
Much of what reads as “the 2D look” is reachable from a 3D solver. The physics is not.
| Property of a 2D simulation | Reachable from a 3D solver? |
|---|---|
| Flat, screen-aligned silhouette | Yes — orthographic camera plus an unlit shading path |
| Reconstruction field finer than the solver grid | Yes, but the field grows with the cube |
| Temporal stabilizer on the density field | Only by paying for a second volume |
| Grid cost that scales with the square of resolution | No — a 3D grid allocates the full cube |
| Zero out-of-plane strain (plane strain) | No — only by constraining the third axis every step, which is a 2D solver in disguise |
| Area-based volume ratio in the constitutive law | No — it follows from the same constraint |
Read the table as a decision aid. If the reasons you want a 2D simulation are all in the first three rows, fix the 3D renderer instead. If they are in the last three, build the 2D path.
5. The cost that usually decides it
Node count is the dominant term, and it is the square versus the cube of the resolution. Assuming a modest per-node footprint — a small fixed-point accumulator plus a velocity vector, around 32 bytes — a single-field grid costs:
| Cells across the play area | 2D grid nodes | 3D grid nodes | 3D grid memory |
|---|---|---|---|
| 96 | 9.4k | 913k | ~28 MB |
| 128 | 17k | 2.1M | ~66 MB |
| 192 | 37k | 7.2M | ~219 MB |
| 256 | 66k | 17.0M | ~518 MB |
Multi-field contact, which keeps several velocity fields per node, multiplies the 3D column by the slot count.
Two traps hide inside those numbers:
- Allocation is cubic even when occupancy is not. A grid indexed over a cube is normally allocated over the whole cube, so a slab occupying three percent of the domain still reserves all of it. Sparse-grid indexing thins the dispatch — and is the right answer for a slab, whose occupancy is genuinely low — but it does not by itself shrink a densely allocated accumulator. Check which of the two a given implementation actually does before budgeting.
- Particle count scales too. For the same in-plane particle spacing, a slab needs one layer of particles per cell of thickness. A slab eight cells thick is roughly eight times the particles as well as the larger grid.
Against that, the 2D stencil is nine nodes per particle rather than twenty-seven, so scatter traffic — the usual bottleneck, see Implementation §1 — falls by a further factor of three per particle.
6. Building a 2D solver
6.1 The particle row
Position, velocity, the affine velocity matrix, and the deformation gradient, plus whatever scalars the constitutive law needs: fifteen floats for an elastoplastic material with one plastic scalar. Pad to a power-of-two stride. Put position at offset zero so a renderer can read positions without knowing the rest of the row — that one convention is what lets the rendering layer depend on a published row contract rather than on solver internals.
6.2 Grid and stencil
Nodes number the square of resolution plus one per axis. The quadratic B-spline gives three weights per axis, so nine nodes per particle. Keep the weight computation identical to the 3D form; only the loop nesting changes.
6.3 The transfer cycle
Clear, particle-to-grid, grid update, grid-to-particle — structurally identical to 3D. Use fixed-point integer atomics for the scatter rather than float atomics: integer atomics are the portable choice across graphics backends, and they make the scatter deterministic, which float atomics are not. Accumulate momentum components and mass per node, and choose the fixed-point scale so the worst-case accumulated momentum stays inside the signed 32-bit range.
6.4 Rotation extraction — the one place the code is genuinely different
Everywhere else, 2D code is 3D code with a narrower loop. Here it is a different algorithm, and the simplification is large.
The polar decomposition of a 2×2 matrix is closed form. For a matrix with entries a, b on the first row and c, d on the second, the best-fit rotation has angle ; normalizing the vector gives its cosine and sine directly, and the rotation follows in about a dozen lines with no iteration. A signed 2×2 singular value decomposition is then that rotation followed by a symmetric eigendecomposition of the remaining stretch — call it twenty-five lines in total.
The 3×3 equivalent is an iterative Jacobi routine, typically well over a hundred and fifty lines, with its own convergence and sign-handling concerns. Do not attempt to share the two behind a common interface; write the closed form and move on.
6.5 Boundary conditions
Identical to 3D with one axis removed. The distinction that matters is the same one: separating boundaries kill only the into-boundary velocity component, sticky boundaries kill the whole vector, and granular material wants the former on every wall.
6.6 What to share with a 3D solver, and what not to
This decision determines the long-run cost, and the intuitive answer is wrong.
Share the stable parts. Constitutive laws are dimension-parameterized formulas — fixed corotated elasticity, Drucker-Prager, snow-style hardening, a Tait equation of state — and one source can serve both dimensions if it is written against a matrix type with determinant, transpose, and polar-decomposition overloads. Material parameter definitions, and the buffer lifecycle around a solver (upload, seeding, disposal), share cleanly too.
Do not share the transfer kernels. It is technically easy — a dimension macro and a conditional loop nest — and it is a trap whenever one of the two solvers is under active development and the other is not. Duplication between an active path and a dormant one is precisely what makes the dormant one free to ignore. A shared kernel makes every change to the active path a change to the dormant one by construction, so every edit must then be verified in both dimensions whether or not anyone wanted the second dimension that day. An ignorable cost becomes an unavoidable one. Duplicate what is changing; share what has settled.
7. Rendering a 2D simulation
Both reconstruction paths start from the same particle buffer and diverge immediately:
The structure is the same — clear, splat, composite — and the reconstruction rules from Rendering apply unchanged: match the splat kernel to the solver’s transfer kernel, and let density decide in-or-out only, never opacity.
What differs is affordability, and it decides one design point in particular. A temporal stabilizer costs a second copy of the field. In the plane that is nothing: at twice the solver resolution a 256-square field is about a quarter of a megabyte, so keeping a previous frame to blend against is free, and an exponential moving average over the density field is the cheapest available cure for boundary shimmer in sparse regions. In 3D the same field is 64 MB and the history is 64 MB more, so the blend tends to be built as opt-in and left off, with its history volume allocated only when someone switches it on; a spatial blur over the composited field is the cheaper habit and the one that ships enabled.
The stabilizer is therefore available in both — what changes is whether you can leave it on without thinking about it. That is the general shape of the 2D advantage in reconstruction: not that a technique is unavailable in 3D, but that the cube makes you ration it.
Two smaller notes for anyone driving a 3D renderer at a planar scene:
- Raymarching is usually perspective-only. A fullscreen reconstruction pass that builds each pixel’s ray from a single camera position is correct only under perspective projection. Under an orthographic camera every ray shares a direction and the origins vary across the image plane, so the origin must come from unprojecting the pixel at the near plane instead. This is a few lines and is correct for both projections; it is also easy to miss, because the perspective form looks right until the day someone switches the camera.
- Gradient-normal shading is what makes a slab read as three-dimensional. If the flat silhouette is the goal, an unlit path that keeps the threshold and discards the normal gets most of the way there.
8. Choosing
| If you need… | Choose |
|---|---|
| Plane-strain physics, or comparison against 2D benchmarks in the literature | 2D |
| Large in-plane resolution — more than roughly 128 cells across | 2D; the cube is decisive |
| To ship on integrated graphics or mobile | 2D |
| A temporal stabilizer on the reconstructed field | 2D, unless you can spare the second volume |
| A flat silhouette and nothing else from the list above | 3D, with an orthographic camera and unlit shading |
| Cutting, multi-field contact, fracture, or any out-of-plane flow | 3D |
| Small domains where the cube is affordable | 3D; one solver is cheaper to maintain than two |
The last row is the one most often underweighted. Two solvers is not twice the code, it is twice the verification: every feature needs a design, an implementation, a gate, and a regression baseline in each dimension. Build the 2D path when a row above it applies — not because it looks like the easier place to start.
Sources for this document
The continuum-mechanics background — plane strain and plane stress, the deformation gradient, polar and singular value decomposition — is in Prerequisites. The transfer cycle and B-spline weighting are in Theory. The scatter bottleneck and fixed-point atomics are treated at length in Implementation §1–2, the reconstruction method space in Rendering §1, and the 2D incumbents a continuum solver has to beat in Games §4.1. Papers for all of it are in References.