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).
