Interpolating Functions

In each time step of the Material Point Method (MPM), particles transfer their mass and momentum to the background grid, and later retrieve updated velocities from the grid to perform advection. Both transfers rely on interpolating functions to determine how each particle interacts with nearby grid nodes.

We define the interpolating function associated with grid node $i$ as $N_{i} (x)$ . Here, $i$ is a multi-index - $(i, j)$ in 2D or $(i, j, k)$ in 3D - and $x$ is an arbitrary spatial position. When evaluating this function at a particle location $x_{p}$ , we use the shorthand:

$w_{i p} = N_{i} (x_{p}) .$

This weight determines how strongly particle $p$ influences grid node $i$ : the closer the particle is to the node, the larger $w_{i p}$ becomes.

Example 27.2.1 (Linear Interpolation). In 1D, the simplest choice of interpolating functions is the linear (tent) function: $N (x) = {1 - ∣ x ∣ 0 if 0 \leq ∣ x ∣ < 1, otherwise .$ To apply this in a grid-based setting, we scale the spatial coordinates relative to the grid spacing $h$ (i.e., the distance between adjacent grid nodes along one axis). The scaled version becomes: $N (\frac{x _{p} - x _{i}}{h}) .$ In higher dimensions (2D, 3D), we use the tensor product of 1D functions: $N_{i} (x_{p}) = N (\frac{x _{p} - x _{i}}{h}) \cdot N (\frac{y _{p} - y _{j}}{h}) \cdot N (\frac{z _{p} - z _{k}}{h}) .$

Linear interpolation is computationally efficient and easy to implement, and is widely used in FLIP [Brackbill et al. 1988] [Bridson 2015] for fluid simulation. However, it suffers from two serious issues when applied in MPM:

Its gradient $\nabla N (x)$ is discontinuous, leading to unstable and noisy force evaluations.
When particles are near cell boundaries, $w_{i p}$ becomes small while $\nabla w_{i p}$ remains large, causing numerical instabilities known as cell-crossing instability.
The function’s support is too narrow, making it prone to numerical fracture when neighboring particles do not share enough grid nodes.

Because of these issues, linear interpolation is typically avoided in MPM, especially in solid simulation.

Example 27.2.2 (Quadratic B-Spline Interpolation). A better choice is the quadratic B-spline interpolating functions, which is $C^{1}$ -continuous (with continuous gradients) and has wider support: $N (x) = ⎩ ⎨ ⎧ \frac{3}{4} - x^{2} \frac{1}{2} (\frac{3}{2} - ∣ x ∣)^{2} 0 if ∣ x ∣ < \frac{1}{2}, if \frac{1}{2} \leq ∣ x ∣ < \frac{3}{2}, otherwise .$

Quadratic B-splines have a support region of width $3 h$ and provide smooth, stable interactions between particles and the grid. They are computationally efficient and require less memory than higher-order B-splines.

Example 27.2.3 (Cubic B-Spline Interpolation). For even better smoothness and broader support, we can use the cubic B-spline interpolating functions: $N (x) = ⎩ ⎨ ⎧ \frac{1}{2} ∣ x ∣^{3} - ∣ x ∣^{2} + \frac{2}{3} \frac{1}{6} (2 - ∣ x ∣)^{3} 0 if 0 \leq ∣ x ∣ < 1, if 1 \leq ∣ x ∣ < 2, otherwise .$

The cubic B-spline has support over $[- 2 h, 2 h]$ , making it more robust to numerical fracture and instabilities. However, the wider support also means that each particle interacts with more grid nodes, which increases the computational overhead during transfers (e.g., mass, momentum, and force). This trade-off between smoothness and efficiency should be considered when choosing an interpolation kernel.

Gradient of the Interpolating Function

Gradients are needed during force computation. For high-dimensional interpolating functions defined using tensor products, we compute their gradient using the chain rule: $\nabla N_{i} (x_{p}) = \frac{1}{h} N^{'} (\frac{x _{p} - x _{i}}{h}) N (\frac{y _{p} - y _{j}}{h}) N (\frac{z _{p} - z _{k}}{h}) N (\frac{x _{p} - x _{i}}{h}) \frac{1}{h} N^{'} (\frac{y _{p} - y _{j}}{h}) N (\frac{z _{p} - z _{k}}{h}) N (\frac{x _{p} - x _{i}}{h}) N (\frac{y _{p} - y _{j}}{h}) \frac{1}{h} N^{'} (\frac{z _{p} - z _{k}}{h}) .$ Here, $N^{'} (x)$ is the derivative of the 1D interpolating function $N (x)$ .