The Material Point Method (MPM) is a powerful and interesting
technique for physics simulation. MPM is great for continuum substances
like fluids and soft bodies, as
well as for dealing with complex problems like fracture and self
collisions. It has received a lot of recent attention from Disney, most
notably the

snow simulation and

phase change/heat transfer papers that appeared in SIGGRAPH.

Recently I coded up optimization based implicit MPM to test out some research ideas. My ideas didn't work out, but I put the

code up on github instead of letting it collect dust on my hard drive. There are many graphics-oriented MPM implementations online, but very few of them use implicit integration, and even less (any?) use optimization.

I'm
going to skip a lot of the mathematical details in this post. I feel
they are better covered by the referenced papers. Instead, I'll give just enough to understand the code.

##
Material Point Method

The basic idea is that particles move
around the domain and carry mass. When it's time to do time integration,
the particle mass/momentum is mapped to a grid. New velocities are
computed on the grid, then mapped back to the particles, which are then
advected. The grid essentially acts as a scratch pad for simulation that
is cleared and reset each time step. For a (better) background on MPM,
read the introduction chapter from the

SIGGRAPH course notes.

###
Time integration

The best starting point for computing new velocities is explicit time integration.
Though explicit integration is fast to compute and easy to understand,
it has problems with stability and overshooting. If we take too large of
a time step, the particle goes past the physically-correct trajectory.
Intermediate techniques like midpoint method can alleviate such
problems, but the silver bullet is implicit integration.

With implicit integration we figure out what the velocity needs to be to make the particles to
end up where they should. The hard part is that we now have both
unknowns in the velocity calculation, which requires solving a system of
nonlinear equations. This is computationally expensive and time
consuming, especially for large and complex systems.

Fortunately, many
important forces (like elasticity) are conservative. Mathematically,
that means we can represent the forces as the negative gradient of
potential energy. This lets us do implicit integration through a technique called

*optimization*. This makes it a bit easier to formulate and can have better performance than a lot of other implicit solvers.

###

By representing our forces as energy
potentials, we can do implicit time integration by minimizing total
system energy. The concept itself is not new, but was recently applied
to Finite Elements/MPM in an

excellent paper by Gast et al. Section 6 of the paper covers the algorithm for optimization-based MPM. One thing to keep in mind: in MPM the velocity calculations happen on grid nodes, not grid boundaries.

###
PIC/FLIP versus APIC

Mapping to/from the grid is an important part of the material point method. MPM was originally an extension of Fluid Implicit Particle (FLIP), which in turn was an extension of the Particle in Cell (PIC) method. PIC and FLIP are the standard in fluid simulation, and were used in the above Gast paper. Mike Seymour wrote a

good article on fxguide that covers what PIC/FLIP is and why they are useful. However, they suffer from a few well known limitations. Namely, PIC causes things to "clump" together, and FLIP can be unstable.

A technique called the

Affine Particle In Cell (APIC) was recently introduced to combat these limitations. It's more stable and does a better job at conserving momentum. In short, it's worth using. To do so, you only need to make a few subtle changes to the MPM algorithm. The SIGGRAPH course notes explain what to change, but the Gast paper does not.

##
Implementation

Everything you need to know to code up MPM can be found in the SIGGRAPH course notes. If you've ever done fluid simulation with PIC/FLIP, it should be straight forward. If not, here are some things I found useful:

1. It helps with stability to have a cell size large enough to have 8-16 particles.

2. Keep a fully-allocated grid as the domain boundary, but do minimization on a separate "active grid" list made up of pointers to grid nodes.

3. A "make grid loop" function for particle/grid mapping: given a particle, return the indices of grid nodes within the kernel radius, as well as kernel smoothing values/derivatives. This was probably the trickiest part of the whole method. Much of this can be preallocated if you have a lot of memory, but I found (in 3D) allocating/deallocating on the spot was necessary.

Another thing to consider is what algorithm to
use for minimization. I'm too lazy to deal with Hessians, so I generally
stick to L-BFGS and Conjugate Gradient. Optimization is widely used across disciplines, so there are some open-implementations you may want to
consider using. I've found that

Pat Wie's CppOptimizationLibrary works well as a go-to for most of my test projects, because I like the interface and the Eigen3 compatibility. It also includes a lot of nice things like finite difference gradient/hessian, line search algorithms, etc... However, there are some issues with its L-BFGS and CG. For some of my applications, they failed to converge (although that may have been fixed in recent revisions, I haven't checked). Instead I use an L-BFGS implementation by

Ioannis Karamouzas that I modified for Pat Wie's library (included in my MPM github repo).

##
Rendering

Since the code was just to test out some ideas, I didn't put any effort into rendering. Instead, I just visualized the particle locations using

VDB. If you haven't used VDB, it's a pretty cool tool for drawing out primitives in fresh simulation code. All you have to do is run the viewer, include the header, then call a few functions (like "vdb_point(x,y,z)") to draw stuff.