Skip to content

Material Point Method

Reference

The Material Point Method for Simulating Continuum Materials, Chenfanfu Jiang, SIGGRAPH 2016 Course Notes Version 1 (May 2016).

Preliminary

Contimuum Motion

Here \(d\) denotes dimension of the space, usually \(d=2,3\).

We denote \(\pmb{x}\in \Omega^t\subset \mathbb{R}^d\) world (deformed) space(coordinates), current position. Physically, we focus on a fixed position in the space and measurethe velocity of whichever particla that is passing by the position.

Then we denote \(\pmb{X}\in \Omega^0\subset \mathbb{R}^d\) material (undeformed) space(coordinates), or initial position. Physically, we measure velocity from a fixed particle, which has its mass and occupies some volume since the beginning.

Define

\[ \pmb{x}=\phi(\pmb{X},t) \]

which is usually a bijection. This is associated with the assumption that no two different particles of material ever occupy the same space at the same time. So

\[ \forall \pmb{x}\in \Omega^t, \exists ! \pmb{X}\in \Omega^0, \text{ s.t. } \phi(\pmb{X},t)=\pmb{x}. \]

Example

For a given initial position \(\pmb{X}\), we have some common function of \(\pmb{x}\).

(i) \(\pmb{x}=\pmb{X}+tv(t)\vec{N}\)

discrete form at time \(n\):

\[ \pmb{x}^{(n)}=\pmb{x}^{(n-1)}+\Delta t v^{(n)} \vec{N}^{(n)} \]

(ii) \(\pmb{x}=R\pmb{X}+\vec{b}\)

discrete form at time \(n\):

\[ \pmb{x}^{(n)}=R^{(n-1)}\pmb{x}^{(n-1)}+\vec{b}^{(n)} \]

The velocity of a given material point \(\pmb{X}\) at time \(t\) is a mapping \(V(\cdot,t):\Omega^0\rightarrow \mathbb{R}^d\) defined by

\[ \pmb{V}(\pmb{X},t)\overset{\Delta}{=}\frac{\partial \phi(\pmb{X},t)}{\partial t} \]

and acceleration \(A(\cdot,t):\Omega^0\rightarrow \mathbb{R}^d\)

\[ \begin{align*} \pmb{A}(\pmb{X},t)&\overset{\Delta}{=}\frac{\partial \pmb{V}(\pmb{X},t)}{\partial t}\\ &=\frac{\partial^2 \phi(\pmb{X},t)}{\partial t^2} \end{align*} \]

Deformation

For every small unit in a material, it has a deformation mapping :

\[ \pmb{F}(\pmb{X},t)\overset{\Delta}{=}\frac{\partial \pmb{x}}{\partial \pmb{X}} \in \mathbb{R}^{d\times d} \]

Note that \(F\) is also related with both \(\pmb{X}\) and \(t\). And \(J=\det{(F)}\) characterizes infinitesimal volume change.

In Example we have \(\pmb{F}=I\) (i) and \(\pmb{F}=R\) (ii).

The above is rigid transformation, and local.

Push forward and Push back

Push forward of a function, is often referred to as Eulerian (a function of \(\pmb{x}\)). That is, given \(G: \Omega^0\rightarrow \mathbb{R}\), the push forward \(g(\cdot, t):\Omega^t\rightarrow \mathbb{R}\) is defined

\[ g(\pmb{x})=G(\phi^{-1}(\pmb{x}),t) \]

Push back function is often referred to as Lagrangian (a function of \(\pmb{X}\)). That is, given \(g:\Omega^t\rightarrow \mathbb{R}\), the push back \(G(\cdot,t):\Omega^0\rightarrow \mathbb{R}\) is defined

\[ G(\pmb{X})=g(\phi(\pmb{x}), t) \]

But we usually do not use push back but push forward more frequently since we can always have access to the grid velocity and acceleration based on the particles around it.

So it is useful to define Eulerian counterparts. The velocity

\[ \pmb{v}(\pmb{x},t)=\pmb{V}(\phi^{-1}(\pmb{x},t),t) \]

and the acceleration

\[ \pmb{a}(\pmb{x},t)=\pmb{A}(\phi^{-1}(\pmb{x},t),t) \]

The following result is really important:

\[ \begin{align*} a_i(\pmb{x},t)&=A_i(\phi^{-1}(\pmb{x},t), t)\\ &=\frac{\partial V_i}{\partial t}(\phi^{-1}(\pmb{X}, t),t)\\ &=\frac{\partial v_i}{\partial t}(\pmb{x}, t)+\sum_{j=1}^d\frac{\partial v_i}{\partial x_j}(\pmb{x},t)v_j(\pmb{x},t)\\ &\overset{\Delta}{=}\frac{D}{Dt}v_i(\pmb{x},t) \end{align*} \]

So

\[ \pmb{a}=\frac{D}{Dt}\pmb{v} \]

More generally, for a general Eulerian function \(f(\cdot, t):\Omega^t\rightarrow \mathbb{R}\), we can use this same notation to mean

\[ \frac{D}{Dt}f(\pmb{x},t)=\frac{\partial f}{\partial t}(\pmb{x}, t)+\sum_{j=1}^d\frac{\partial f}{\partial x_j}(\pmb{x},t)v_j(\pmb{x},t) \]

which is called material derivative, and the first item of derivative is called local rate of change. The above one is also the push forward of \(\partial F/\partial t\) where \(F\) is a Lagrangian function with \(F(\cdot,t):\Omega^0\rightarrow \mathbb{R}\).

Relationship between two Deformation Gradients

For deformation gradient, most of the time in the physics of a material, the Lagrangian view is the dominant one. There is however a useful evolution of the Eulerian (push forward) of \(\pmb{F}(\cdot,t):\Omega^0\rightarrow \mathbb{R}^{d\times d}\). If we denote \(\pmb{f}(\cdot, t):\Omega^t\rightarrow \mathbb{R}^{d\times d}\) be the push forward of \(\pmb{F}\), then

\[ \frac{D}{Dt}\pmb{f}(\pmb{x},t) = \frac{\partial \pmb{v}}{\partial \pmb{x}}\pmb{f} \]

or

\[ \dot{\pmb{F}}=(\nabla \pmb{v})\pmb{F} \]

cause we have

\[ \begin{align*} \frac{\partial}{\partial t}F_{ij}(\pmb{X}, t)&=\frac{1}{\partial t }\frac{\partial \phi_i}{\partial X_j}(\pmb{X},t)\\ &=\frac{1}{\partial X_j}\frac{\partial \phi_i}{ \partial t }(\pmb{X},t)\\ &=\frac{\partial V_i}{\partial X_j}(\pmb{X},t)\\ &=\frac{\partial v_i}{\partial X_j}(\phi(\pmb{X},t),t)\\ &=\sum_{k=1}^d\frac{\partial v_i}{\partial x_k}(\phi(\pmb{x},t),t)\cdot \frac{\partial \phi_k}{\partial X_j}(\pmb{X},t)\\ &=\sum_{k=1}^d\frac{\partial v_i}{\partial x_k}(\phi(\pmb{x},t),t)\cdot F_{kj}(\pmb{X},t) \end{align*} \]

The above equation will play an important role in deriving the discretized deformation gradient \(F\) update on each MPM particle.

Volume and Area Change

  • Volume

Consider \(dV\) being defined over the standard basis vectors \(\vec{e}_i\), \(i=1.2.3\), with \(dV=dL_1\vec{e}_1\cdot (dL_2\vec{e}_2 \times dL_3\vec{e}_3)\). If we denote \(d\pmb{L}_i=dL_i\vec{e}\), then

\[ dV=dL_1dL_2dL_3, \quad \]

cause we have \(d\pmb{l}_i=\pmb{F}d\pmb{L}_i\) (like \((\pmb{x}_2-\pmb{x}_1)=\pmb{F}(\pmb{X_2}-\pmb{X}_1)\)), so it can be shown that \(dl_1dl_2dl_3=JdL_1dL_2dL_3\) or

\[ dv=JdV. \]

Based on the above property we can have a great weapon in proof. Given function \(G(\pmb{X})\) or \(g(\pmb{x},t)\) we have

\[ \int_{B^t}g(\pmb{x})d\pmb{x}=\int_{B^0}G(\pmb{X})J(\pmb{X},t)d\pmb{X} \]

Actually it is the variable substitution formula of multiple integral.

  • Areas

And similar analysis can be done for areas.

\[ \begin{align*} dv&=JdV\\ \pmb{n}ds\cdot d\pmb{l}&=J\pmb{N}dS\cdot d\pmb{L}\\ \pmb{n}ds\cdot \pmb{F}d\pmb{L}&=J\pmb{N}dS\cdot d\pmb{L}\quad \text{using $d\pmb{l}=\pmb{F}d\pmb{L}$}\\ \pmb{n}ds\cdot \pmb{F}&=J\pmb{N}dS \end{align*} \]

So

\[ \begin{equation} \pmb{n}ds=J\pmb{F}^{-T}\pmb{N}dS.\label{eq-area} \end{equation} \]

which is also called Nansom's formula.

Piola-Kirchhoff Stress

Consider a vector element of surface in the material world(Lagrangian), \(\pmb{N}dS\) and after deformation, the material particles making up this area now occupy the element defined by \(\pmb{n}ds\), where \(ds\) is the area and \(\pmb{n}\) is the normal vector in the current configuration(Eulerian).

Then by definition of the Cauchy stress

\[ d\pmb{f}=\pmb{\sigma}\pmb{n}ds \]

The first Piola-Kirchhoff stress tensor \(\pmb{P}\) (PK1 stress, Nominal Stress tensor) is defined by

\[ d\pmb{f}=\pmb{P}\pmb{N}dS \]

which relates the force acting in the current configuration to the surface element in the reference configuration.

So similarly with (Cauchy) traction vector was defined by

\[ \pmb{t}=\frac{d\pmb{f}}{ds} \]

we can introduce a PK1 traction vector

\[ \pmb{T}=\frac{d\pmb{f}}{dS}, \quad \pmb{T}=\pmb{P}\pmb{N} \]

which is a fictitious quatity.

Note that \(d\pmb{f}=\pmb{t}ds=\pmb{T}dS\), so \(\pmb{t}\) and \(\pmb{T}\) at the same area has the same direction but different magnitudes.

  • Relation between the Cauchy and PK1 Stresses

From the above definition,

\[ \begin{align*} \pmb{P}\pmb{N}dS&=\pmb{\sigma}\pmb{n}\\ \pmb{P}\pmb{N}dS&=\pmb{\sigma}J\pmb{F}^{-T}\pmb{N}dS\quad\text{using equation $\ref{eq-area}$}\\ \pmb{P}&=J\pmb{\sigma}\pmb{F}^{-T}\\ \pmb{\sigma}&=J^{-1}\pmb{P}\pmb{F}^T \end{align*} \]

Hyperelasticity

Stress is related to strain(Deformation gradient \(F\)) through "Constitutive Relationship".

For perfect hyperelastic materials, the relation is defined through the potential energy, which increases with non-rigid deformation from the initial state. That is, The elastic solids whose first Piola-Kirchoff stress \(\pmb{P}\) can be derived from an strain energy density function \(\Psi(\pmb{F})\) (a scalar function) via

\[ \pmb{P}=\frac{\partial \Psi}{\partial \pmb{F}} \]

and Cauchy stress with respect to \(\Psi(\pmb{F})\)

\[ \sigma=\frac{1}{J}\pmb{P}\pmb{F}^T=\frac{1}{J}\frac{\partial \Psi}{\partial \pmb{F}}\pmb{F}^T \]

Energy density function

  • Neo-Hookean
\[ \Psi(\pmb{F})=\frac{\mu}{2}(\text{tr}{(\pmb{F}^T\pmb{F})}-d)-\mu\log(J)+\frac{\lambda}{2}\log^2(J) \]

where \(d\) denoted dimension, \(2\) or \(3\) in practice and

\[ \mu=\frac{E}{2(1+v)},\quad \lambda=\frac{Ev}{(1+v)(1-2v)} \]

in which \(E\) is Young's modulus and \(v\) is Poisson's ratio.

  • Fixed corotated model

This is derived from the Singular Value Decomposition (SVD). Assuming the polar SVD \(\pmb{F}=U\Sigma V^T\), then the energy for fixed corotated model is

\[ \Psi(\pmb{F})=\hat{\Psi}(\Sigma(\pmb{F}))=\mu\sum_{i=1}^d(\sigma_i-1)^2+\frac{\lambda}{2}(J-1)^2 \]

where \(J=\prod\limits_{i=1}^d\sigma_i\) and

\[ \pmb{P}(\pmb{F})=\frac{\partial\Psi}{\partial \pmb{F}}(\pmb{F})=2\mu(\pmb{F}-\pmb{R})+\lambda(J-1)J\pmb{F}^{-T}. \]

Governing Equations

Let

\[ \pmb{V}(\pmb{X},t)=\frac{\partial \phi(\pmb{X},t)}{\partial t}=\frac{\partial{\pmb{x}}}{\partial t} \]

be the velocity defined over \(X\). Then from Lagrangian view of point, the equations are

\[ \begin{cases} \displaystyle R(\pmb{X},t)J(\pmb{X},t)=R(\pmb{X},0) \quad &\text{Conservation of mass},\\ \displaystyle R(\pmb{X},t)\frac{\partial \pmb{V}}{\partial t}=\nabla^{\pmb{X}}\cdot \pmb{P}+R(\pmb{X},0)g \quad &\text{Conservation of momentum}. \end{cases} \]

where \(R\) is the Lagrangian mass density which is related to the more commonly used Eulerian mass density \(\rho\). Note that mass conservation can also be written as

\[ \frac{\partial }{\partial t}[R(\pmb{X},t)J(\pmb{X},t)]=0 \]

In Eulerian view, the governing equations are

\[ \begin{cases} \displaystyle \frac{D}{Dt}\rho(\pmb{x},t)+\rho(\pmb{x},t)\nabla^{\pmb{x}}\cdot \pmb{v}(\pmb{x},t)=0,\quad &\text{Conservation of mass},\\ \displaystyle \rho(\pmb{x},t)\frac{D\pmb{v}}{Dt}=\nabla^{\pmb{x}}\cdot \sigma +\rho(\pmb{x},t)g, \quad &\text{Conservation of momentum}.\\ \end{cases} \]

where

\[ \frac{D}{Dt}=\frac{\partial}{\partial t}+\pmb{v}\cdot \nabla^{\pmb{x}} \]

Conservation of Mass

To be more specific, we have two ways to get the Continuity Equation.

By using multiple integral for density.

Notice that

\[ \rho(\pmb{x},t)=\lim_{\varepsilon\rightarrow 0^+}\frac{\text{mass}(B_{\varepsilon}^t)}{\int_{B_{\varepsilon}^t}d\pmb{x}} \]

We choose a fixed volume, in which the rate of increase of mass must equal the rate at which mass is flowing into the volume through its bounding surface.

The rate of increase mass in a fixed volume \(v\) is

\[ \frac{\partial m}{\partial t}=\frac{\partial}{\partial t}\int_{v}\rho(\pmb{x},t)dv=\int_v \frac{\rho(\pmb{x},t)}{\partial t}dv \]

while the mass flux out through the surface is given by

\[ \int_s \rho \pmb{v}\cdot \pmb{n} ds \]

here \(s\) can denote the overall outer surface of the material(we consider in and out).

so combine the above two we get

\[ \int_v \frac{\rho(\pmb{x},t)}{\partial t}dv + \int_s \rho \pmb{v}\cdot \pmb{n} ds=0 \]

using divergence theorem, we get

\[ \begin{align*} \int_v \frac{\rho(\pmb{x},t)}{\partial t}dv + \int_v \nabla^{\pmb{x}}\cdot(\rho \pmb{v}) dv&=0\\ \frac{\rho(\pmb{x},t)}{\partial t}+ \nabla^{\pmb{x}}\cdot(\rho \pmb{v}) &=0\\ \frac{\partial \rho(\pmb{x},t)}{\partial t}+ (\nabla^{\pmb{x}}\cdot\rho)\pmb{v} +(\nabla^{\pmb{x}}\cdot\pmb{v})\rho &=0\\ \frac{D\rho}{Dt} +\rho \nabla^{\pmb{x}}\cdot(\pmb{v})&=0 \end{align*} \]

We check the fixed volume in terms of material point. Note \(R(\pmb{X},t)=\rho(\phi(\pmb{X},t),t)\), then we see the mass from initial time to time \(t\) must be the same:

\[ \begin{align*} \text{mass}(v)&=\text{mass}V_0\\ \int_{v}\rho(\pmb{x},t)d\pmb{x}&=\int_{V}R(\pmb{X},0)d\pmb{X}\\ \Rightarrow \int_{V}R(\pmb{X},t)J(\pmb{X},t)d\pmb{X}&=\int_{V}R(\pmb{X},0)d\pmb{X} \end{align*} \]

So

\[ R(\pmb{X},t)J(\pmb{X},t)=R(\pmb{X},0), \quad \forall \pmb{X}\in \Omega^0, t\geq 0 \]

Note that \(J(\pmb{X},0)=1\), so \(R(\pmb{X},t)J(\pmb{X},t)=R(\pmb{X},0)J(\pmb{X},0)\), that is,

\[ \frac{\partial }{\partial t}[R(\pmb{X},t)J(\pmb{X},t)]=0. \]

Conservation of Momentum

Also, we have two ways to consider.

Use the similar logic in proof of mass conservation.

This is also the spatial form.

We know that

\[ \pmb{t}(\pmb{x},\pmb{n},t)=\pmb{\sigma}(\pmb{x},t)\pmb{n} \]

We consider a fixed mass, so the space occupied by this matter may change over time.

If \(\pmb{v}\) denotes the Eulerian velocity, then the linear momentum of Euler can be denoted as

\[ \pmb{L}(t)=\int_{B_{\varepsilon}^t} \rho\pmb{v}d\pmb{x}, \]

where

then by \(\pmb{f}^{ext}=\frac{d(m\pmb{v})}{dt}\) formulated by Euler,

\[ \begin{align*} \int_{\partial B_{\varepsilon}^t}\pmb{\sigma}\pmb{n}ds+\int_{B_{\varepsilon}^t}\pmb{f}^{ext}d\pmb{x}&=\frac{d}{dt}\int_{B_{\varepsilon}^t}\rho \pmb{v}d\pmb{x}\\ \int_{B_{\varepsilon}^t}\nabla^{\pmb{x}}\cdot \sigma+\pmb{f}^{ext} d\pmb{x}&=\frac{d}{dt}\int_{B_{\varepsilon}^0}R\pmb{V}Jd\pmb{X}\\ &=\int_{B_{\varepsilon}^0}R\pmb{A}Jd\pmb{X}\\ &=\int_{B_{\varepsilon}^t}\rho\pmb{a}d\pmb{x} \end{align*} \]

Note \(\pmb{V}(\pmb{X},t)=\pmb{v}(\phi(\pmb{X},t),t)\), so we have a fixed volume \(V\) with momentum

\[ \pmb{L}(t) = \int_V R(\pmb{X},t)V(\pmb{X},t)dV \]

So

\[ \begin{align*} \frac{\partial}{\partial t}\int_V R(\pmb{X},t)V(\pmb{X},t)dV &=\int_S\pmb{T}dS +\int_V \pmb{F}dV\\ \int_V R(\pmb{X},t)\frac{\partial}{\partial t}V(\pmb{X},t)dV&= \int_S \pmb{P}\cdot\pmb{N}dS+\int_V\pmb{F} dV\\ \int_V R(\pmb{X},t)\frac{\partial}{\partial t}V(\pmb{X},t)dV&= \int_V \nabla^{\pmb{X}}\cdot \pmb{P}+\pmb{F} dV\\ \Rightarrow R(\pmb{X},t)\frac{\partial}{\partial t}V(\pmb{X},t) &=\nabla^{\pmb{X}}\cdot \pmb{P}+\pmb{F} \end{align*} \]

Material particles

Recall that the material point method is Lagrangian in the sense that we track actual particles of material. That is we keep track of mass (\(m_p\)), velocity (\(v_p\)) and position (\(x_p\)) for a collection of material particles \(p\).

However, all stress based forces are computed on the Eulerian grid, so we have to transfer the material state to the Eulerian configuration to incorporate the effects of material forces.

Then, we transfer these effects back to the material particles and move them in the normal Lagrangian way. The Lagrangian nature makes advection very trivial compared to pure Eulerian methods (such as grid-based fluid simulation).

Eulerian Interpolating Functions

We can denote the interpolation function at grid node \(\pmb{i}=(i,j,k)\) evaluated at a particle location \(\pmb{x}_p\) with

\[ N_{\pmb{i}}(\pmb{x}_p)=N\left(\frac{1}{h}(\pmb{x}_p-\pmb{x}_{\pmb{i}})\right)N\left(\frac{1}{h}(\pmb{y}_p-\pmb{y}_{\pmb{i}})\right)N\left(\frac{1}{h}(\pmb{z}_p-\pmb{z}_{\pmb{i}})\right) \]

where \(h\) is the grid spacing. We can define diffenrent kernel \(N:\mathbb{R}\rightarrow \mathbb{R}\).

Common Kernel \(N\)

  • cubic kernel. It is more expensive but provide wider coverage, thus less sensitive to numerical errors.
\[ N(x)=\begin{cases} \displaystyle \frac{1}{2}|x|^3-|x|^2+\frac{2}{3}, \quad &0\leq |x| < 1\\ \displaystyle \frac{1}{6}(2-|x|)^3,\quad &1\leq |x|< 2\\ \displaystyle 0,\quad &2\leq|x| \end{cases} \]
  • quadratic kernel. It is more computational efficient and memory saving.
\[ N(x)=\begin{cases} \displaystyle \frac{3}{4}-|x|^2,\quad & \displaystyle 0\leq |x| < \frac{1}{2}\\ \displaystyle \frac{1}{2}\left(\frac{3}{2}-|x|\right)^2,\quad &\displaystyle \frac{1}{2}\leq |x|<\frac{3}{2}\\ \displaystyle 0,\quad &\displaystyle \frac{3}{2}\leq |x| \end{cases} \]

Then the gradient of function \(N_{\pmb{i}}(\pmb{x}_p)\) is

\[ \nabla N_{\pmb{i}}(\pmb{x}_p) =\sum_{k=1}^d\left[N'\left(\frac{1}{h}(x_k-x_{\pmb{i}})\right)\prod_{j=1\atop j\neq k}^d N\left(\frac{1}{h}(x_{j}-x_{\pmb{i}})\right)\right] \]

where \(x_k\), \(x_j\) denote the component index of \(\pmb{x}_p\), i.e. \(x_k, x_j\in\{x_p,y_p,z_p\}\).

Eulerian/Lagrangian Mass

Mass of the particle

\[ m_p^n=\int_{B_{\Delta x,p}^{t^n}}\rho(\pmb{x},t^n)d\pmb{x} \]

and define the mass from Eulerian perspective

\[ m_{\pmb{i}}=\sum_pm_p\cdot N_{\pmb{i}}(\pmb{x}_p) \]

Easy to see that

\[ \sum_{\pmb{i}}m_{\pmb{i}}=\sum_p m_p \]

since the weight function \(N_\pmb{i}(\pmb{x}_p)\) is normalized to \(1\).

Eulerian Momentum

Similarly, we transfer particle monentum \(m_p\pmb{v}_p\) to the grid

\[ (m\pmb{v})_{\pmb{i}}=\sum_p m_p \pmb{v}_p N_{\pmb{i}}(\pmb{x}_p) \]

Also easy to that

\[ \sum_{\pmb{i}}(m\pmb{v})_{\pmb{i}}=\sum_{p}m_p\pmb{v}_p \]

and the Eulerian velocity \(\pmb{v}_\pmb{i}\) is defined as

\[ \pmb{v}_{\pmb{i}}=\frac{(m\pmb{v})_{\pmb{i}}}{m_{\pmb{i}}} \]

Eulerian to Lagrangian Transfer

We do not need to transfer mass from the grid to the particles since Lagrangian particle mass never changes. But the Velocity is simply interpolated as

\[ \pmb{v}_p =\sum_{\pmb{i}}\pmb{v}_{\pmb{i}}N_{\pmb{i}}(\pmb{x}_p) \]

Easy to see that

\[ \sum_p m_p \pmb{v}_p=\sum_{\pmb{i}}m_i\pmb{v}_{\pmb{i}} \]

Discretization

Explicite time Integration

APIC Transfers

Transfer from particle to grid

\[ \begin{cases} m_i=\sum_p w_{ip}m_p\\ m_i\pmb{v}_i=\sum_p w_{ip}m_p(\pmb{v}_p+\pmb{B}_p(\pmb{D}_p)^{-1}(\pmb{x}_i-\pmb{x}_p)) \end{cases} \]

where \(\pmb{B}_p\) is a matrix quatity stored at each particle(like mass, position and velocity), \(\pmb{D}_p\) is given by

\[ \pmb{D}_p=\sum_i w_{ip}(\pmb{x}_i-\pmb{x}_p)(\pmb{x}_i-\pmb{x}_p)^T \]

which has a simple form \(\frac{1}{4}\Delta x^2\pmb{I}\) for quadratic and \(\frac{1}{3}\Delta x^2\pmb{I}\) for cubic interpolation stencils.

Then from grid to particle

\[ \begin{cases} \pmb{v}_p=\sum_i w_{ip}\pmb{v}_i\\ \pmb{B}_p=\sum_i w_{ip}\pmb{v}_i (\pmb{x}_i-\pmb{x}_p)^T \end{cases} \]

Deformation Gradient Update

Given \(f_p\), we can update the position and velocity of the grid

\[ \begin{cases} \pmb{v}_i^{n+1}=\pmb{v}_i^n+\Delta t f_i(\pmb{x}_i^n)/m_i\\ \pmb{x}_{i}^{n+1}=\pmb{x}_i^n+\Delta t \pmb{v}_i^{n+1} \end{cases} \]

Given a grid velocity field \(\pmb{v}_i^{n+1}\), we can update \(F\) as

\[ F_p^{n+1}=\left(I+\Delta t \sum_i \pmb{v}_i^{n+1}(\nabla w_{ip}^n)^T\right) F^n_p \]

Forces

MPM Forces are defined on grid nodes. If we assume a deformation gradient based hyperelastic energy density, then the total elastic potential energy is then

\[ e=\sum_p V_p^0\Psi_p(F_p) \]

where \(V_p^0\) is the material space volume of particle.

Nodal elastic force is the negative gradient of the total potential energy evaluated at nodal positions. So the MPM spatial discretization of the stress-based forces is given as

\[ f_i(x_i^n)=-\frac{\partial e}{\partial x_i}(x)=-\sum_p V_p^0\left(\frac{\partial \Psi_p}{\partial F}(F_p(x_i^n))\right)(F^n_p)^T\nabla w_{ip}^n \]

which fully depends on the existing particle/grid weights and particle attributes.