# Optimization based integration for MPM

I'm considering implementing (just for simplicity) the unconstrained implicit optimization based integration for Material Point Method as described in Chenfanfu Jiang's thesis on MPM (the minimization algorithm starts on page 101 – "A.3.1 Unconstrained minimization").

I have difficulties understanding how would be the algorithm implemented.

## Are the statements below (in)correct?

The goal of the algorithm is obtaining new grid velocities in response to external and internal forces of the material represented by material points.

The function to be minimized is $E(\mathbf{v}) = \sum_i\frac{1}{2}m_{\mathbf{i}} \lVert \mathbf v_i - \mathbf {v_i}^n \rVert + \Phi(\mathbf {x_i}^n + \Delta t\mathbf v_i)$. (from Jiang's thesis)

It's derivative is $\mathbf h(\mathbf v) = \mathbf M\mathbf v - \Delta t \mathbf f(\mathbf x^n + \Delta t \mathbf v) - \mathbf M\mathbf v^n$.

Based on the thesis I assume the approach to minimizing $E$ is by the Newton method. The Newton method requires inverse of hessian $\mathbf{\mathit{H}}$. However we're in fact interested in $\Delta \mathbf v$ that would minimize (the linearization of) $E$, so we can instead solve $\mathbf{\mathit{H}}\Delta \mathbf v=\mathbf h$.

When solving the system via Conjugate gradient, we only need to know the product $\mathbf{\mathit{H}}\Delta \mathbf v$, which can be thought of as differential $\delta \mathbf h = \mathbf M\Delta \mathbf v - \Delta t^2 \delta \mathbf f(\mathbf x^n + \Delta t \mathbf v) = \mathbf{\mathit{H}}\Delta \mathbf v$.

From [Stomakhin et al. MPM snow paper] we know how to compute $\delta\mathbf{f_i}$ (page 6 – "6 Stress-based forces and linearization") – i.e. force differential for the grid velocity at node $\mathbf i$.

My assumptions:

1. When solving $\mathbf{\mathit{H}}\Delta \mathbf v=\mathbf h$ via Conjugate Gradient (step 4. in the algorithm)
1. ${\Delta \mathbf v}$ is a column vector consisting of component-wise grid velocity differentials – $[{\Delta v}_{1x}, {\Delta v}_{1y}, {\Delta v}_{1z},\,\,\, {\Delta v}_{2x}, {\Delta v}_{2y}, {\Delta v}_{2z}, \,\,\,\ldots]$ (length of $\Delta \mathbf v$ is therefore $3n$, where $n$ is the number of grid nodes (containing velocity))
2. When multiplying $\mathbf{\mathit{H}}\Delta \mathbf v$, we're in fact always processing 3 rows at the same time - i.e. we process one grid node at time (which consists of x,y,z velocity components).
3. When computing residual $\mathbf r = \mathbf h - \mathbf{\mathit{H}}\Delta \mathbf v$, we can compute its consecutive 3 components corresponding to grid node $\mathbf i$ (like in assumption 1.1.) by writing $\mathbf{r_i} = \left(m_{\mathbf i}\mathbf v_{\mathbf i} - \Delta t \mathbf f_{\mathbf i}(\mathbf x^n_{\mathbf i} + \Delta t \mathbf v_{\mathbf i}) - m_{\mathbf i}\mathbf {v_{\mathbf i}}^n\right) - \left(m_{\mathbf i}\Delta\mathbf v_{\mathbf i} - \Delta t^2 \delta \mathbf {f_i}(\mathbf x^n_{\mathbf i} + \Delta t \mathbf v_{\mathbf i})\right)$.
4. When computing $\mathbf r \cdot \mathbf r$, I ignore the fact that each 3 consecutive elements of vector $\mathbf r$ correspond to 1 grid node, and just compute the dot product like with any other vector.