implementation.md

Welcome to the zone of technical details! Here we'll walk through the implementation of the simulation from the ground up.

Tip

Think something can be explained better? Open a pull request! :)

Important

Prerequisites Vector calculus, transformation matrices

Particles and external forces

We'll start with a simple particle simulation. We'll spawn a certain number of particles making up our material. They each have their own position $\mathbf x$ in meters $\text m$, velocity $\mathbf v$ in meters per second $\dfrac {\text m}{\text s}$, and mass $m$ in kilograms $\text{kg}$.

struct Uniforms {
    simulation_timestep: f32,
}

struct ParticleData {
    position: vec3f,
    velocity: vec3f,
    mass: f32,
}

Every timestep (of a fixed duration $\Delta t$ in seconds $\text s$), we'll calculate the force $\mathbf f$ (in newtons $\text N$) on each particle. For now, the only force we'll calculate is a constant gravity force $\mathbf f = \mathbf f_g = m \cdot \langle 0, 0, -9.81 \rangle \text{ N}$ (where +z is our "up" vector).

Note

Lowercase $\mathbf f$ for force here. Uppercase $\mathbf F$ is used for the deformation matrix later on!

Using Newton's second law of motion, we can obtain the acceleration on each particle in meters per second per second $\dfrac {\text m}{\text s^2}$: $\mathbf f = m \mathbf a \implies \mathbf a = \dfrac{\mathbf f}{ m} = \langle 0, 0, -9.81 \rangle \text{ }\dfrac {\text m}{\text s^2}$. (If we wanted, we could also add some user-controlled forces so we can grab and throw the particles from the UI, to make sure our simulation works with all kinds of forces.)

We now have the rate of change of position (velocity) as well as the rate of change of velocity (acceleration). To update our particles' states, we'll use Euler integration, or simply multiplying the rates of change by our timestep $\Delta t$, to obtain the velocity and position of the particle in the next frame:

$$\begin{align*} \mathbf v_{\text{next}} &:= \mathbf v + \mathbf a \cdot \Delta t \\ \mathbf x_{\text{next}} &:= \mathbf x + \mathbf v \cdot \Delta t \end{align*}$$

If we do a timestep or a few every frame, our particles should now fall downward!

Material point method

Now our particles are all moving independently of each other, so they're all acting like separate objects rather than a continuous material.

Material point method (MPM) is a technique that will help us model the cohesion of our material. With MPM, we'll treat each particle as a material point, or a sample of how our imagined continuous, deformable material is moving or deformed at a certain point in space.

MPM will also introduce a grid that we'll modify with our forces, rather than the particles directly. Notably, grid cells are built to not have gaps in between them... which makes a grid strategy fitting for a continuous material! The main drivers of the simulation are still the particles, though, which gives us the flexibility to move around freely without being constrained by the grid.

We can divide the main MPM algorithm into 4 steps:

1. Particle-to-grid (P2G). Transfer the momentum $\mathbf p = m\mathbf v$ and mass of each particle to grid cells that are near the particle.
2. Grid update. Update each grid cell's momentum and mass using our external forces.
3. Grid-to-particle (G2P). Transfer the momentum of each grid cell back to the particles in its vicinity.
4. Grid clear. Zero out all the grid cell momentums and masses, so the next frame can do this process over again.

Before we jump in, let's construct our MPM grid first.

struct GridUniforms {
    grid_min_coords: vec3f, // bottom-left-front corner of the grid
    grid_cell_dims: vec3f, // length/width/height of each cell
    n_cells_per_grid_axis: vec3u, // number of cells in each axis

    fixed_point_scale: f32, // fixed-point scale factor
}

struct CellData {
    momentum_x: atomic<i32>,
    momentum_y: atomic<i32>,
    momentum_z: atomic<i32>,
    mass: atomic<i32>,
}

Note the use of atomics here as well as the inclusion of a fixed-point scale factor. We'll come back to those after describing the particle-to-grid step.

Tip

We started seeing some nice-looking snow breakage effects at a grid cell resolution of 128 in a 4-meter-large grid and 50 000+ particles! If your grid resolution is too low, your snow will move more slowly and deform more like soap bubbles.

Particle-to-grid

First, we need to transfer the momentum and mass of each particle to grid cells that are near the particle. To do this transfer, we'll weight the particle's influence on those grid cells based on the distance of the cell's center from the particle.

We'll need to calculate which grid cell $\verb|cell_number: vec3u|$ each particle belongs to, as well as the fractional positions in the cell (how far along in each direction the particle is within the cell, as a percentage) $\verb|frac_pos: vec3f|$. This can be done knowing the position $\verb|grid_min_coords: vec3f|$ where the grid starts, as well as the size $\verb|grid_cell_dims: vec3f|$ of each cell in each direction $x,y,z$ (which we'll write as a subscript $a$ here):

$$\begin{align*} \verb|cell_number|_a &:= \left\lfloor \frac{\mathbf x_a - \verb|grid_min_coords|_a}{\verb|grid_cell_dims|_a} \right\rfloor \\ \verb|frac_pos|_a &:= \dfrac{\mathbf x_a - \verb|grid_min_coords|_a - \verb|cell_number|_a \cdot \verb|grid_cell_dims|_a}{\verb|grid_cell_dims|_a} \end{align*}$$

For the weights themselves, we'll use a computationally cheap piecewise function for the weights, a quadratic B-spline based on the particle's fractional position in the cell.

To do this, we'll want to take the 27 nearest cells to the particle (or the 3 nearest 2D planes of cells in each axis). Let's then calculate the amount of influence $\verb|weight|_{b, a}$ on the neighboring planes, where $a$ is the axis and $b$ is how many planes over it is in that axis, relative to the particle. We'll store these in $\verb|weight: array<vec3f, 3>|$.

Tip

For example, in the diagram above, $\verb|weight|_{0, y}$ is the weight for the column of cells containing the particle, $\verb|weight|_{-1, y}$ is the weight for the column of cells to the left, and $\verb|weight|_{1, y}$ is the weight for the column of cells to the right.

$$\begin{align*} \verb|weight|_{-1, a} &:= 0.5 \cdot (1 - \verb|frac_pos|_a)^2 \\ \verb|weight|_{0, a} &:= 0.75 - (\verb|frac_pos|_a - 0.5)^2 \\ \verb|weight|_{1, a} &:= 0.5 \cdot (\verb|frac_pos|_a)^2 \end{align*}$$

(All planes farther than 1 plane away—that is, $b < -1$ or $b > 1$—are assigned a weight of $0$.)

For a given plane of cells, these formulae result in the following weights based on the particle's distance from that plane's center in each axis:

(click to open in Desmos)

This should give us 3 weight values $w_x, w_y, w_z$ for each cell (1 for each axis), which we can multiply together to obtain the overall weight $w = w_x \cdot w_y \cdot w_z$ of the cell.

Finally, we can calculate the amount of momentum and mass to transfer to each grid cell:

$$\begin{align*} \verb|grid_cell_momentum| &\verb| += | w \cdot \mathbf p \\ \verb|grid_cell_mass| &\verb| += | w \cdot m \end{align*}$$

Sum up these values for every particle, and we'll be ready to apply forces!

Since we're working with the GPU here, we'll handle one particle per thread in a compute shader. The grid and particle data is handled in storage memory, shared between all threads. However, we'll run into a problem if we run this algorithm as written: we'll need to read from and write to the same grid cells from many different threads, resulting in a race condition.

Recall in our GridCell struct we used atomic<i32>s to store the momentum and mass. Atomics ensure that only one thread will read, write, or read-and-write to a memory location at a time. Since we're accumulating mass and momentum here, we'll use the atomicAdd function to read-and-write to the grid cell without running into conflicts.

Note that WebGPU devices won't allow floating-point-type atomics like atomic<f32> unless it is requested and it is available on the device. To handle this, we'll use fixed-point arithmetic, where we multiply our floats by a constant factor, like $1000$, before converting into an i32 and passing it to the atomic function.

All in all, our algorithm looks something like:

let particle = &particle_data[thread_index];

let mass = (*particle).mass;
let momentum = (*particle).vel * mass;
let pos = (*particle).pos;

let cell_number = vec3i((pos - grid_uniforms.grid_min_coords) / grid_uniforms.grid_cell_dims);
let fractional_pos = pos - grid_uniforms.grid_min_coords - vec3f(cell_number) * grid_uniforms.grid_cell_dims;

for (var offset_z = -1i; offset_z <= 1i; offset_z++) {
    for (var offset_y = -1i; offset_y <= 1i; offset_y++) {
        for (var offset_x = -1i; offset_x <= 1i; offset_x++) {
            let cell_index = linearizeCellIndex(cell_number + vec3i(offset_x, offset_y, offset_z));
            let cell = &grid_data[cell_index];

            let weights = calculateQuadraticBSplineCellWeights(fractional_pos); // see above
            let weight = weights[u32(offset_x + 1)].x
                * weights[u32(offset_y + 1)].y
                * weights[u32(offset_z + 1)].z;

            let momentum_contribution = momentum * weight * grid_uniforms.fixed_point_scale;
            let mass_contribution = mass * weight * grid_uniforms.fixed_point_scale;

            atomicAdd(&(*cell).momentum_x, i32(momentum_contribution.x));
            atomicAdd(&(*cell).momentum_y, i32(momentum_contribution.y));
            atomicAdd(&(*cell).momentum_z, i32(momentum_contribution.z));
            atomicAdd(&(*cell).mass, i32(mass_contribution));
        }
    }
}

Grid update

In this step, we'll apply forces to the grid cells based on their stored momentum and mass. If we already know what those forces $\mathbf f$ are, this step is easy: we just need to add the force to the momentum, since it turns out that $\dfrac{\partial \mathbf p}{\partial t} = \mathbf f$, which means we can Euler-integrate using:

$$\mathbf p_\text{next} := \mathbf p + \mathbf f \cdot \Delta t$$

And that's it!

let cell = &cell_data[thread_index];

let mass = atomicLoad(&(*cell).mass) / grid_uniforms.fixed_point_scale;


let force = vec3f(0, 0, -9.81) * mass;
let momentum_contrib = force * uniforms.simulation_timestep * grid_uniforms.fixed_point_scale;

atomicAdd(&(*cell).momentum_x, i32(momentum_contrib.x));
atomicAdd(&(*cell).momentum_y, i32(momentum_contrib.y));
atomicAdd(&(*cell).momentum_z, i32(momentum_contrib.z));

Grid-to-particle

In this step, we're going to accumulate the new momentums onto the particles based on the grid cells, and then we'll update the particle's position based on that momentum.

To do this, we're going to use our B-spline weights from before as well as the grid loop to determine the influence each grid cell now has on the current particle. After we've added up the velocity contributions from all the grid cells, we can then update the position.

let particle = &particle_data[thread_index];

let mass = (*particle).mass;
let momentum = (*particle).vel * mass;
let pos = (*particle).pos;

let pos_from_grid_min = pos - grid_uniforms.grid_min_coords;
let cell_number = vec3i(pos_from_grid_min / grid_uniforms.grid_cell_dims);
let fractional_pos = pos_from_grid_min - vec3f(cell_number) * grid_uniforms.grid_cell_dims;

let weights = calculateQuadraticBSplineCellWeights(fractional_pos); // see above

var new_particle_velocity = vec3f(0);
for (var offset_z = -1i; offset_z <= 1i; offset_z++) {
    for (var offset_y = -1i; offset_y <= 1i; offset_y++) {
        for (var offset_x = -1i; offset_x <= 1i; offset_x++) {
            let cell_index = linearizeCellIndex(cell_number + vec3i(offset_x, offset_y, offset_z));
            let cell = &grid_data[cell_index];

            let weight = weights[u32(offset_x + 1)].x
                * weights[u32(offset_y + 1)].y
                * weights[u32(offset_z + 1)].z;

            let grid_momentum = vec3f(
                atomicLoad(&(*cell).momentum_x),
                atomicLoad(&(*cell).momentum_y),
                atomicLoad(&(*cell).momentum_z),
            ) / grid_uniforms.fixed_point_scale;
            let grid_mass = atomicLoad(&(*cell).mass);


            new_particle_velocity += grid_momentum / grid_mass * weight;
        }
    }
}

(*particle).vel = new_particle_velocity;
(*particle).pos += (*particle).vel * uniforms.simulation_timestep;

If we'd like, we can also add some boundary conditions here to handle any particles that end up leaving the grid.

Grid clear

Since we add together all the momentums and masses together in the particle-to-grid step, we need to zero out the grid cells before we start accumulating them again for the next timestep. We'll just use some simple atomicStore calls:

let cell = &cell_data[thread_index];

atomicStore(&(*cell).momentum_x, 0);
atomicStore(&(*cell).momentum_y, 0);
atomicStore(&(*cell).momentum_z, 0);
atomicStore(&(*cell).mass, 0);

Next steps

After running the 4 steps repeatedly in a simulation loop, our MPM implementation is complete! Notably, if you add forces other than gravity or set the initial velocities in different directions, you can start to notice the particles bunching up. The velocities of particles will influence nearby particles, making the material look more cohesive, with a goopy or stringy look.

At this point, though, we've only modeled external forces. In fact, if we only have gravity, then we'll just see all the particles fall straght down through the material, which isn't very interesting and doesn't really model any material all that well. We can make things more interesting!

Internal forces

Recall how our material points represent samples of how a continuous material is moving or deformed at each of our particles' positions. It turns out that deformation is the key to adding internal forces!

Deformation

The deformation matrix $\mathbf F$ represents the local transformation of the material at a given material point. By using a matrix, we have the ability to represent shearing, stretching, and rotation of the material.

Note

It'll be helpful here to know how matrices represent linear transformations! In the diagrams below, we're going to draw the deformation matrix as how it transforms the basis vectors of a coordinate system.

Deformation is a property of our particles, so let's update our ParticleData struct:

struct ParticleData {
    pos: vec3f;
    vel: vec3f;
    mass: f32;
    deformation: mat3x3f; // !! NEW
}

When initializing our particles, we'll want to set this to the identity matrix $\mathbf I$, representing no deformation. Over time, deformation will accumulate in this matrix as the material deforms at that point.

Suppose we have an initial state of a continuous material at rest. Let's start out by "baking" or saving, into each point of the material, the current world position of that point as a reference. We'll call this saved reference position $\mathbf x_\text{material}$.

This way, we can keep track of how much any given point in the material has moved at some time later in our simulation. With a continuous field of points, we can even tell at any point in space how much a material has separated, sheared, or squished by looking at nearby points!

Our tool of choice for this is the derivative. In particular, we'll want to differentiate the saved positions $\mathbf x_\text{material}$ by the current world position $\mathbf x$.

$$\mathbf F = \frac{\partial \mathbf x_\text{material}}{\partial \mathbf x}$$

Now, our simulated material isn't exactly continuous (so we can't differentiate it), and it'd be nice to be able to access the deformation without having to do a bunch of expensive sampling (you'd have to loop through a lot of particles to find the closest ones!). So, similar to position, we'll track deformation progressively by accumulating changes in deformation over time, $\dfrac{\partial\mathbf F}{\partial t}$.

Looks like this is a second-order partial derivative:

$$\frac{\partial\mathbf F}{\partial t} = \frac{\partial^2\mathbf x_\text{material}}{\partial\mathbf x\partial t}$$

Now, that looks an awfully lot like we can just say $\dfrac{\partial\mathbf F}{\partial t} = \dfrac{\partial\mathbf v_\text{material}}{\partial\mathbf x}$, where $\mathbf v_\text{material}$ tells us how a given point in the material is moving over time (and $\mathbf v_\text{material}$ is a good choice for us, because we can interpolate that value from our grid). That's most of the way there, actually! We'll just have to be careful since the material has already been transformed a certain way, so we'll need to locally transform something we get by $\mathbf F$ later down the line to compensate for it.

The first thing we'll do is modify the grid-to-particle step to calculate the current change in deformation $\dfrac{\partial \mathbf F}{\partial t}$ with respect to time at each particle. One way to think about this change in deformation is to consider how much the material's velocity $\mathbf v_\text{material}$ varies on opposite sides of the particle.

We can tell that the velocity field above results in a shearing effect, where the area above the particle is moving rightward and the area below is moving leftward. To get our desired change in deformation, we can write the particle's basis vectors for the diagram on the right side, and then take their difference with the basis vectors on the left side:

$$\frac{\partial \mathbf F}{\partial t} = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} - \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}$$

Also note that $\mathbf F$ is not a homogeneous transformation matrix, so it doesn't encode translation. Therefore, even if the velocity vectors result in a net displacement, we only encode the stretching, shearing, and rotation in the deformation matrix:

$$\frac{\partial \mathbf F}{\partial t} = \begin{bmatrix} 1 & 0 \\ 0.25 & 1 \end{bmatrix} - \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & 0 \\ 0.25 & 0 \end{bmatrix}$$

And one more example, this time with compression instead of shearing:

$$\frac{\partial \mathbf F}{\partial t} = \begin{bmatrix} 1 & 0 \\ 0 & 0.6 \end{bmatrix} - \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & 0 \\ 0 & -0.4 \end{bmatrix}$$

We can probably intuit now that this change in deformation is somehow dependent on how the material's velocity vector varies with respect to position along each axis. In calculus terms, we might say we want to take the derivative of the material's velocity field $\mathbf v_\text{material}$ with respect to the position $\mathbf x$. We'll call the result of this the velocity gradient $\dfrac{\partial \mathbf v_\text{material}}{\partial \mathbf x}$. It turns out that the tool for differentiating a vector with respect to another vector is the Jacobian matrix $\mathbf J$, which is simply a matrix such that $\mathbf J_{i,j} = \dfrac{\partial\mathbf v_{\text{material},i}}{\partial\mathbf x_j}$ represents the derivative of the $i\text{th}$ component of $\mathbf v_\text{material}$ with respect to the $j\text{th}$ component of $\mathbf x$. For our 2D case:

$$\frac{\partial \mathbf v_\text{material}}{\partial \mathbf x} = \begin{bmatrix} \dfrac{\partial\mathbf v_{\text{material},x}}{\partial\mathbf x_x} & \dfrac{\partial\mathbf v_{\text{material},x}}{\partial\mathbf x_y}\\ \dfrac{\partial\mathbf v_{\text{material},y}}{\partial\mathbf x_x} & \dfrac{\partial\mathbf v_{\text{material},y}}{\partial\mathbf x_y} \end{bmatrix}$$

...which, for both of the examples above, gives us exactly the values we had written for $\dfrac{\partial \mathbf F}{\partial t}$! So now we have a way to calculate the change in deformation given a velocity field, $\dfrac{\partial \mathbf F}{\partial t} = \dfrac{\partial \mathbf v_\text{material}}{\partial \mathbf x}$.

One major assumption we made above is that we're starting from the identity matrix as our deformation. We'll get different results with a different starting matrix. In general, we'll want to locally transform the change in deformation above by our particle's existing deformation:

$$\begin{align*} \dfrac{\partial \mathbf F}{\partial t} = \dfrac{\partial \mathbf v_\text{material}}{\partial \mathbf x} \cdot \mathbf F & \implies \mathbf F_\text{next} = \mathbf F + \dfrac{\partial \mathbf v_\text{material}}{\partial \mathbf x} \cdot \mathbf F \cdot \Delta t \\ & \implies \mathbf F_\text{next} = \left(\mathbf I + \dfrac{\partial \mathbf v_\text{material}}{\partial \mathbf x} \cdot \Delta t\right) \cdot \mathbf F \end{align*}$$

One problem, though: we don't exactly have a continuous velocity field to differentiate, but a set of discrete velocities at the center of each grid cell. We need some way to interpolate those grid cell velocities into a continuous field.

Let's return to our 3D simulation and try writing out the velocity field $\mathbf v_\text{material}(\mathbf x)$ as a function of $\mathbf x$ so we know what we should be trying to differentiate. Recall that, in the grid-to-particle step, we analyzed how much velocity each single cell contributed to our particle's velocity. Also keep in mind that our material points represent samples of a continuous material's velocity and deformation; by definition, the material's velocity $\mathbf v_\text{material}(\mathbf x)$ at any particle's position is just the velocity of that particle.

So, let's take another look at how we calculated each particle's velocity in the grid-to-particle step to get some more insight:

let weight = weights[u32(offset_x + 1)].x
    * weights[u32(offset_y + 1)].y
    * weights[u32(offset_z + 1)].z;

// ...

new_particle_velocity += grid_velocity * weight; // grid_velocity = grid_momentum / grid_mass

Think about what's going on here: given a specific configuration of grid momentums and masses, a particle's velocity is going to be entirely determined by where it is in the grid. That means the grid_velocity * weight ($=: w \cdot \mathbf v_\text{cell}$) calculation we're doing, summed altogether across all 27 neighbor cells, is what we're defining as the material's velocity $\mathbf v_\text{material}(\mathbf x)$ at any given position $\mathbf x$. So we have a way to express the velocity field as a big sum of 27 cell contributions, and by the addition rule of derivatives, the derivative is just each sum of all the derivatives of all 27 addends. (Let $\mathbf v_\text{cell}$ be the velocity stored in a cell.)

$$\begin{align*} \mathbf v_\text{material} &= \sum_{\Delta x = -1}^1 \sum_{\Delta y = -1}^1 \sum_{\Delta z = -1}^1 w \cdot \mathbf v_{\text{cell}}\text{ (of the cell at $(\Delta x, \Delta y, \Delta z)$)} \\ \dfrac{\partial \mathbf F}{\partial t} = \dfrac{\partial \mathbf v_\text{material}}{\partial \mathbf x} &= \sum_{\Delta x = -1}^1 \sum_{\Delta y = -1}^1 \sum_{\Delta z = -1}^1 \dfrac{\partial [w \cdot \mathbf v_{\text{cell}}]}{\partial \mathbf x} \text{ (of the cell at $(\Delta x, \Delta y, \Delta z)$)} \end{align*}$$

The B-spline functions we used to compute the weights, it turns out, additionally has the purpose of interpolating our velocities here in a smooth way. This makes our function easy to differentiate with respect to position, and we can even get cheap, exact values for the derivatives since the functions are all simple quadratics.

Our B-spline functions were piecewise, so their derivatives are also piecewise. Let's differentiate the 3 pieces analytically with respect to the fractional position in each axis $\mathbf x_a$:

$$\frac{\partial }{\partial \mathbf x_a}\begin{bmatrix} 0.5 \cdot (1 - \mathbf x_a)^2 \\ 0.75 - (\mathbf x_a - 0.5)^2 \\ 0.5 \cdot (\mathbf x_a)^2 \end{bmatrix} = \begin{bmatrix} \mathbf x_a - 1 \\ -2 \cdot (\mathbf x_a - 0.5) \\ \mathbf x_a \end{bmatrix}$$

Just as we selected which weight $w_x, w_y, w_z$ to use for each plane of cells based on the position of that plane from a particle, we'll select which weight derivative in the same manner.

Now back to the overall weight for a cell. After we multiplied the three weights $w = w_x \cdot w_y \cdot w_z$, we were able to calculate its final velocity contribution $w \cdot \mathbf v_\text{cell}$ to each particle. $\mathbf v_\text{cell}$ for a given cell is constant during the entire grid-to-particle step, so we'll mainly need to worry about differentiating $w$ with respect to position, giving us what we'll call the weight gradient $\dfrac{\partial w}{\partial \mathbf x}$. We'll need some more Jacobians for this:

$$\begin{align*} \frac{\partial w}{\partial \mathbf x} &= \begin{bmatrix} \dfrac{\partial w}{\partial \mathbf x_x} & \dfrac{\partial w}{\partial \mathbf x_y} & \dfrac{\partial w}{\partial \mathbf x_z} \end{bmatrix} \\ \frac{\partial [w \cdot \mathbf v_\text{cell}]}{\partial \mathbf x} &= \begin{bmatrix} \dfrac{\partial w}{\partial \mathbf x_x}\cdot \mathbf v_{\text{cell},x} & \dfrac{\partial w}{\partial \mathbf x_y}\cdot \mathbf v_{\text{cell},x} & \dfrac{\partial w}{\partial \mathbf x_z}\cdot \mathbf v_{\text{cell},x} \\ \dfrac{\partial w}{\partial \mathbf x_x}\cdot \mathbf v_{\text{cell},y} & \dfrac{\partial w}{\partial \mathbf x_y}\cdot \mathbf v_{\text{cell},y} & \dfrac{\partial w}{\partial \mathbf x_z}\cdot \mathbf v_{\text{cell},y} \\ \dfrac{\partial w}{\partial \mathbf x_x}\cdot \mathbf v_{\text{cell},z} & \dfrac{\partial w}{\partial \mathbf x_y}\cdot \mathbf v_{\text{cell},z} & \dfrac{\partial w}{\partial \mathbf x_z}\cdot \mathbf v_{\text{cell},z} \end{bmatrix} \\ &= \begin{bmatrix} \dfrac{\partial w}{\partial \mathbf x_x}\cdot \mathbf v_\text{cell} & \dfrac{\partial w}{\partial \mathbf x_y}\cdot \mathbf v_\text{cell} & \dfrac{\partial w}{\partial \mathbf x_z}\cdot \mathbf v_\text{cell} \end{bmatrix} &\text{(shorthand)} \end{align*}$$

It'll help us here to consider the overall weight in terms of the 3 individual weight factors in each axis $w_x, w_y, w_z$. Note that only $w_x$ will vary with $\mathbf x_x$, only $w_y$ will vary with $\mathbf x_y$, and only $w_z$ will vary with $\mathbf x_z$. So we can further rewrite our Jacobian as:

$$\begin{align*} \frac{\partial w}{\partial \mathbf x} &= \begin{bmatrix} \dfrac{\partial w_x}{\partial \mathbf x_x} \cdot w_y \cdot w_z & w_x \cdot \dfrac{\partial w_y}{\partial \mathbf x_y} \cdot w_z & w_x \cdot w_y \cdot \dfrac{\partial w_z}{\partial \mathbf x_z} \end{bmatrix} \\ \frac{\partial [w \cdot \mathbf v_\text{cell}]}{\partial \mathbf x} &= \begin{bmatrix} \dfrac{\partial w_x}{\partial \mathbf x_x} \cdot w_y \cdot w_z \cdot \mathbf v_\text{cell}& w_x \cdot \dfrac{\partial w_y}{\partial \mathbf x_y} \cdot w_z \cdot \mathbf v_\text{cell}& w_x \cdot w_y \cdot \dfrac{\partial w_z}{\partial \mathbf x_z} \cdot \mathbf v_\text{cell} \end{bmatrix} \end{align*}$$

...where $\dfrac{\partial w_x}{\partial \mathbf x_x}, \dfrac{\partial w_y}{\partial \mathbf x_y}, \dfrac{\partial w_z}{\partial \mathbf x_z}$ are just the weight derivatives we calculated above!

At last, we now have all the pieces we need to calculate the new deformation gradient, again doing some more Euler-integration:

$$\mathbf F_\text{next} := \left(\mathbf I + \frac{\partial \mathbf v_\text{material}}{\partial \mathbf x} \cdot \Delta t\right) \cdot \mathbf F$$

let weights = calculateQuadraticBSplineCellWeights(fractional_pos); // see above
let weight_derivs = calculateQuadraticBSplineCellWeightDerivatives(fractional_pos); // see above

var total_velocity_gradient = mat3x3f();
for (var offset_z = -1i; offset_z <= 1i; offset_z++) {
    for (var offset_y = -1i; offset_y <= 1i; offset_y++) {
        for (var offset_x = -1i; offset_x <= 1i; offset_x++) {
            let cell_index = linearizeCellIndex(cell_number + vec3i(offset_x, offset_y, offset_z));
            let cell = &grid_data[cell_index];

            let weight = weights[u32(offset_x + 1)].x
                * weights[u32(offset_y + 1)].y
                * weights[u32(offset_z + 1)].z;

            // ...

            let grid_velocity = grid_momentum / grid_mass;

            let weight_gradient_contrib = vec3f(
                weight_derivs[u32(offset_x + 1)].x * weight[u32(offset_y + 1)].y * weight[u32(offset_z + 1)].z,
                weight[u32(offset_x + 1)].x * weight_derivs[u32(offset_y + 1)].y * weight[u32(offset_z + 1)].z,
                weight[u32(offset_x + 1)].x * weight[u32(offset_y + 1)].y * weight_derivs[u32(offset_z + 1)].z,
            );

            total_velocity_gradient += mat3x3f(
                weight_gradient_contrib.x * grid_velocity,
                weight_gradient_contrib.y * grid_velocity,
                weight_gradient_contrib.z * grid_velocity,
            );
        }
    }
}

(*particle).deformation += (IDENTITY_MAT3 + total_velocity_gradient * uniforms.simulation_timestep) * (*particle).deformation;

Even after all that, though, our deformation isn't actually doing anything to the simulation. Let's visualize the fruits of our labor by adding a force that comes as a direct result of deformation.

Elasticity and stress force

We noted before that all of our particles, with only a gravity force, are going to sink through the material, all the way to the bottom of the grid. The stress force will counteract this by attempting to restore the material's deformation to the original, undeformed state.

Note

TBD

Plasticity

Note

TBD

Hardening

Note

TBD

Rigid colliders

Note

TBD

Position-based material point method (PBMPM)

Note

TBD

Raymarching

Note

TBD