On the use of domain-based material point methods for problems involving large distortion

Challenging solid mechanics problems exist in areas such as geotechnical and biomedical engineering which require numerical methods that can cope with very large deformations, both stretches and torsion. One candidate for these problems is the Material Point Method (MPM), and to deal with stability issues the standard form of the MPM has been developed into new domain-based techniques which change how information is mapped between the computational mesh and the material points. The latest of these developments are the Convected Particle Domain Interpolation (CPDI) approaches. When these are demonstrated, they are typically tested on problems involving large stretch but little torsion and if these MPMs are to be useful for the challenging problems mentioned above, it is important that their capabilities and shortcomings are clear. Here we present a study of the behaviour of some of these MPMs for modelling problems involving large elasto-plastic deformation including distortion. This is carried out in a unified implicit quasi-static computational framework and finds that domain distortion with the CPDI2 approaches affects some solutions and there is a particular issue with one approach. The older CPDI1 approach and the standard MPM however produce physically realistic results. The primary aim of this paper is to raise awareness of the capabilities or otherwise of these domain-based MPMs.


Introduction
The Material Point Method (MPM) [1][2][3], originally proposed as an extension of a similar method known as the Fluid-Implicit Particle or FLIP method [4], itself an extension of the Particle-in-Cell method [5], is a numerical method combining advantages of Eulerian and Lagrangian approaches to solve solid mechanics problems. In the MPM, a body is described by a number of material points which store physical field variables such as stiffness, stress and displacement. At the beginning of each load step, the information held at material points is mapped to nodes  Fig. 1 for a simple shear problem, the total deformation and other state variables are stored at the material points, while the background mesh is reset and extended with the incremental displacement, thus avoiding the mesh distortion seen with the standard Finite Element Method (FEM) for large deformation problems. In fact, the background mesh can be any mesh at the beginning of each load step. Because of this attraction the MPM has been applied to many large deformation problems, e.g. [6][7][8][9][10][11][12][13][14][15].
The method of mapping information between the material points and background mesh nodes can significantly influence the stability of the MPM. In the standard MPM (sMPM) [1], a material point only influences (and is only influenced by) its parent element, (i.e. the background element in which it is currently located), and the conventional FEM shape functions are used to map information between the nodes and the material points. However, the sMPM has instability issues from various sources: (i) the transition of a material point between elements leads to a sudden change in stiffness and internal force; (ii) under deformation fields involving large stretch, material points can end up separated by more than one element, resulting in artificial fracturing of the physical domain [16]; and (iii) partially filled elements with material points close to an element boundary can lead to an ill-conditioned global stiffness matrix. In order to reduce these problems, several extensions to the sMPM have been proposed. Each of the extensions replaces the discrete material point with a deformable material point domain, called a particle domain in some of the approaches [16][17][18]. The basis functions of these domain-based methods are based on an approximation of the integral of the background finite element shape functions over the material point's domain. The most notable of these approaches are the Generalised Interpolation Material Point (GIMP) method [19] and the Convected Particle Domain Interpolation (CPDI) approaches. The latter comprise CPDI1 [17], Second-order CPDI with quadrilateral particle domains (CPDI2q) [16] and Second-order CPDI with triangular particle domains (CPDI2t) [18]. Using a deformable particle domain for a material point enables the influence domain of a material point to overlap more than one element, thus reducing or eliminating the problems cited above.
Today, there is considerable interest in the potential use of these MPMs on very challenging problems in biomechanical and geotechnical engineering which involve large deformation and material nonlinearity. An early mechanobiological application can be found in [20], where a cell scaffold consisting of growing microvascular fragments embedded in a collagen gel is modelled using a variant of the sMPM with over 13.6 million material points. Much more recently, a review also identifies the MPM as of substantial interest for biomechanical modelling [21] and in [22], the MPM is used to model an entire human head. The standard MPM and CPDI2t methods have also been used to simulate thin-walled tubes under lateral compression and validated against experiments in [23]. The motivation for our contribution here arises from a project in geotechnical engineering which aims to model the soil response to the installation of a screwed-in pile foundation, as the first step to providing a computer-aided design tool for engineers to optimise pile design for offshore wind turbine foundations [24,25]. During installation of a screw pile, it is pushed and rotated into the ground, deforming the soil in complex patterns, a problem for which the MPM seems ideal.
The nature of the majority of these problems prompts the use of unstructured representations of the problem domain and for the background computational mesh, and of the domain-based MPMs only the CPDI approaches can support an unstructured background mesh. CPDI2 approaches can have unstructured particle domains that can more accurately discretise a general complex body shape but the remaining methods will, in general, result in gaps or overlaps between particle domains. Previously, each of the new developments of the sMPM have been demonstrated in the individual papers referenced above and have been verified on selected problems. The purpose of this contribution is to present an investigation of the behaviours of some of the domain-based MPMs on a selection of specific problems that test their predictive abilities when modelling problems involving large stretch, shear and torsion, thus hopefully providing useful guidance to those wishing to employ the methods on real-world applications of the types mentioned above. This paper also casts all methods within a common implicit quasi-static elasto-plastic computational framework so that any differences in the results are purely due to differences in their basis functions and domain updating methods under large deformation.
The paper is organised as follows. In Section 2, the formulations for the continuum problem and the basic MPM discretisation are introduced, including finite stain theory for characterising the elasto-plastic deformation, the basis functions for the different MPM approaches, a formulation for implicitly solving the global system of nonlinear equations and details of the moving mesh concept. In Section 3, the implicit computational framework is presented with details of the essential differences between the different methods explained. The framework is verified and used to investigate the methods' capacities in various examples with representative deformation fields in Section 4, including a discussion of our findings via these examples in Section 4.7. Conclusions are drawn in Section 5. In this paper we use a combination of tensor and matrix notation in order to ensure that both the continuum formulation and the steps to numerical implementation are as clear as possible.

Material point continuum formulation
This section details the quasi-static implicit finite deformation elasto-plastic material point method formulation adopted in this paper. The approach is largely based on the elasto-plastic updated Lagrangian material point formulation of Charlton et al. [26] and [27] but with the extension to non-uniform background meshes and domain-based material point approaches.

Problem statement and kinematics
The weak statement of equilibrium for the adopted formulation can be expressed in the updated configuration as ∫ where ϕ(Ω ) indicates the deformed domain of the material body, Ω , which is subjected to tractions, t i , on the boundary (with surface, s), ∂Ω , and body forces, b i , acting over the volume, v. These external loads result in a Cauchy stress field, σ i j , through the body. The weak form is derived in the current frame assuming a field of admissible virtual displacements, η i . The deformation gradient, F i j , then provides the fundamental link between the original and deformed configurations where X i are the original (or reference) coordinates and x i = ϕ(X i , t) are the updated coordinates in the current (or deformed) body. Following the work of [28][29][30], the deformation gradient is multiplicatively decomposed into elastic and plastic components where the superscripts e and p denote elastic and plastic components, respectively. This multiplicative decomposition of the deformation gradient is combined with a logarithmic strain-Kirchhoff stress formulation. This is a powerful combination as it allows existing small strain constitutive formulations to be used directly rather than reformulating them for the particular choice of stress and strain measures used in the large deformation formulation. The elastic logarithmic strain is defined as is the left elastic Cauchy-Green strain. The Kirchhoff stress, τ i j is linked to the Cauchy stress, σ i j , through where J = det(F i j ) is the volume ratio between the deformed and reference states. It is assumed that the Kirchhoff stress is linearly dependent on the logarithmic strain through where D e i jkl is the isotropic linear elastic stiffness matrix. In this work, this linear relationship between Kirchhoff stress and logarithmic strain is combined with an implicit elastic predictor-plastic corrector constitutive algorithm for elasto-plastic material behaviour (see [31], amongst others, for details of these algorithms).

Discrete material point implementation
Starting from (1), the Galerkin form of the weak statement of equilibrium for a single element in the background mesh can be expressed as where E indicates an element in the background mesh, [G] is the strain-displacement matrix containing derivatives of the basis functions with respect to the updated coordinates, and [S] is the matrix form of the basis functions. The first term in (7) is the internal force within an element and the combination of the second and third terms is the external force vector. The residual for each element, { f R E }, can be assembled into a global residual, { f R }, which is non-linear in terms of the unknown nodal displacements, {d}. One option for solving the assembled global system of equations is the standard Newton-Raphson (NR) procedure. The nodal displacements within a load step, {∆d}, can be obtained by iteratively updating the nodal displacements until (7) is satisfied within a given tolerance using where k is the current iteration within the loadstep, [K ] is the global stiffness matrix and {δd} k is the iterative increment in the displacements, { f R } k−1 is the global residual out of balance force vector associated with the previous displacement value. The current displacement in a loadstep can be obtained by summing the iterative increments within the loadstep, that is {∆d} = ∑ {δd} k . In material point methods it is more convenient to express the global equilibrium equation (7) in terms of material point, rather than element, contributions. In this context, the global residual vector can be assembled through where { f int } and { f ext } are internal and external forces. The internal force vector where V p is the current volume associated with the material point and A is the standard assembly operator. The external force vector, which is constant over the load step, is given by The global stiffness matrix, [K ], can be obtained by linearising (9) with respect to the unknown nodal displacements. This gives the stiffness associated with a single material point as For the standard MPM, V p = J V 0 p , where V 0 p is the original volume associated with the material point. For domain-based MPMs, V p is computed as the volume of the updated domain associated with material point p.
In (12), [A] is the spatial consistent tangent modulus, which can be most conveniently (and compactly) expressed in tensor notation as where and D alg i jmn is the small-strain algorithmic tangent obtained from the constitutive model. 1 An algorithm for the derivative of the logarithm of the elastic Cauchy-Green strain tensor with respect to its argument is given in [33]. The map between tensor and matrix forms can be found in many textbooks, e.g. [34].

Pseudo-time discretisation
To advance the non-linear solution algorithm, the finite deformation equations are discretised in pseudo-time by imposing the deformation over a number of load (or pseudo-time) steps. This allows the current deformation gradient to be defined using where [∆F] is the increment in the deformation gradient over the current loadstep and [F n ] is the deformation gradient from the previously converged (or initial) state. Following the work of Charlton et al. [26], the increment in the deformation gradient component is obtained from where [I ] is a three by three identity matrix, {∆u v } is the displacement increment of a background grid node (or vertex) within the current loadstep, {x n } are the coordinates at the start of the load step and n in is the number of nodes that influence the material point. 2 In order to obtain the updated stress state for the current deformation gradient, the adopted constitutive algorithm requires an initial estimate (or trial) of the elastic strain state. In this approach the trial elastic Cauchy-Green strain and logarithmic strain tensors are given by where the subscript t denotes a quantity defined in the trial state. The constitutive model should then return the updated Kirchhoff stress, [τ ], and the associated algorithmic consistent tangent, [D alg ]. Once equilibrium has been obtained over the current pseudo-time step, the material point positions are updated through where {x p } and {x p n } are coordinates of the material point, p, after and before updating, respectively.

Basis functions
The major difference between the different material point methods used in this paper (sMPM, CPDI1, CPDI2q and CPDI2t) is the basis functions, S vp , and their spatial gradients, ∇ S vp . In the following, the basis functions for the four methods are provided in a common format: • The basis functions and their spatial derivatives for the standard MPM (sMPM) are and where S v is the standard FEM basis function of node v. In the case of regular grids, it is straightforward to express the basis functions in terms of the global coordinates of the background element nodes and the material points. In general, S vp and ∇ S vp can be computed after determining the material point's local coordinates within its parent element.
• For the CPDI1 method [17], the basis functions and their spatial derivatives can be expressed as where {x c }, c = 1, 2, 3, 4 are coordinates of particle domain corners and {s} and {t} are two vectors forming the parallelogram particle domain, as shown in Fig. 2. The coordinates of the corners of the domain can be expressed in terms of the material point position, {x p } and the parallelogram vectors as • For the CPDI2q method [16] the basis functions and their spatial derivatives can be expressed as where a = ( , and x c , y c are components of coordinate {x c } for the four corners of the quadrilateral particle domain. It is obvious that both S vp and ∇ S vp only depend on the coordinates of the corners, rather than the coordinates of the material point. Therefore, in this method the positions of material points do not need to be stored.
• For the CPDI2t method [18] the basis functions and their spatial derivatives can be expressed as and As with the CPDI2q method, S vp and ∇ S vp only depend on the coordinates of the corners so, once again, it is not necessary to store the locations of the material points.

Computational framework
Algorithm 1: Computational framework 1 Set up problems: generate computational mesh, material points and particle domains, and initial material point volumes, V 0 p ; 2 for each load step do 3 for each material point p do 4 Find influence elements of p in the computational mesh; Initialise the out of balance force,   f R   ← 2 tol to ensure that the while loop is entered; 11 for each material point p do 12 assembling stiffness matrix k p (∆u) into the global stiffness matrix [K ]; 13 Update position of material points or particle domains; 19 Update volumes of material points, V p ; 20 Update the computational mesh with the moving mesh strategy; 21 end The previous section has detailed the continuum formulation, non-linear solution procedure and basis functions for the material point approaches analysed in this paper. This section provides a unified computational procedure for the implicit solution of quasi-static problems for the four MPMs. This framework guarantees that the simulation conditions, except for the basis functions related to different methods, are exactly the same. Therefore, any differences in the simulation results can only be attributed to the different methods used. The computational procedure is summarised in Algorithm 1. The following specific steps are needed for particular methods: • Line 1: For each material point, p, the sMPM only needs its coordinate, {x p }, and the CPDI1 also needs two vectors, ({s 0 }, {t 0 }), of the initial parallelogram particle domain. However, the CPDI2 approaches need corners of the particle domain and its connectivity.
• Line 4: The influence element is which encloses the material point for the sMPM. However, for CPDI approaches there are multiple influence elements, which enclose corners of the particle domain.

Numerical examples
A unified computational framework for the four MPMs was introduced in the previous section. In order to verify this computational framework, this section starts by solving three large deformation problems and compares the numerical predictions from the four MPMs with analytical or FEM solutions. The capabilities of these methods are then investigated via the solution of several problems involving large stretch, shear and torsion. We show that the CPDI2 approaches are less accurate than the standard MPM, under certain deformation fields, due to particle domain distortion. As the focus on this paper is a comparison of different MPMs, the problems have been carefully selected to avoid issues such as volumetric locking for elasto-plastic analysis. For detailed treatment of volumetric locking with different MPMs, see [27].
The presented simulations are conducted in two dimensions with a plane strain assumption, and all variables have compatible units. All of the analyses use a linear elastic, perfectly plastic constitutive model with the von Mises yield surface and associated plastic flow. The von Mises yield function has the following form

Confined column
The first example is compression of a column subject to body force. An analytical solution is available for this problem (see [26] for details), which can be used to verify and investigate the convergence behaviour of each method. Roller boundary conditions were applied to the base and sides of the column and an incremental body force ∆b = 10 4 was applied per load step. The material parameters were: a Young's modulus of E = 10 6 , a Poisson's ratio of ν = 0, and a yield strength of ϱ y = 2 × 10 5 . The initial height of the column was 20, and the width was set to 1. The mesh and material points at the start of the simulation and end of the first, second and third load steps are shown in Fig. 3 with 5 elements and 4 material points per initially populated element. The Cauchy stress components through the column at the end of the third load step are also shown in this figure. In the analytical solution, the stress component σ yy = σ x x and all shear stress components are zero. For the CPDI approaches, the Cauchy stress components agree with the analytical values, as shown in Fig. 3. However, the sMPM simulation suffers from cell crossing instabilities resulting in oscillations in the stress field in both the σ zz and σ x x components. In this problem, the particle domains in the CPDI1 method stay as rectangles during the simulation, so the results from the CPDI2q method are identical to those from the CPDI1 method.
To study the effect of mesh size and material point density on the convergence, relative stress error is defined as where superscript a indicates the analytical values, and p indicates the material point values from the sMPM and CPDI approaches, and ∥(·)∥ is the Euclidean L 2 norm of (·). The problem domain was discretised vertically by 5,  Fig. 4. The error in results from the sMPM is higher than that from the CPDI approaches because of the cell-crossing instability described above. The error is reduced by increasing the number of material points per element with the sMPM, while the error is reduced by increasing the number of elements with the CPDI approaches. We observe that increasing the number of elements cannot reduce the error in the sMPM. This is because of the increase in the number of instances that material points cross element edges. However, increasing the number of material points per element will reduce the volume associated with each material point, therefore the size of the discrete transfer of stiffness due to grid-crossing will be reduced. In the CPDI approaches, the alternative basis functions reduce the stiffness jump when material points cross element edges, therefore increasing material point density only has a minor impact on the maximum error.
The relative error is also plotted against either the number of material points per element or the number of elements in Fig. 5. With an increase in the number of material points per element, the error with the sMPM decreases, but increases slightly with CPDI approaches. With an increase in the number of elements, the error from the sMPM decreases in the initial refinement and then stays almost constant. By contrast, convergence (i.e. the rate of error reduction) with CPDI approaches has an approximately linear relationship with the number of elements.

Simple stretch
Simple stretching of a square domain is the next simulation to investigate the performance of these methods. A 2 × 2 square domain is used, comprised of a material with the following properties: Young's modulus E = 1000, Poisson's ratio ν = 0 and yield strength ϱ c = 400. In this case a finite element analysis (with four bi-linear quadrilateral elements integrated using 2 by 2 Gauss quadrature) is taken as the reference solution as the problem does not involve any mesh distortion. For the material point analyses, the same mesh as used in the finite element simulation is adopted and the physical material is discretised by four material points per element. Fig. 6 shows the geometry, mesh and boundary conditions for the CPDI2q method in both initial and deformed configurations. A roller boundary condition was applied on the bottom and left sides of the domain, with a horizontal displacement of 0.2 per load step applied on the right, and the simulation was run up for 20 steps. The moving mesh strategy [24,35] was utilised by horizontally stretching the mesh so that the edge of the physical domain was always aligned with the right end, where the displacement boundary condition is applied.
The convergence of the CPDI2q method is shown in Fig. 7(a) by plotting the residual, i.e. L 2 norm of (9), against the cumulative NR iteration steps. Each of the discrete lines is a single load step and each marker the iterations within each load step. The first six load steps are elastic and converge in two iterations whereas the remaining load   steps include elasto-plastic deformation and take four iterations to find equilibrium. Examination of the norm of the residual in an iterative step against that in the previous step in Fig. 7(b) indicates that the average convergence rate is 1.996 which is very close to the theoretical asymptotic quadratic convergence rate of the NR method. The sMPM and other CPDI approaches also have this correct quadratic convergence rate. This correct convergence rate indicates that the tangent stiffness matrix for large deformation elasto-plasticity is consistent within the out of balance force computation, verifying the correct implementation of the methods.
The reaction forces on the right end of the domain from the sMPM and FEM models are plotted in Fig. 8, where the markers indicate the load steps. The force increases nonlinearly in the first six steps due to the large deformation mechanics, i.e. the geometric nonlinear finite strain measure used in the computational framework. In the seventh step, the material yields and so this reaction force starts to decrease gradually. The results from the sMPM and CPDI approaches are the same as the FEM, as this simulation is not affected by mesh distortion. Also, as this is a constant stress problem with zero normal stress in the vertical direction the sMPM does not suffer from cell-crossing instabilities.

Manufactured solution of torsional deformation
This section presents an example using the Method of Manufactured Solutions (MMS) to verify the implemented numerical methods on a problem involving large torsional deformation using a problem similar to that presented in [36]. In the MMS, the first step is to postulate a deformation field and then, based on the governing equations, calculate analytically the body forces required to enforce this deformation field. In the numerical analysis, these body forces are imposed at the material points along with appropriate Dirchlet or Neumann conditions on the boundary of the domain. The numerical method can be verified by comparing the numerically computed deformation field to the postulated one.
In this example a circular annulus, as shown in Fig. 9, with inner (R i = 0.75) and outer (R o = 1.25) radii is subjected to a smooth rotational deformation field applied over 10 equal load steps. The material parameters are E = 10 6 , ν = 0 and a large yield stress value to ensure that the material remains in the elastic regime. The inner and outer radial surfaces of the annulus are fully fixed and the material points are subjected to the appropriate body forces. The mathematical expression of the deformation field and associated stresses are detailed in the Appendix. With this manufactured deformation field and the relevant stress, the body force for a quasi-static problem can be computed via (A.2).
The comparison of the numerical results from the four MPMs with the manufactured exact solution can be used to analyse the convergence of each method and quantify any difference in accuracy of the methods. The relative error in the displacement field is computed as where {u} This error is computed for several background meshes with different levels of refinement for each MPM. As seen in Fig. 9(b), the annular domain is discretised by 5 × 16 elements in radial and circumferential directions,  . Manufactured torsional deformation: the relative error in the displacement field for all four methods as indicated by markers: sMPM (circle), CPDI1 (diamond), CPDI2q (square) and CPDI2t (triangle), with the mesh including (5 × r f ) and (16 × r f ) elements in radial and circumferential directions, respectively. Note that the variance in the curves with r f = 8 cannot be seen due to the large scale in the vertical axis and please see the variance in Fig. 12(d).
respectively. Each element is initially populated with 4 × 4 material points. Notably, there is one layer extra of empty elements along both inner and outer radial surfaces for dealing with corners of some particle domains which may be outside of the mesh due to the curved boundary. The background mesh shown in Fig. 9(b) is the coarsest mesh, indicated by a refinement factor r f = 1. The finer meshes comprise (5 × r f )×(16 × r f ) elements, where r f is set to 1, 2, 4 and 8. The relative error with these meshes in all simulation steps shown in Fig. 10 shows that all of the methods have an error less than 3% with a mesh when r f ⩾ 2. A comparison of the convergence rate of the Fig. 11. Manufactured torsional deformation: the comparison of convergence with mesh refinement for four methods by plotting the arithmetic mean of ERR u across all of 10 load steps in Fig. 10 against the refinement level. mean error across the 10 load steps of each method with mesh refinement is shown in Fig. 11. All of the methods initially converge quadratically, however the figure shows that the CPDI1 method is starting to plateau whereas the other results are continuing to converge.
The comparison of accuracy of these methods for this problem across all load steps is shown in Fig. 12. With a coarse mesh, r f = 1, the CPDI1 produces the most accurate solution, followed by the CPDI2q, and the sMPM produces the solution with the maximum error. It is interesting to observe that the error with the CPDI2t method initially decreases and then increases with increasing deformation. The comparison is similar for r f = 2 and r f = 4. With a fine mesh, r f = 8, the error with the CPDI1 becomes the largest, while with the CPDI2q is the smallest, for α > 4 • . Examination of these four figures leads to the conclusion that there is no clear evidence that one method is better than others for this manufactured deformation field.   The error values are similar to that shown in Fig. 12, however the oscillations in the error for the sMPM are more clearly shown in Fig. 13. The cause of these oscillations can be attributed to cell-crossing as the other (domainbased) methods show a smooth variation in error with progressive rotation. It should be notated that the errors for the CPDI2 methods have a general trend of increasing error with rotation, which could be due to distortion of the particle domains, whereas the error associated with the CPDI1 method remains approximately constant.

Corner stretch
In order to test the performance of these methods for problems involving large shear, the next example is the stretching of a square domain along its diagonal. The material behaviour is purely elastic with Young's modulus E = 10 3 , Poisson's ratio ν = 0.4 and a very large value for yield strength ϱ c . Fig. 14 shows the initial configuration and loading of the problem domain which is discretised by a mesh of four quadrilateral finite elements with four material points per element, see Fig. 15. Roller boundary conditions are applied on both the left and bottom edges, while the top-right corner is stretched by 0.2 per load step in both directions. The moving mesh strategy was adopted such that the top and right edges moved with the top-right corner.
The deformed configurations for a displacement of u = 4 in the both directions are shown in Fig. 15. It is obvious that the movement of the top right material point is the smallest with the sMPM, followed by the CPDI1 method, while the CPDI2 approaches have the largest deformations. This is because the deformation is tracked by the continuous particle domains in the CPDI2 approaches, while the sMPM and CPDI1 use discrete material points to represent the physical domain. The loss of information during the mapping of displacement between material points and mesh nodes is less with the continuous particle domains than with the discrete material points.
Plots of the nodal reaction force at top-right corner against displacement with different densities of material points are shown in Fig. 16. The reaction force reaches a peak and then decreases in all simulations. With 4 material points  per element, the CPDI2t approach reaches the peak with the greatest gradient and predicts the largest reaction force. The other three methods lead to similar pre-peak results, but CPDI2q predicts a higher value than the others. This appears to be because the particle domains in the CPDI2q can track the deformation more accurately, as shown in Fig. 15(c). With 64 material points per element, the results from the four methods are similar, although the reaction force from CPDI2t is slightly larger than the others. It should also be noted that the deformed shape of the material point domains is significantly different for the CPDI2t method as compared to the other methods. This is because the derivatives of the basis functions ∇ S vp in the CPDI2t have some non-physical values, which is explained as follows.
Consider the standard reference element shown in Fig. 17 where the indices of nodes are denoted, v and material points are denoted, p. If we consider node 4 and material point 1, from (23), since the shape function S 4 at first two corners with coordinates (−1, −1) and (1, −1) are equal to zero, we can obtain the spatial derivatives of the basis function of this pair as where V p = 1 and S 4 (0, 0) = 1/4. We find that ∂ S 41 /∂ x ≡ 0, wherever the third corner is located. That means that the evaluation of the shape function S 4 (x, y) is independent of the first variable x, however this contradicts the definition of basis functions for both the underlying linear finite element grid and the method itself, as from (23) we can see that S 41 = 1 3 S 4 (x 1 , y 1 ). In addition, this means that the material point has zero contribution to the stiffness in the horizontal direction as the internal force in the x direction is not dependent on the normal stress in the horizontal direction, σ x x . This is a fundamental deficiency in the CPDI2t method and explains why its results are so different from the other MPMs.

Column collapse
The collapse of a column under self weight is the next example to be presented in this paper. The initial height of the column was 20 and the width 4, and was comprised of an elasto-plastic material with the following material parameters: Young's modulus E = 10 6 , Poisson ratio ν = 0 and yield strength ϱ c = 2 × 10 4 . The density was set to ρ = 80 and a gravity loading of ∆g = 20 per step was applied with 15 steps, giving a total gravitational load of g = 300. Due to symmetry only half of the column was analysed and a roller boundary condition imposed on the symmetry line. The base is also applied by a roller boundary condition.
The deformed profiles from using different methods are shown in Fig. 18 where the material points are shown as blue dots (for the CPDI methods the centre of the material point domain is plotted). The response of the sMPM, CPDI1 and CPDI2q methods are similar but the CPDI2t method shows very different material behaviour with the deformed profile reminiscent of the collapse of a purely frictional material. This response is at odds with the plasticity model used in the simulation and again potentially highlights a serious deficiency of the CPDI2t method.
The horizontal displacement of the lower-right corner of the column during the collapse is plotted in Fig. 18. In the first step, material response is elastic and all methods predict a very similar response. In the following plastic regime, the displacement-gravity response is more linear using CPDI2t than the other methods.

Azimuthal shear of an annulus
The final simulation to be presented in this paper is that of an annulus of elasto-plastic material subjected to internal twist with a fixed outer boundary. This example has been included to test the methods' performance under problems with large torsional deformation, which is related to our project for modelling of the soil response to the installation of a screwed-in pile foundation. Two different configurations are examined: circular internal and external boundaries, and an elliptical internal boundary and a circular external boundary. The problem modelled here is closely related to that studied in [37].

Circular annulus
The geometry, boundary conditions and finite element background mesh for this problem are shown in Fig. 19. The inner radius was R i = 5 and the outer radius R o = 10. In these simulations, incremental rotation was applied directly via a displacement boundary condition on the inner surface and the moving mesh strategy for rotation [24] was adopted to fix the background mesh to the inner circle.
To verify setup of this model, a small rotation with α = 10 • is firstly considered, and the numerical results are compared to the FE results. In this case the deformed finite element mesh suffers only minor distortion, it is therefore assumed that the finite element results provide a reasonable reference solution. The FEM [38] was used to analyse the problem with both fine (160 × 1152 elements) and coarse (10 × 72 elements) discretisations. 3 The agreement of the results between fine and coarse meshes (as shown in Fig. 20) suggests the coarse mesh is good enough to model this problem. Therefore, this coarse mesh was used as the computational mesh in the simulation with the sMPM and CPDI approaches. In all simulations, four material points per element were used. All simulations were run with a rotation increment of ∆α = 5 • in the clockwise direction up to α = 10 • in two load steps.
As the deformation field is axisymmetric, the magnitude of displacements along a radial line through the problem domain is compared in Fig. 20. The mesh nodes along a radial line with θ = 0 • were selected as sampling points, indicated by the markers in Fig. 19. In the FEM simulation, these nodal displacements were obtained directly. However, sampling points had to be added as infinitesimal-volume material points in the sMPM and CPDI1 methods, in order to obtain comparable total displacements. In contrast, the CPDI2 approaches do not need these material points. Instead the displacements at the corners of particle domains were used. The comparison in Fig. 20 shows that all of MPM variants can produce an accurate solution under small levels of rotational deformation.
All of the methods are now applied to simulate the annulus problem up to α = 80 • . The magnitude of reaction force along the inner circle against the rotation α is plotted in Fig. 21 for each method. The reaction forces are almost the same for all of the methods in the elastic region (α < 25 • ) but there is significant variation in the methods' elasto-plastic responses.
When the material yields, the reaction force is expected to decrease due to elastic unloading with continued plastic deformation of the yielded region close to the inner bore. This is observed in the results of the sMPM and CPDI approaches (at least initially). The reaction force in the FEM erroneously increases due to mesh distortion errors. The reaction force from the CPDI2t approach increases immediately after the post-yield drop, showing it to again face some issues in modelling. We have found that the CPDI2t approach performs differently to other CPDI methods because of the issue in the gradient of basis functions, as explained earlier (see Section 4.4). At around α = 35 • , a larger peak is observed in the sMPM than in the CPDI approaches. The drop following the peak appears to be due to stress relaxation when the inner material yields while the external material is still elastic, as shown in Fig. 22. The incremental deformation is shown in the deformed mesh for α = 30 • , 45 • and 60 • . It should be observed that at the α = 45 • step the first layer mesh nodes inside the physical domain rotate in the opposite direction, counter-clockwise, to the boundary condition which is applied clockwise on the inner circle, showing the stress relaxation.   Both the sMPM and CPDI1 approaches predict the same constant reaction force for α > 50 • , but the CPDI2q method has a slight increase on this response and the CPDI2t method is actually closer to the finite element response than the other MPMs. This erroneous increase in the CPDI2q is because of the distortion of particle domains. Recall that in the CPDI2 approaches the particle domains exactly follow the deformation as in a finite element mesh, but the distortion causes more serious errors in the FEM than the CPDI2 approaches.

Elliptical annulus
The performance of the material point methods was also investigated through the simulation of azimuthal shear with an elliptical inner boundary. In these simulations, the mesh at the ends of the long axis was locally refined because of the higher gradients in the pattern of deformation gradient and stresses expected at these positions. The deformed mesh with incremental displacement and material points in loadsteps α = 20 • , 30 • and 60 • is shown in Fig. 23. Due to the stress concentration at the ends of the long axis of the ellipse, the material yields earlier than for the circular annulus. As the ellipse rotates in the plastic material, two regions with fewer material points are created after the ends of long axis pass, as shown in Fig. 23 when α = 60 • . We also observed that the distortion of particle domains is significant in this simulation, as seen in Fig. 24. The magnitude of torque against α plotted in Fig. 25 shows that the CPDI2 approaches predict an erroneous increase in the torque due to particle domain distortion and that the degradation in the CPDI2t method is more severe than the CPDI2q. Both the CPDI1 method and the sMPM predict a more physically realistic response in the plastic region.

Summary of findings
This section has presented a number of numerical examples to verify implementation of the four different material point approaches and investigate their performance. In these examples, we have considered various deformation modes, including large stretch, shear and torsion, and different loading conditions, including body force and displacement boundary conditions.
Our implementations of these MPM and CPDI approaches are first verified in the column compression and simple stretch examples, by comparing the numerical solutions to analytical or FEM solutions. In addition, we have shown the convergence behaviour of these approaches, for the first time to our knowledge. It was found that the standard MPM will only converge with increasing numbers of material points per element (causing a reduction in cell-crossing instabilities), while the CPDI approaches will converge with increasing numbers of background elements for the same number of material points per element. The expected quadratic convergence rate of the NR iterative solver is also shown in the simple stretch example.
The MMS is also used to verify these methods and our implementation. The internal torsion of an elastic ring is simulated by applying body force, which is computed analytically from the manufactured deformation. However, the comparison of the methods for this problem shows that there is no clear evidence for that one method is more accurate than the others for this deformation field. Notably, the elastic material response will have a gradually varying deformation, but can not include the very large distortion for the elasto-plastic deformation, e.g. see Figs. 22 and 24. There is an indication that the CPDI1 method is starting to plateau with continued mesh refinement whereas the others continue to converge at a quadratic rate.
The other examples serve to illustrate some important features of the methods which will affect their ability to model certain deformation fields. The corner stretch example clearly demonstrates both the benefits and drawbacks of the CPDI2 approaches. They can clearly track the deformation of the physical domain better than the standard MPM and CPDI1 approach. In particular, the particle domains in the CPDI2 approaches can exactly represent the deformed domain. However, this characteristic is also a source of erroneous modelling prediction when using CPDI2 methods. This has been shown in the annulus example, where we observe the erroneous increase in the torque predicted using the CPDI2 approaches after the material has yielded. In contrast, the standard MPM and CPDI1 approaches predict physically reasonable responses -almost constant reaction force and torque after the material yields. This error is because of the distortion of the particle domains, as shown in Fig. 24 and might be surprising, i.e. that the later developed method (CPDI2) fails in comparison to earlier methods. Raising this awareness is very important, because this pattern of very large distortion in a region with material yielding is very common in geotechnical engineering, e.g. the installation of a screw pile, a shear vane test, etc.
We have also found that the results from the CPDI2 method with triangular particle domains are considerably different to those from the other methods in some problems. In the example of corner stretch, we have shown that there is an inherent error in calculation of the gradient of basis functions so that some components of the basis function gradients are zero, which leads to an unrealistic response compared to other approaches. It appears that there is work to do to remedy the CPDI2t approach.

Conclusion
This paper has presented a unified computational framework for the standard MPM and its latest particle domain based extensions, the CPDI approaches. The framework has been verified and then applied to problems involving large elasto-plastic deformation with different dominant deformation modes, namely: stretch, shear and torsion. CPDI2 approaches can increase the stability of material point analyses as they reduce cell-crossing problems and provide an accurate representation of the physical domain. However, they effectively include an unstructured mesh of particle domains that, like finite element methods, can suffer from erroneous results due to domain distortion, especially under torsional deformation. In addition to this the spatial gradients of the CPDI2t method's basis function degenerate under certain conditions leading to physically unrealistic results. The sMPM and CPDI1 are preferred to the CPDI2 approaches for simulating large elasto-plastic deformation problems involving torsional deformation modes.
where Div is the divergence with respect to the Cartesian basis ({E 1 }, {E 2 }, {E 3 }), [P] the 1st Piola-Kirchhoff or PK1 stress, ρ 0 is the initial density of the material, and {b} is the body force. For the problem of axisymmetric torsional deformation in Fig. 9, it is more convenient to compute the body force via the polar coordinate system than direct computation in the Cartesian system as shown in [36]. The polar coordinate of a point is denoted by (R, Θ) in the reference configuration and (r, θ ) in the deformed configuration.
The body force is given by Following equations (54)-(55) in [36] and neglecting the inertia term due to our focus on quasi-static analysis, the body force can be computed through where P F is the PK1 stress associated with the deformation F, and ϵ is the shear strain which is defined as and F is an angle-independent 'baseline' deformation representing simple shear without superimposed rotation. t is a pseudo-time for indicating particular loading in this quasi-static analysis. h is function of R, and g is function of t. α is the rotation angle which is dependent on the radial coordinate of the material point, R, and the imposed degree of deformation through α(R, t) = h(R)g(t), (A. 8) where h(R) controls the radial deformation field and g(t) the deformation magnitude In the above equation, R i and R o are inner and outer radii and α 0 is the maximum imposed rotation. The deformation gradient expressed in the Cartesian frame is where Θ is the circumferential coordinate of a point in the reference configuration.

A.2. Stress and its derivatives
The formulation used in this paper adopts a linear relationship between Kirchhoff stress and logarithmic strain which is a function of shear strain, ϵ. Finally, we use a finite difference technique to approximate the derivative of [P F ] which respect to the imposed shear strain, that is, where δϵ is an infinitesimal increment in the shear strain and was taken to be 1 × 10 −6 . These derivatives are used in (A.4)-(A.5) to determine the body force.