View source for MPM Methods and Simulation Timing

The [[#Input Commands|input commands]] described below set details about the MPM method and control the simulation timing. Before describing the commands, the first few sections give background theory about these options.

== Theory: Shape Functions ==

Many tasks in MPM involve extrapolations from particles to the grid or from the grid to particles. These extrapolations are controlled by "Shape Functions," and the various MPM methods correspond to various methods for evaluating the shape functions. In the most generalized description of MPM (called GIMP for General Interpolation Material Point<ref name='GIMP'>S. G. Bardenhagen and E. M. Kober, "The Generalized Interpolation Material Point Method," ''Computer Modeling in Engineering & Sciences'', '''5''', 477-496 (2004).</ref>), the shape function (S<sub>pi</sub>) and shape function gradient (G<sub>pi</sub>) for extrapolation from node i to particle p are:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>S_{pi} = {\int_{\Omega_p} \chi_p(\vec x)N_i(\vec x) dV\over \int_{\Omega_p} \chi_p(\vec x) dV}\qquad{\rm and}\qquad G_{pi} = {\int_{\Omega_p} \chi_p(\vec x)\nabla N_i(\vec x) dV\over \int_{\Omega_p} \chi_p(\vec x) dV}</math>

where &Omega;<sub>p</sub> is the domain for particle p, &chi;<sub>p</sub>(x) is the characteristic function for particle p, and N<sub>i</sub>(x) is grid shape function for node i. The various MPM methods depend on choices for &chi;<sub>p</sub>(x) and N<sub>i</sub>(x) and on the technique used for evaluating the shape function integrals. The chart on the right shows an MPM "genealogy" or shows how the various MPM methods descend from GIMP. The methods colored in blue are the ones supported in [[NairnMPM]].
[[File:MPMTree2.png|right]]

=== Particle Shape Functions ===

One choice is for  &chi;<sub>p</sub>(x) or the particle shape functions &mdash; the chart shows three options. The choice of the Dirac delta function (&delta;(x)) leads to classic MPM methods. The choice of a Heaviside function (''i.e.'', &chi;<sub>p</sub>(x) = 1 in &Omega;<sub>p</sub> and 0 elsewhere) leads to a family of Heaviside GIMP methods. Finally the choice of some other function (''e.g.'', a Gaussian centered on the particle) could lead to other MPM methods, but is currently an unexplored area.

=== Spline Grid Shape Functions ===

A second choice is for grid shape functions can be derived from B<sup>n</sup> spline basis functions, which are C<sup>(n-1)</sup> continuous spline functions on the grid. Most MPM methods choose B<sup>1</sup> order spline functions or linear shape functions.  In 1D, linear shape functions have the form:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>N_i(x) = \left\{ \begin{array}{cc} 1+{x-x_i\over \Delta x} & {\rm for\ } x_i-\Delta x < x < x_i \\
     1-{x-x_i\over \Delta x} & {\rm for\ } x_i \le x < x_i+\Delta x \\
    0 & {\rm otherwise}
\end{array} \right.</math>

where x<sub>i</sub> is coordinate of node i and &Delta;x is cell size (assumed same on both sides of node i).

Another choice supported in [[NairnMPM]] is to use B<sup>2</sup> spline functions or quadratic spline functions. In 1D, quadratic spline shape functions have the form:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>N_i(x) = \left\{ \begin{array}{cc} {1\over 8}\left(3+{2(x-x_i)\over \Delta x}\right)^2 & {\rm for\ } x_i-{3\Delta x\over 2} < x < x_i-{\Delta x\over 2} \\
      {3\over 4} - \left({x-x_i\over \Delta x}\right)^2 & {\rm for\ } x_i-{\Delta x\over 2} \le x \le x_i+{\Delta x\over 2}  \\
      {1\over 8}\left(3-{2(x-x_i)\over \Delta x}\right)^2 & {\rm for\ } x_i+{\Delta x\over 2} < x < x_i+{3\Delta x\over 2} \\
    0 & {\rm otherwise}
\end{array} \right.</math>

== Classic MPM  Methods ==

If &chi;<sub>p</sub>(x) = &delta;(x), the shape functions reduce to the grid shape function (''e.g.'', S<sub>ip</sub> = N<sub>i</sub>(x<sub>p</sub>)), which recovers "Classic MPM" or the original derivation of MPM.<ref name='SS'>D. Sulsky, Z. Chen, and H. L. Schreyer, "A Particle Method for History-Dependent Materials," ''Comput. Methods Appl. Mech. Engrg.'', '''118''', 179-186 (1994).</ref><ref name='SZS'>D. Sulsky, S. -J. Zhou, and H. L. Schreyer, "Application of a Particle-in-Cell Method to Solid Mechanics," ''Comput. Phys. Commun.'', '''87''', 236-252 (1995).</ref><ref>S. -J. Zhou, "The Numerical Prediction of Material Failure Based on the Material Point Method," Ph.D. Thesis, University of Mexico (1998).</ref> (when using linear N<sub>i</sub>(x)).  Note that "Classic GIMP" predates GIMP even though it descends from GIMP in the MPM "genealogy." The "Dirac" method options are

* <tt>Classic</tt> (or <tt>Dirac</tt>) chooses "Classic MPM" with linear N<sub>i</sub>(x) on the grid.
* <tt>B2SPLINE</tt> chooses Dirac particle functions with quadratic spline N<sub>i</sub>(x) on the grid (do not use yet for axisymmetric MPM).

Although Classic MPM is the most efficient method, it is rarely used because it is not very robust for large displacement problems. Whenever particles cross cell boundaries, they can cause noise that soon leads to unsatisfactory results. These cell-crossing artifacts are reduced by using <tt>B2SPLINE</tt> (because first derivatives are continuous), but problems often still persist. The better approach to handling cell cross artifacts in to choose a different &chi;<sub>p</sub>(x) (''e.g.'', the Heaviside methods).

== Heaviside GIMP Methods ==

MPM results are greatly improved by choosing &chi;<sub>p</sub>(x) to be a Heaviside function. With this function, the shape function integrals reduce to integrals over the particle domain of the grid-based shape functions, which are normalized to particle volume (see chart). Many MPM methods descend from this selection. The specific descendant methods  depend on how they integrate over the particle domain.

=== Undeformed GIMP ===

The simplest approach is to assume the particle domain has the same shape as the initial domain and merely translates with particle motion. This method, which is called "Undeformed GIMP", has the greatest efficiency because shape functions can be explicitly calculated for all possible particle locations relative to a regular grid. Undeformed GIMP works well from small to fairly large deformation and greatly reduces cell-crossing artifacts. The two undeformed GIMP methods can be selected using:

* <tt>uGIMP</tt> or undeformed GIMP with linear N<sub>i</sub>(x) on the grid.
* <tt>B2GIMP</tt> or undefomed GIMP with quadratic spline N<sub>i</sub>(x) on the grid (do not use yet for axisymmetric MPM).
 
Very large tension strains (>50%), however, can cause undeformed particle domains to become numerically separated, which can lead to non-physical fracture of the object. When this occurs, one option is to use a different method that accounts for deformed particle domain.

=== Deformed GIMP ===

The following GIMP methods attempt to account for particle deformation in evaluation of the shape function integrals. The methods differ in how they do the domain integration and in the grid shape functions that are integrated.

==== 1. Finite GIMP ====

In finite GIMP, the shape function integrals are fully integrated over the deformed particle domain. This option can be selected using <tt>Finite</tt> for shape function method. This exact integration is available only in [[OSParticulas]], only for 2D calculations, and only for linear grid shape functions. The full integration compared to approximate integral is similar in efficiency, but surprisingly, the results are not improved over other methods. The reasons are not yet clear. One hypothesis is that while finite GIMP exactly integrates the full particle domain, the particle shape function is still treating particle properties (such as velocity) as constant over the domain (due to Heaviside &chi;<sub>p</sub>(x) function). These contrary views may be inconsistent.

==== 2. Convected Partical Domain Integration ====

A useful improvement for large to massive deformation is Convected Particle Domain Integration approach or CPDI.<ref name='CPDI'>A. Sadeghirad, R. M. Brannon, and J. Burghardt, "A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations," <i>Int. J. Numer. Meth. Engng.</i>, '''86''', 1435–1456 (2011)</ref> In CPDI, the integral over the deformed particle domain is approximated. In brief:

<ol>
<li>The deformed particle domain is discretized by isoparametric finite elements. In practice, the discretization uses a single element, but it could use more.</li>
<li>The Integrands for S<sub>pi</sub> and G<sub>pi</sub> are then approximated by expanding them in the finite element shape functions:<br>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>N_i(\vec x) \approx \sum_j N_i(c_j) M_j(\vec x)\qquad{\rm and}\qquad \nabla N_i(\vec x) \approx \sum_j N_i(c_j) \nabla M_j(\vec x)</math><br>
where c<sub>j</sub> are nodes of the deformed particle domain (''e.g.'', usually corners of the domain) and M<sub>j</sub>(x) are isoparametric finite element shape functions for the elements used to discretize the deformed particle domain.</li>
<li> Integrals over the deformed domain can then be evaluated to give S<sub>pi</sub> and G<sub>pi</sub> in terms for N<sub>i</sub>(c<sub>j</sub>).</li>
</ol>

The version of CPDI available here finds the corners of the deformed domain using the deformation gradient on the particle.<ref name='CPDI'/> If r<sub>j0</sub> is a vector from the particle to any edge location j in the undeformed domain, then r<sub>j</sub> = F<sub>p</sub>r<sub>j0</sub>, where F<sub>p</sub> is the deformation gradient on the particle, is the vector to that location in the deformed domain. When using this approach, an initial 2D rectangle domain deforms to a parallelogram and an initial 3D orthogonal box deforms to a parallelepiped. The subclasses of this method depend on grid shape functions being used and on the finite element types used to map that particle domain:

* <tt>lCPDI</tt> for "linear" CPDI: this method uses linear grid shape functions and maps the corners of the domain to a single linear element (4 node element in 2D and 8 node element in 3D) resulting in an isoparametric element with linear shape functions.
* <tt>qCPDI</tt> for "quadratic" CPDI: this method uses linear grid shape functions and additionally maps midpoints of edges to nodes (''e.g''., 8 node element in 2D) resulting in an isoparametric element with quadratic shape functions. This method is only available for 2D simulations.
* <tt>B2CPDI</tt> for "B2 Spline" CPDI: this method uses B<sup>2</sup> spline grid shapes functions and maps the deformed domain as in lCPDI.

CPDI methods are less efficient the undeformed GIMP because of the extra overhead required to locate corner nodes needed for evaluating shape functions. Nevertheless, the use of CPDI is crucial for handling problems with deformation beyond the capabilities of undeformed GIMP.

==== 3. Convected Particle Domain Integration 2 ====

A second version of CPDI, known as CPDI2, has been developed.<ref name='CPDIQ'>A. Sadeghirad, R. M. Brannon,and J. E. Guilkey, "Second-order convected particle domain interpolation (CPDI2) with enrichment for weak discontinuities at material interfaces," <i>Int. J. Numer. Meth. Engng.</i>, '''95''', 928–952 (2013)</ref> In CPDI2 the corner locations are separately tracked during the analysis and their current locations are used to define the deformed particle domain. For example, a 2D domain deforms to an arbitrary quadrilateral element and a 3D domain to an arbitrary 8-cornered brick element. The efficiency should be similar to CPDI, but CPDI2 adds the additional overheard of tracking and updating corner particles. CPDI2 options are currently not implemented in either  [[NairnMPM]] or [[OSParticulas]].

==== 4. cpGIMP ====

In cpGIMP, or convected particle GIMP, the deformed domain is stretched in orthogonal directions but is never sheared. In other words, particle domains that start as squares (in 2D) or cubes (in 3D) will deform to rectangles (in 2D) or orthogonal boxes (in 3D). In 1D, cpGIMP is a true [[#Finite GIMP|finite GIMP]] method, but in 2D or 3D, it does not account for arbitrary deformation that includes shear deformation. Compared to undeformed GIMP, cpGIMP may handle uniaxial deformation better, but may not be any better for general deformation modes. This approach is currently not implemented in either  [[NairnMPM]] or [[OSParticulas]].

== Theory: Stress and Strain Updates ==

The basic tasks in an MPM time step are:

# Extrapolate particle mass and momentum to the grid
# Find internal and external forces
# Update momenta on the grid using those forces
# Update particle position and velocity using updated grid information

To correctly model internal forces, this loop needs to track stresses and strains on the particles. This tracking is done by calculating velocities on the grid (from the grid momenta). Once these velocities are found, they can be extrapolated to the particle using gradient shape functions to get velocity gradients on the particle. Particle velocity gradients provide the key input for updating particle strains and through constitutive laws to update particle stresses. In the first MPM work,<ref name='SS'/> the strain and stress update was done using the updated momenta on the grid. It is easy to show, however, that this approach is unstable if the update particle position interacts with a node having zero mass. This numerical challenge can be eliminated by two variants to MPM:

* <tt>USL</tt> - follow the recommendation of Sulsky, Zhou, and Schreyer<ref name='SZS'/> which was to do a second extrapolation of momenta after updating momenta on the grid and use those new momenta to update strains and stresses. This approach adds step 5 to MPM time step above to reextrapolate momenta to the grid and then update stress and strain after that step.
* <tt>USF</tt> - update strains and stresses first or to update them after the initial extrapolation in step 1 of above time step tasks.<ref name="USF'>S. G. Bardenhagen, "Energy Conservation Error in the Material Point Method," <i>J. Comp. Phys.</i>, <b>180</b>, 383-403 (2002).</ref>

These two approaches give similar results. The <tt>USF</tt> method is more efficient than <tt>USL</tt> because it avoids the need to do a second extrapolation of momenta.

A unique option in [[NairnMPM]] combines these two methods:

* <tt>USAVG</tt> - update strains after the first extrapolation and after a second extrapolation of updated momenta on the grid.

The <tt>USAVG</tt> method seems to behave like a midpoint rule. In observations of temporal convergence, the <tt>USAVG</tt> method for a given time step gives results close to <tt>USF</tt> or <tt>USL</tt> using half that time step. In other words,  compared to <tt>USF</tt> and <tt>USL</tt>, the <tt>USAVG</tt> is analogous to using half the time step without the cost of doubling the number of calculations. In most simulation, the <tt>USAVG</tt> method also exhibits improved energy conservation<ref name='CRAMP'>J. A. Nairn, "Material Point Method Calculations with Explicit Cracks," <i>Computer Modeling in Engineering &amp; Sciences</i>, <b>4</b>, 649-664 (2003). ([http://www.cof.orst.edu/cof/wse/faculty/Nairn/papers/MPMCracks.pdf See PDF])</ref>.

Although  Sulsky, Zhou, and Schreyer<ref name='SZS'/> recommend a second momentum extrapolation prior to update, that second extrapolation is not needed provided one is careful in calculation of shape functions. Skipping the second extrapolation is always more efficient and can result in better simulations. Whether or not to use the second extrapolation, however, seems to depend on the problem. Some problems with contact or cracks seem to require it while problems with liquids seem improved without it. The default is to include the second extrapolation, but [[#Input Commands|commands below]] let you pick from the following options:

* <tt>USL+</tt> - update strains after momentum update and after a second extrapolation (same as <tt>USL</tt> above).
* <tt>USL-</tt> - update strains after momentum update, but do not do a second extrapolation.
* <tt>USAVG+</tt> - update strains after the first extrapolation and after a second extrapolation of updated momenta on the grid (same as <tt>USAVG</tt> above).
* <tt>USAVG-</tt> - update strains before and after the momentum update, but do not do a second extrapolation.

The key to skipping the second extrapolation compared to early MPM work,<ref name='SS'/> is to extrapolate velocity gradient from updated momenta by using gradient shape function based on the original particle position rather then the updated position.

=== Incremental Deformation Gradient ===

During strain updating, some [[Material Models|materials]] (''e.g.'', hyperelastic materials) use the full incremental deformation gradient, ''dF'', which is defined by

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>{dF\over dt} = \nabla\vec v F</math>

where &Delta;t is the time step and &nabla;v is the velocity gradient calculated from the current nodal velocities. If velocity gradient is constant for the time step, the exact solution is

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math> F(t+\Delta t) = \exp(\Delta t \nabla\vec v) F(t) = dF\cdot F(t)</math>

where

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>dF = \exp(\Delta t \nabla\vec v)</math>

is the incremental deformation gradient. An interesting review article<ref name='expM'>C. Moler and C. Van Loan, "Nineteen Dubious Ways to Compute the Exponential of a Matrix. Twenty Five Years Later." ''SIAM Review'', '''46''', (2003)</ref> on finding the exponential of a matrix, which is needed here, concludes no single way is stable and efficient for all problems, but in MPM with sufficiently small time steps, an expansion method should work well

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\exp(\Delta t \nabla\vec v) = \sum_{k=0}^{k_{max}} {(\nabla u)^k \over k!} 
            =  I + \nabla u + \sum_{k=2}^{k_{max}} {(\nabla u)^k \over k!}</math>

where &nabla;u = &Delta;t &nabla;v is the incremental displacement gradient for the current time step. Many computational mechanics methods truncate after the first term (the linear term), because it avoids the need to take powers of a matrix. But, using just the linear term can have issues for problems with lots of rotational motion. Fortunately the extra terms can all be evaluated with no matrix multiplications (in 2D) and with only the square of a matrix (in 3D) by using the "Cayley-Hamilton" theory (see method #8 in the review article<ref name='expM'/>). You can select the number of terms to use in the expansion. The default is 2, which seems to work well for almost all simulations and to reduce errors caused by large rotations.

== Theory: MPM Time Step ==

In explicit computation mechanics, the Courant–Friedrichs–Lewy condition (CFL condition) is a necessary condition for convergence. It states that the time step, &Delta;''t'', must satisfy:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\Delta t \le { \Delta x\over v}</math>

where &Delta;''x'' is the minimum dimension of any cell in the background grid and ''v'' is the maximum wave speed in the material (''i.e.'', the highest wave speed in anisotropic materials). It is usually better to be below this minimum by a input factor of less than one. The factor is also known as the "Courant" factor:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math> C = { v\Delta t\over \Delta x} \qquad {\rm or} \qquad\Delta t =C { \Delta x\over v}</math>

The entered ''C'' must be less then 1.0. The default in [[NairnMPM]] is ''C'' = 0.5. Lower values may be are needed to resolve instabilities around [[Setting Velocity and Transport Values|fixed-velocity boundary conditions]] on the grid or created by [[Rigid Material|rigid particles]] when modeling [[Multimaterial MPM|contact]], and when modeling brittle materials using [[Material Models#Softening Materials|damage mechanics]]. An alternative to very low ''C'' is to try [[Restarting Time Steps|restarting time steps]] that develop high accelerations on the grid.

When transport tasks are activated ([[Thermal Calculations|conduction]], [[Diffusion Calculations|diffusion]], or [[Poroelasticity Calculations|poroelasticity]]), convergence of transport calculations requires the time step to satisfy

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\Delta t_T \le { \Delta x^2\over 2}{C_T\over \max(\kappa)}</math>

where ''C<sub>T</sub>'' and <math>\kappa</math> depend on the type of transport being modeled. For [[Thermal Calculations|conduction]], <math>C_T=\rho C_v</math> and <math>\kappa</math> is thermal conductivity. For [[Diffusion Calculations|diffusion]], ''C<sub>T</sub>=1'' and <math>\kappa</math> is diffusion coefficient. For [[Poroelasticity Calculations|poroelasticity]], ''C<sub>T</sub>=1/Q'' and <math>\kappa</math> is Darcy's law permittivity divided by fluid viscosity. It is usually better to be below this minimum by a input factor of less than one. The factor is also known as the "Transport Courant" factor or transport time step is:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\Delta t_T = F { \Delta x^2\over 2}{C_T\over \max(\kappa)}</math>

The entered ''F'' must be less then 1.0. The default in [[NairnMPM]] is ''F'' = 0.5.

When running coupled calculations, the simulation time step must use the minimum of mechanics and transport time steps or:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>\Delta t_{final} = \min(\Delta t,\Delta t_T)</math>

For realistic properties, the minimum is usually <math>\Delta t</math> by several orders of magnitude, which means entered value for ''F'' has no effect. Notice, however, that transport time step scales with <math>\Delta x^2</math> while mechanics time step is linear in <math>\Delta x</math>. The minimum value could change if using very small cells.

=== Restarting Time Steps ===

Some problems might develop high accelerations on the grid that cause simulations to break down (such as particles leaving the grid). The first steps to fixing such problems should be:

# Use a smaller [[#Theory: MPM Time Step|time step ''C'' factor]]
# Add an [[AdjustTimeStep Custom Task]] - this option only works for materials that support it
# Increase number of terms used to find the [[#Incremental Deformation Gradient|incremental deformation gradient]] - this option works for all materials except [[Material Models|small-strain materials]] when they are not using the [[Common Material Properties#Basic Properties|largeRotation=1]] option.

If none of these work, you can restart time steps as needed based on nodal velocities and accelerations. Each MPM time step extrapolates mass (''m''), momentum (''p''), and force (''f'') to the grid. When time-step restarting is activated, each step checks the virtual nodal displacement on each node implied by extrapolated:

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
<math>d = \| \vec p + {1\over2}\vec f\Delta t\|{\Delta t\over m}</math>

where &Delta;''t'' is the current time step. If any node's ''d'' relative to its cell size is greater than an entered rescaling factor ''S'', the time step is discarded and restarted using a smaller time step. This option is invoked using the [[#Input Commands|<code>&lt;RestartScaling&gt;</code> command]]. To prevent the time step from permanently remaining small (and make a simulaiton too slow), this option should always be used in conjunction with an [[AdjustTimeStep Custom Task]].

== Input Commands ==

The [[#Theory: Shape Functions|MPM shape functions]], the [[#Theory: Stress and Strain Updates|strain updating method]], the number of points per element, and [[#Theory: MPM Time Step|simulation timing]] are set with one or more of the commands in this section. In script input files, the commands are:

 MPMMethod (strainUpdate),<(MPM_method)>
 PtsPerElement (axisNumber)
 TimeStep (timeStep),<(maxTime)>,<(timeFactor)>
 MaximumTime (maxTime)
 CFLFactor (timeFactor),<(transTimeFactor)>
 CPDIrcrit (rcrit)

In <tt>XML</tt> input files, the settings are made with one or more of the following commands (all of which must be within the [[MPM Input Files#MPM Header|<tt><MPMHeader></tt> element]]):

 <MPMMethod>(strainUpdate number)</MPMMethod>
 <SkipPostExtrapolation/>
 <GIMP type='(MPM_method)'/>
 <CPDIrcrit>(rcrit)</CPDIrcrit>
 <MatlPtsPerElement>(totalNumber)</MatlPtsPerElement>
 <TimeStep>(timeStep)</TimeStep>
 <TimeFactor>(timeFactor)</TimeFactor>
 <TransTimeFactor>(transTimeFactor)</TimeFactor>
 <MaxTime>(maxTime)</MaxTime>
 <DefGradTerms>(kmax)</DefGradTerms>

where

* <tt>(strainUpdate)</tt> - the method used to [[#Theory: Stress and Strain Updates|update stresses and strains]] on the particles. The options are (scripted files can use name or number while <tt>XML</tt> files must use number):
** <tt>USF</tt> (or 0) - update first
** <tt>USAVG</tt> or <tt>USAVG+</tt> (or 2) - update strain first and last. This option is the default option.
** <tt>USAVG-</tt>  (no number) - update strain first and last, but last is [[#Notes|modified<sup>1</sup>]]. Because there is no number option, <tt>XML</tt> input files should add the <tt><SkipPostExtraolation></tt> command to get this method.
** <tt>USL</tt> or <tt>USL+</tt> (or 3) - update strain last (<tt>SZS</tt> is old name for this method that is still accepted)
** <tt>USL-</tt>  (no number) - update strain last, using [[#Notes|modified<sup>1</sup>]] method. Because there is no number option, <tt>XML</tt> input files should add the <tt><SkipPostExtraolation></tt> command to get this method.
* <tt>(MPM_method)</tt> - the [[#Theory: Shape Functions|MPM method]] used for finding shape functions. The options are:
** <tt>Classic</tt> (or <tt>Dirac</tt> or 0) - classic MPM. This method is the default if no option is specified or if no <tt><GIMP></tt> command is used in an <tt>XML</tt> input file.
** <tt>uGIMP</tt> (or <tt>GIMP</tt> or 1)[[#Notes|<sup>2,3</sup>]] - undeformed GIMP. If <tt><GIMP></tt> command is used in <tt>XML</tt> files without a <tt>type</tt> attribute, this option is used.
** <tt>lCPDI</tt> (or <tt>CPDI</tt> or 2)[[#Notes|<sup>3</sup>]] - CPDI using linear shape functions.
** <tt>qCPDI</tt> (or 3)[[#Notes|<sup>3</sup>]] - CPDI using quadratic shape functions (available in 2D only).[[#Notes|<sup>5</sup>]]
** <tt>Finite</tt> (or 4)[[#Notes|<sup>2,3</sup>]] - exactly integrated finite GIMP (available in 2D only and only in [[OSParticulas]]).
** <tt>B2GIMP</tt> (or 5)[[#Notes|<sup>2,4</sup>]] - undeformed GIMP methods using quadratic spline grid shape functions.
** <tt>B2SPLINE</tt> (or 6)[[#Notes|<sup>2,4</sup>]] - Classic method using quadratic spline grid shape functions.
** <tt>B2CPDI</tt> (or 7)[[#Notes|<sup>2,4</sup>]] - CPDI methods using quadratic spline grid shape functions.
* <tt>(axisNumber)</tt> (scripted files only) - the number of material points along each axis in the grid.[[#Notes|<sup>6</sup>]] The options are
** 1 to 5 for 2D - these settings results in 1, 4, 9, 16, or 25 points per element. The default is 2 or 4 points per element.
** 1 to 3 for 3D - these settings results in 1, 8 or 27 points per element. The default is 2 or 8 points per element.
* <tt>(rcrit)</tt> - one problem with CPDI shape functions is that if a particle deforms too much, its' corners may extend beyond ghost rows into other patches used in parallel coding. If that becomes a problem, this parameter will restrict corner deformation to keep them within a single patch. Enter maximum allowed particle radius in units of cell size. Truncation occurs when half of any diagonal of the deformed particle exceeds <tt>(rcrit)</tt>. The particle domain truncation implements the methods described in Homel ''et al.''<ref>Homel, M. A., Brannon, R. M., and Guilkey, J. (2016). Controlling the onset of numerical fracture in parallelized implementations of the material point method (MPM) with convective particle domain interpolation (CPDI) domain scaling. ''Int. J. Numer. Meth. Engng'', '''107''':31–48.</ref> (see Eqs. (4) to (11)).
* <tt>(totalNumber)</tt> (<tt>XML</tt> files only) - the number of material points in each element of the background grid.[[#Notes|<sup>6</sup>]] The options are
** 1, 4, 9, 16, or 25 for 2D. The default is 4 points per element.
** 1, 8 or 27  for 3D. The default is 8 points per element.
* <tt>(timeStep)</tt> - the MPM time step in [[ConsistentUnits Command#Legacy and Consistent Units|alt time units]] (in <tt>XML</tt> files, the time is in [[ConsistentUnits Command#Legacy and Consistent Units|time units]] or determined by a [[Units Attribute|units attribute]]). This input will be ignored unless it is less that time step calculated for input (or default) values for <tt>(timeFactor)</tt> and <tt>(transTimeFactor)</tt>.
* <tt>(maxTime)</tt> - the time that the calculations will stop in [[ConsistentUnits Command#Legacy and Consistent Units|alt time units]] (in <tt>XML</tt> files, the time is in [[ConsistentUnits Command#Legacy and Consistent Units|time units]] or determined by a [[Units Attribute|units attribute]]). In scripted files, this value can be entered as second argument to the <tt>TimeStep</tt> command or in a separate <tt>MaximumTime</tt> command (whichever is last will be the one that is used) (it can be an [[Entity Command|entity]]).
* <tt>(timeFactor)</tt> and <tt>(transTimeFactor)</tt> - the  [[#Theory: MPM Time Step|Courant factors]] to use when determining the MPM time step. The time step used in the analysis will be the minimum of the one determined using these factors and the one entered in <tt>(timeStep)</tt>. The default [[#Theory: MPM Time Step|Courant factors]] are 0.5. In scripted files, The <tt>(timeFactor)</tt> can be entered as third argument to the <tt>TimeStep</tt> command or first argument in a <tt>CFLFactor</tt> command
* <tt>(kmax)</tt> - the number of terms to use when evaluating the [[#Incremental Deformation Gradient|incremental deformation gradient]].[[#Notes|<sup>7</sup>]] The default is 2.

== Notes ==

<ol>
<li>The <tt>USAVG+</tt> and <tt>USL+</tt> methods re-extrapolate momenta to the grid and re-impose contact conditions prior to updating the strain after the momenta update. The <tt>USAVG-</tt> and <tt>USL-</tt> methods skip the second extrapolation. The method to use may depend on the problem, but some problems with contact seem to require the "+" methods. If both are equal, the "-" methods are more efficient.
<li><tt>uGIMP</tt>, <tt>Finite</tt>, <tt>B2GIMP</tt>, <tt>B2SPLINE</tt>, and <tt>B2CPDI</tt> all require use of a regular grid (''i.e.'', a grid with constant size and orthogonal elements).
<li>These methods are based on linear N<sub>i</sub>(x) grid shape functions.
<li>These methods are based on quadratic spline N<sub>i</sub>(x) grid shape functions.
<li><tt>qCPDI</tt> is less efficient than <tt>lCPDI</tt> and in trial runs appears to have little or no benefit. It is available for comparison and potential future research into MPM shape function methods. 
<li>The number of points per element must be set before defining material points in your analysis. Once you start defining material points, it cannot be changed. In other words elements with different numbers of material points cannot be used in a single analysis. Note that script files and <tt>XML</tt> files set this number differently. Script files set the number for each axis (which is squared for 2D or cubed for 3D) while <tt>XML</tt> files gives the total number for the cells.
<li>The default number of terms for the incremental deformation gradient is 2 for 2D simulations and 1 for 3D simulations (the latter for efficieny). Because these settings work well for most calculations, there is currently no script command to change this setting. If necessary, you can change it using an [[XMLData Command]]:
<pre>XMLData "MPMHeader"
    <DefGradTerms>(kmax)</DefGradTerms>
EndXMLData
</pre>
</ol>

== References ==
<references/>