information on MLS-MPM material point method from nialltl.neocities.org/articles/mpm_guide https://github.com/nialltl/incremental_mpm
MPM is a modern general approach to the simulation of substances like fluids, gasses, and deformable objects, where each point can be described by consistent equations. This page describes MLS-MPM which is simpler. # References and terms first. Then notes. text pasted from page. MPM Course (SIGGRAPH 2016 notes): http://mpm.graphics absolute behemoth of a walkthrough. probably the best bottom-up explanation of MPM i’ve found! i’ll reference this throughout my examples. APIC original paper: https://www.math.ucla.edu/~jteran/papers/JSSTS15.pdf a useful reference for the actual implementation of grid-particle transfers that are used in MLS-MPM. Taichi MPM implementation, slides & notes https://github.com/yuanming-hu/taichi_mpm/ the PDF slides and slides with notes are really helpful and give a great side-by-side comparison of MPM and MLS-MPM. MLS-MPM original paper: https://www.yzhu.io/projects/siggraph18_mlsmpm_cpic/index.html this might be a bit overwhelming to launch straight into without implementing standard MPM first, but i’d recommend the section “from MPM to MLS-MPM”. MLS-MPM implemented in 88 lines of c++ https://github.com/yuanming-hu/taichi_mpm/blob/master/mls-mpm88-explained.cpp this is the comments-expanded version of Yuanming Hu’s original 88-line implementation. similar in structure to our Part 2 implementation, but with a Disney-style constitutive model for snow. useful to see how they differ! a maintained curated list of current MPM research https://www.seas.upenn.edu/~cffjiang/mpm.html ## appendix – simulation methods and their general vibes Lagrangian SPH: Smoothed Particle Hydrodynamics particle-based fluid simulation method. computes field quantities by looking at surrounding particle neighbourhoods. has a hilarious amount of variations differing mostly in acronym length (ISPH, IISPH, PCISPH, DCSPH, … generally trying to improve incompressibility). PBD: Position Based Dynamics general simulation methodology, handles constraints geometrically using direct positional corrections. not that physically motivated but generally fast and great for real-time applications. super stable. often suffers from damping, but higher order PBD integrators help with this. PBD: Position Based Fluids like SPH, but using a PBD integrator. enforces incompressibility by putting a constraint on particles’ densities. Hybrid Lagrangian-Eulerian PIC: Particle-In-Cell general family of methods using a mixture of grids and particles for simulation. an absolute classic, been around a long time. used for astronomical simulations and liquids / gases generally. FLIP: Fluid Implicit Particle like PIC, but with an improved way of advecting particles. typically used for fluids. common in a lot of 3D simulation packages, e.g. Blender, Houdini. MPM: Material Point Method another PIC method, more versatile than something like FLIP. good for deformable soft matter and fluids. supports multiple phases (different materials interacting with each other) fairly effortlessly. typically requires a very small timestep, comparatively high memory requirements. easily exhibits numerical fracture (stuff splitting apart), which you might not always want. GIMP: Generalised Interpolation Material Point Method a particular flavour of MPM that emerged early on that solved early issues with linearly weighted MPM, to stop particles oscillating erroneously at cell boundaries. APIC: Affine Particle-In-Cell very recent significant improvement on PIC, one-upping a fix for PIC’s dissipation problems that FLIP tried to solve. PolyPIC: Polynomial Particle-In-Cell the wider family of methods that APIC belongs to. basically a higher-order generalisation. MLS-MPM: Moving Least Squares Material Point Method like MPM, but taking ideas from APIC and bringing them in state o’ the art. a straight up improvement on MPM. simplifies the whole algorithm and is more computationally efficient too. # To describe the domain people considered different ways of representing the material -- e.g. whether to discretize particles or space. But the approaches that seemed most effective used different representations for different parts of the algorithm. These are called Hybrid. The page describes an earlier hybrid approach called PIC and shows a demo animation. PIC made fluids seem too rubbery and viscous ("mushy"). The page says this is because of numerical accuracy lost during a back-and-forth transformation in the step update. The next section displays an example of FLIP, which modifies PIC by something like preserving a hardcoded percentage of velocity in an area. The FLIP simulation looks much more watery with little droplets and splashes going around. The stated problem with FLIP is that it is not very stable, the particles and surfaces jitter and simulations contain noise. The article states that these problems of FLIP have been resolved in APIC and PolyPIC which associate a context information matrix with each particle. To clarify, it sounds like the problems these algorithms encounter relate to holding both a grid-based and particle-based set of information at the same time, and how that causes drift. (sounds like a potentially fun problem) MPM is a generalization of PIC/FLIP . MLS-MPM adds an approach similar to APIC which also simplifies the implementation by removing grid-based gradients. The author says to reference the SIGGRAPH 2016 course for further information. ## The page says an MLM "hello world" will be made, a 2D simulation of an elastic material. Fun big words in the paragraph: "i’m going to be building a 2D simulation of an elastic material, using MLS-MPM with quadratic b-spline weights under an explicit time integration scheme." The author chooses "Unity's data-oriented language, High Performance C#" (strange sentence construction, gives one small curiosity as to this language). The author says there isn't really a good reference for MLS-MPM aside from the paper (at time of writing). So I should take notes on it! Near minumum: - The simulation has both a particle set and a grid - Each particle has position, velocity, and mass (3 fields) - Each grid cell has velocity and mass (2 fields) (grid attributes are a subset of particle attributes) The convention is to have the particle positions be scaled to grid indices, so that a particle at (1,1) is at [1,1], and to use a power of 2 for the total simulation size. The grid is a temporary per-frame construct used for intermediate data. Realistic simulation: - Each particle also has a deformation gradient matrix, and information regarding its volume. - For MLS-MPM, each particle also needs an APIC-style affine momentum matrix. - In traditional MPM, each grid cell needs a force vector. In MLS-MPM, regarding the force vector, you can "instead fuse it together with the velocity update". ## The author says for a slower walkthrough, refer first to Chapter 10.5, “MPM Scheme: Full Algorithm” in the MPM course for the standard MLM procedure, and then Section 4, “From MPM to MLS-MPM” in the original MLS-MPM paper. Pseudocodeish: initialization: construct particles and a grid, each particle has an identity matrix for its deformation update: - clear grid, mass=0 velocity=0 - particle-to-grid. for each particle: - use an interpolation function to calculate weights for all neighboring cells - calculate quantities such as stress based on the substance equations - apply the particle's momentum to its cell and every neighboring cell based on the interpolation weights - grid update: calculate velocity in each grid cell based on its momentum. enforce boundary conditions. - grid-to-particle: for each particle: - update particle's deformation gradient using MLS-MPM paper equation 17, the velocity gradient estimate - recalculate the interpolated cell weights if needed - use the weights to set the particle's velocity from all neighboring cells - update the particle positions There ya have it. Sounds simple and intuitive aside from 2-3 missing heuristic/equations. The author says the source code contains inline references to parts of papers and that reviewing these will make the approach quite clear. The hello world demo is at https://github.com/nialltl/incremental_mpm . It sounds likely that the particles do not enforce volume and all fall to the bottom of the grid because there is no deformation implemented yet. It looks like the code is in the Assets subfolder. I can now infer here that "quadratic b-spline weights" is likely the use of a quadratic function for neighbor interpolation, whereas time integration might relate to the position update. ## I'm now looking at the C# source code for this first demo. A single scalar called "padding" has been added to each particle and each grid cell, commented as being for performance. Note on time steps from code: for the simulation to be stable, a particle should travel no more than 1 grid cell per update. The author hardcodes a constant small timestep. (opportunities here) Here's their quadratic interpolation in particle-to-grid: // quadratic interpolation weights uint2 cell_idx = (uint2)p.x; float2 cell_diff = (p.x - cell_idx) - 0.5f; weights[0] = 0.5f * math.pow(0.5f - cell_diff, 2); weights[1] = 0.75f - math.pow(cell_diff, 2); weights[2] = 0.5f * math.pow(0.5f + cell_diff, 2); They state higher order approaches with larger "stencils" yield more accuracy. I've certainly seen this formula before but I don't recall what it is. It's nonintuitive that a particle's mass and velocity would disperse onto a grid using a blur function. But that appears to be what is happening here. The author says this works much better than linear interpolation, which makes sense. The final weight appears to be the product of the weights indexed for 2 dimensions (weights[x].x * weights[y].y). This could make sense when one considers that each axis has the same fixed sum. - A vector Q is set to the product between the particle's momentum matrix C and its distance from the center of each cell it contributes to. - The interpolation weights scale the mass - Each cell's mass (initialized to 0) is increased by the partial mass. - The velocity is summed with Q. Then each cell's velocity (also initialized to 0) is increased by the product of this sum with the partial mass. So there are two unmentioned things going on there. The simpler one is that the cell velocity is weighted by the partial masses of the particles contributing to it. Since it also has the sum of these partial masses, an average velocity could be found via division. The more complex thing is that the velocity transferred to the grid is summed with a matrix product based on distance from the cell. This matrix is always 0 in this intro code, so if there were a typo here it wouldn't have been discovered. But we can literally infer I suppose that this particle matrix is how to calculate a velocity contribution to each grid cell based on its distance from the particle. Note that momentum (mass * velocity) is what is actually being stored in the velocity field of the grid, and what this matrix is used for. In the grid update, this velocity field is converted to velocity from momentum by dividing by mass. Gravity is calculated via += dt * g. Boundary conditions are enforced by setting axial velocity at the grid edges. (developing understanding of this algorithm shows how the need for grid boundaries might be removed via discretization of the data at a fixed resolution around every particle.) Here's the comment for grid-to-particle: // constructing affine per-particle momentum matrix from APIC / MLS-MPM. // see APIC paper (https://web.archive.org/web/20190427165435/https://www.math.ucla.edu/~jteran/papers/JSSTS15.pdf), page 6 // below equation 11 for clarification. this is calculating C = B * (D^-1) for APIC equation 8, // where B is calculated in the inner loop at (D^-1) = 4 is a constant when using quadratic interpolation functions float2x2 B = 0; Particle velocity is initialized to 0, and then grid velocities summed in via the weights. Here's how B is updated: // APIC paper equation 10, constructing inner term for B var term = math.float2x2(weighted_velocity * dist.x, weighted_velocity * dist.y); B += term; dist is again the distance to the center of the interpolated cell (inside 2 nested loops). So B ends up being a sum of these distances * velocities. ie, B contains the velocities scaled by their distances from the cells. Noting that the upscaled velocities have small components from the interpolation function -- ie the weighting and distance scaling may partly counter each other. p.C = B * 4; What a random constant "4"! can make sense since there are 4 values. But this intro code doesn't use C so hasn't explained it yet. p.x += p.v * dt; Same kind of update as gravity. I don't usually call this integration but maybe that's right. Discrete. It's discrete integration I guess. The code also includes this, although the boundary conditions earlier should make this unneeded I think? // safety clamp to ensure particles don't exit simulation domain p.x = math.clamp(p.x, 1, grid_res - 2); that's it! The full sourcefile is at https://raw.githubusercontent.com/nialltl/incremental_mpm/refs/heads/master/Assets/1. MLS_MPM_Intro_SingleThreaded/MLS_MPM_Intro_SingleThreaded.cs ## Example 2 is isotropic elasticity. The animation looks roughly like a dropped blob of gelatin --it bounces and breaks a little. Unlike the previous example where all the particles collected at the bottom, this object retains volume and structure when dropped (as do the chunks that break off of it). In theory this is may just be a small change from example 1! Pseudocodeish, for each particle in particle-to-grid: - each particle has F a deformation gradient, and volume_0 - 'MPM course page 13, "Kinematics Theory"' J is taken as the determinant of F, and scales volume_0 for particle-to-grid - 'Neo-Hookean model (eq. 48, MPM course)' P is calculated such that P = elastic_mu * (F - inverse(transpose(f))) + elastic_lambda * log(J) * inverse(transpose(F)) elastic_mu scales the F self-difference value, whereas elastic_lambda scales the F inverse value with the log of its determinant - 'Cauchy stress = (1 / det(F)) * P * F_T; equation 38, MPM course' stress is set as P * transpose(F) / J - '(M_p)^-1 = 4, see APIC paper and MPM course page 42' 'this term is used in MLS-MPM paper eq. 16. with quadratic weights, Mp = (1/4) * (delta_x)^2.' 'in this simulation, delta_x = 1, because i scale the rendering of the domain rather than the domain itself.' 'we multiply by dt as part of the process of fusing the momentum and force update for MLS-MPM' eq_16_term_0 = -volume * 4 * dt * stress; note volume is volume_0 * J - for each interpolated cell: 'we calculate the rest of the force term inside the loop' 'fused force/momentum update from MLS-MPM' 'see MLS-MPM paper, equation listed after eqn. 28' float2 momentum = math.mul(eq_16_term_0, cell_dist) * weight; cell.v += momentum; 'total update on cell.v is now: weight * (dt * M^-1 * p.volume * p.stress + p.mass * p.C)' 'this is the fused momentum + force from MLS-MPM. however, instead of our stress' 'being derived from the energy density, i use the weak form with cauchy stress. converted:' 'p.volume_0 * (dΨ/dF)(Fp)*(Fp_transposed)' 'is equal to p.volume * σ' so eq_16_term_0 scales the distance contribution to momentum this doesn't sound that different from the source code i reviewed already but [it's quite confusing, referencing many new equations [.... might take a break unsure. the comment in the page says that Position Based Fluids are actually better at fluid sims, MPM takes tweaking. A web search shows there is now something called Position Based MPM GitHub - electronicarts/pbmpm: A WebGPU implementation of Position Based MPM (PB-MPM), presented at SIGGRAPH 2024. [[ it's really interesting how this approach uses spatial discretization and particles together to approximate things. I'd like to understand more how that helps. I'm curious what a correct particle-only and grid-only solution would look like. It sounds like the use of the grid isn't really engaged in the first example, it seems more useful in the second and third. It's notable that fluids do behave in ways well-described by both a continuous structure and a particle-based structure at the same t--
