On Lagrangian mechanics and the implicit material point method for large deformation elasto-plasticity

The material point method is ideally suited to modelling problems involving large deformations where conventional mesh-based methods would struggle. However, total and updated Lagrangian approaches are unsuitable and non-ideal, respectively, in terms of formulating equilibrium for the method. This is due to the basis functions, and particularly the derivatives of the basis functions, of material point methods normally being defined on an unformed, and sometimes regular, background mesh. It is possible to map the basis function spatial derivatives using the deformation at a material point but this introduces additional algorithm complexity and computational expense. This paper presents a new Lagrangian statement of equilibrium which is ideal for material point methods as it satisfies equilibrium on the undeformed background mesh at the start of a load step. The formulation is implemented using a quasi-static implicit algorithm which includes the derivation of the consistent tangent to achieve optimum convergence of the global equilibrium iterations. The method is applied to a number of large deformation elasto-plastic problems, with a specific focus of the convergence of the method towards analytical solutions with the standard, generalised interpolation and CPDI2 material point methods. For the generalised interpolation method, different domain updating methods are investigated and it is shown that all of the current methods are degenerative under certain simple deformation fields. A new domain updating approach is proposed that overcomes these issues. The proposed material point method framework can be applied to all existing material point methods and adopted for implicit and explicit analysis, however its advantages are mainly associated with the former.


Introduction
Numerical approaches based on Lagrangian mechanics, particularly using the finite element method, dominate the analysis of solid mechanics problems. However there are issues with mesh-based Lagrangian approaches for problems involving very large deformation in that the discretisation used to analyse the physical problem becomes distorted leading to inaccurate results and, in extreme cases, eventual breakdown of the numerics due to element inversion. The material point method (MPM) is an alternative to pure Lagrangian approaches and is well suited to problems involving very large deformations. The method was developed in the 1990s by Sulsky and co-workers [1] as a solid mechanics extension to the FLuid Implicit Particle (FLIP) method [2] which itself was developed from the Particle-In-Cell (PIC) method [3]. In the material point method, equilibrium computations take place on a background grid but the calculations are based on information (mass, volume, stress, state variables, etc.) held at material points that are convected through the background grid as the material deforms. This allows computations to take place on an undistorted background mesh (structured or unstructured) whilst modelling problems involving very large deformations. One way to summarise the material point method is: a finite element method where the integration points (material points) are allowed to move independently of the mesh.
One of the most commonly cited drawbacks of the material point method is its inherent cell-crossing instability caused by the sudden transfer of mass, stiffness, internal force, etc. as a material point crosses between background elements. There have been several approaches proposed to mitigate this issue including the Generalised Interpolation Material Point (GIMP) method [4], the Convected Particle Domain Interpolation (CPDI) method [5], the second order Convected Particle Domain Interpolation (CPDI2) method [6,7] and Dual Domain Material Point (DDMP) methods [8][9][10][11]. Conceptually these methods are similar in that the discrete material point used in the standard material point method is replaced by a deformable material point domain which provides a gradual transition of information as a material point moves between background grid elements. An alternative is to adopt an additional level of smoothing, such as in the moving least squares modified material point method [12,13] or the weighted least squares approach of Wallstedt and Guilkey [14]. Another alternative is to use different basis functions which offer increased smoothness between elements, such as B-splines or cubic splines [15][16][17][18]. In this paper we examine the standard, generalised interpolation and CPDI2 material point methods, however the findings from this research can be applied to other material point methods with only modification of the basis functions. Three different material point methods are considered as, in the authors' experience, different material point variants are well suited to modelling different types of problems and some of the latest material point methods can become unstable under certain deformation fields (such as large distortion, shear and/or rotation), these points are explored in detail in Wang et al. [19].
The majority of material point papers implement the method using an explicit approach [20], however there are benefits of using an implicit approach [21][22][23][24][25] namely: (i) improved stability, (ii) error control and (iii) increased load step size. However, implicit approaches are more complex than explicit approaches in that they require the discrete equilibrium equations to be linearised with respect to the primary grid unknowns (typically incremental displacements). As the focus of this paper is on equilibrium formulations within the material point method, we adopt an implicit approach for quasi-static analysis as this allows precise error control when satisfying the governing equations.
Although the material point method is focused on solving large deformation problems, the precise application of finite deformation Lagrangian mechanics is a neglected area in the development of material point methods to date. In the authors' opinion this is due to a number of reasons: (i) the vast majority of current material point method codes are based on explicit time integration 1 ; (ii) it is assumed that standard finite element technology can be adopted without modification; and (iii) the majority of recently published material point method papers are focused on the application of the method rather than its fundamental development. 2 Lagrangian mechanics for large deformation analysis is broadly split into two approaches; total and updated formulations which use material and spatial descriptions of motion, respectively. In particular they differ in the frame in which equilibrium is satisfied. In total Lagrangian approaches the equilibrium equations are formulated based on the reference, or original, frame whereas in updated Lagrangian formulations equilibrium is satisfied at the current, or deformed, state. Most material point approaches use an updated Lagrangian formulation as it is a more natural fit with the method. This is because, unlike conventional finite element methods, in the material point method information between the original and current background meshes (where equilibrium is actually satisfied) is lost by the fact that at the end of each load/time step the mesh is reset or replaced. Although possible, this makes mapping information, specifically spatial derivatives of the basis functions, back to the reference frame cumbersome for a total Lagrangian approach (this point is expanded upon later in the paper).
This paper first reviews the applicability of total and updated Lagrangian mechanics formulations for use with the material point method. The paper then proposes a new approach which formulates the statement of equilibrium in the previously converged frame. This previously converged Lagrangian formulation is well suited to material point methods as, unlike total/updated Lagrangian approaches, it allows all calculations to take place on an undeformed background mesh from the start of the load/time step. The approach allows the basis functions (and their spatial derivatives) to be calculated once per load step and used until equilibrium is achieved, reducing one of the significant costs in domain-based material point methods (GIMP, CPDI1, CPDI2, etc.). In addition to this, for the generalised interpolation material point method the basis functions are only defined in terms of regular background grids and the use of an updated Lagrangian approach requires the derivatives of these basis functions to be mapped into the current frame [21,26]. The paper also explores different domain updating procedures for the generalised interpolation material point method and shows that the existing approaches are inappropriate for some basic deformation fields. A new approach is proposed which overcomes the limitations of existing methods.
The formulations presented in this paper can be applied to both implicit and explicit approaches but it is with the former where the advantages of the formulation are apparent. This is because explicit approaches advance a simulation through forward prediction based on the state of the analysis at the start of the load (or time) step. In this case the updated Lagrangian and previously converged Lagrangian formulations are identical as the load step has yet to experience any deformation.

Large deformation mechanics
This paper is restricted to quasi-static finite deformation mechanics with isotropic elasto-plastic material behaviour with a multiplicative decomposition of the deformation gradient into elastic and plastic components combined with a linear relationship between elastic logarithmic strains and Kirchhoff stresses. For the sake of brevity details of the large deformation framework are not given here, instead the interested reader is referred to [19,21,26] for more information.

Updated Lagrangian formulation
This section details the quasi-static implicit finite deformation elasto-plastic material point method based on an updated Lagrangian formulation 3 defined by the following weak statement of equilibrium ∫ ϕ t is the motion of the material body with domain, Ω , which is subjected to tractions, t i , on the boundary of the domain (with surface, s), ∂Ω , and body forces, b i , acting over the volume, v of the domain, which lead to 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 .

Discrete material point formulation
The Galerkin form of the weak statement of equilibrium over each background grid element, E, can be obtained from (1) as ∫ where [∇ x S vp ] is the tensorial form of the strain-displacement matrix containing derivatives of the shape functions with respect to the updated coordinates and {σ } is the nine-component vector form of σ i j . The first term in (2) is the internal force within an element and the combination of the second (body forces) and third (tractions) terms is the external force vector. Eq. (2) is non-linear in terms of the unknown nodal displacements and can be efficiently solved using the standard implicit Newton-Raphson (N-R) procedure (see Section 6.3). This requires (2) to be linearised with respect to the nodal displacements to arrive at the stiffness of each element in the background mesh where [a] is the matrix form of a i jkl = (∂σ i j /∂ F km )F lm , which is the spatial consistent tangent modulus for a point within the element which is detailed in Appendix A.1.
In material point methods the physical domain is discretised by a number of material points. These points are used to numerically approximate the stiffness (3) of the elements in the background mesh, essentially replacing the conventional Gauss points (or other integration method). The key difference between material point and finite element methods is that these integration points move relative to the background mesh rather than being directly coupled to the positions of the background grid nodes. The stiffness contribution of a single material point to the background mesh is where V p is the volume associated with the material point in the spatial (updated) frame where V n p and V 0 p are the volume associated with the material point in the previously converged state and the initial configuration, respectively. The internal force contribution of a single material point to the background mesh is Following the work of Charlton et al. [21], and as adopted by Coombs et al. [26], the increment in the deformation gradient is obtained from where ∆u i is the displacement increment within the current load step,X j are the coordinates at the start of the load step and n is the number of nodes that influence the material point. This allows the increment in the deformation gradient to be obtained from derivatives of the basis functions based on the coordinates of the nodes at the start of the load step. This is essential as in material point methods there is no concept of the current (deformed) nodal coordinates as information is lost between incremental steps. Also, the basis functions for material point methods are typically only defined on a regular background grid. The spatial derivatives of the basis functions can subsequently be calculated using the method proposed in Charlton et al. [21], that is It is essential that the spatial derivatives are used in the strain-displacement matrix, [∇ x S vp ], to both ensure the correct order of convergence in the N-R process and convergence towards the correct solution based on the internal force contribution, (6), of each material point [21,26]. For the standard material point method an alternative to mapping the derivatives of the basis functions into the spatial frame using (8) is to calculate the spatial derivatives of the basis functions using the standard finite element Jacobian mapping based on the deformed nodal positions. However, this is only possible for certain material point methods as some domain-based material point formulations, such as the generalised interpolation material point method, rely on a regular background grid.
It is not necessary to map the basis functions between the start of the load step and the current configuration because it is assumed that, during a load step, the displacement of a material point, (∆u p ) i , is linked to the nodal displacements through the basis functions, S vp , that is where n is the number of nodes that influence the material point and (∆u v ) i are the displacement of the nodes over the current load step. Therefore the basis functions evaluated at the start or the end of the load step are identical, as shown in Fig. 1. The figure shows a single material point on a background mesh at the start and end (deformed configuration) of the load step. It also shows the local coordinates, ξ i , of the material point's parent element and gives the basis function for a single node (dark-grey shaded). The local position of the material point is the same throughout the load step and therefore the basis functions are also constant, whereas the spatial derivatives of the basis functions are dependent on the deformation over the load step.

Total Lagrangian formulation
This section details the quasi-static implicit finite deformation elasto-plastic material point method based on a total Lagrangian formulation defined by the following weak statement of equilibrium ∫ where the reference domain, Ω , is subjected to tractions, t i , on the boundary of the domain (with surface, s), ∂Ω , and body forces, b i , acting over the volume, v of the domain, which lead to a first Piola Kirchhoff stress field, P i j = J σ ik F −1 jk , throughout the body. The weak form is derived in the reference (or material) frame assuming a field of admissible virtual displacements, η i , and the derivatives of these virtual displacements in the first term in (10) are taken with respect to the original material coordinates, X i .

Discrete material point formulation
The Galerkin form of the weak statement of equilibrium over each background grid element, E, can be obtained from (10) as where [∇ X S vp ] is the tensorial form of the strain-displacement matrix and contains the derivatives of the basis functions, [S vp ], with respect to the original coordinates. The first term in (11) is the internal force within an element and the combination of the second (body forces) and third (tractions) terms is the external force vector. In material point methods, at the end of each load (or time) step the background mesh is reset or replaced. This causes issues when trying to adopt a total Lagrangian formulation as information is lost regarding the deformation of the background mesh between load steps. For example, the strain-displacement matrix requires the derivatives of the basis functions with respect to the original coordinates, but there is no guarantee that a material point is in the same element as at the start of the analysis. It is possible to map the derivatives of the basis functions back to the reference frame using but these derivatives would be associated with the material point's current element. There are also issues associated with determining the external force vector in the reference frame due to the material points moving between elements. Some of these issues are highlighted by Fig. 2 which shows the motion and deformation of a material point from the reference frame to the current configuration. The material point is initially located in the grey-shaded element and over a number of load steps it has moved toX i and in the current load step has deformed to x i with an associated increment in the deformation gradient. The deformed mesh in the current load step is shown by the grey dashed line and the original mesh by the grey solid lines. Consider the calculation of the nodal body forces associated with the material point as it convects through the background mesh. After the initial load step the material point's local position within an element, and its associated basis functions, will change. Therefore, if we attempted to satisfy equilibrium in the reference frame using (11) we would have to map the basis functions back to the original coordinates. In the authors' opinion this is not practical as it would require mesh deformation throughout the analysis to be stored, destroying one of the key advantages of the material point method. Therefore we feel that it is inappropriate to combine a total Lagrangian formulation with the material point method.

Previously converged formulation
This section provides a new approach for satisfying equilibrium in the material point method. 4 The updated Lagrangian weak statement of equilibrium (1) can be recast within a previously converged formulation where equilibrium is satisfied at the starting point of a load step (or the previously converged state). The statement of equilibrium becomes ∫ ϕ tn (Ω) where ϕ t n is the motion of the material body evaluated at the previously converged (or for the first load step, the initial) state andP i j is the Cauchy stress pulled back to the previously converged state,X i , that is where σ i j is the Cauchy stress, ∆F i j is the deformation gradient increment over the load step and ∆J = det(∆F i j ) is the volume ratio between the current and previously converged state. (14) can be expressed in terms of the Kirchhoff stress, τ i j as where J = det(F i j ). Using the fact that the current deformation gradient can be obtained using where (F n ) i j is the deformation gradient from the previously converged state, and (15) becomes This particular form is useful when linearisingP i j with respect to the deformation gradient increment over a load step as J n is constant over the load step. The next section presents the discrete material point approach for this previously converged formulation.

Discrete material point formulation
The Galerkin form of the weak statement of equilibrium for the approach described above over each background grid element, E, can be obtained from (13) as where [∇X S vp ] is the tensorial form of the strain-displacement matrix containing derivatives of the basis functions with respect to the nodal coordinates at the start of the load step. As with (2), (19) is non-linear in terms of the unknown nodal displacements and can be solved in the same way, the only difference being where equilibrium is satisfied. The N-R solution procedure is detailed in Section 6.3. The element stiffness matrix can be obtained by linearising the discrete statement of equilibrium with respect to the unknown nodal displacements to give whereÃ i jkl = ∂P i j /∂∆F kl is the previously converged consistent tangent modulus and its derivation is detailed in Appendix A.2. The stiffness contribution of a single material point to the background mesh is where V n p is the volume associated with the material point at the end of the previously converged (or the start of the current) load step The internal force contribution of a single material point to the background mesh is The advantage of the proposed formulation is that it does not require the derivatives of the basis functions to be mapped into the spatial (8) or material (12) frames as in the updated and total Lagrangian formulations, respectively. Therefore all equilibrium calculations truly take place on the background mesh at the start of the load step before the material points are convected through the mesh once equilibrium has been obtained.  Table 1 Comparison of key quantities required for total, updated and previously converged Lagrangian approaches.

Lagrangian formulation Total Updated Previous
Coordinates The total, updated and previously converged formulations are compared and contrasted in Fig. 3 and Table 1. Fig. 3 shows a deforming domain, the grey-shaded body, in the reference, previous and deformed configurations along with the mesh associated with these frames. The reference, X i , and previous,X i , configurations are linked by the deformation gradient at the previously converged state, (F n ) i j , and the stress measures that are used to satisfy equilibrium in these frames are P i j andP i j , respectively. The previous and current, x i configurations are linked by the deformation gradient increment over the load step, ∆F i j , and the Cauchy stress, σ i j , is used to satisfy the equilibrium equations in this deformed frame. The dashed shapes shows the domain from the preceding frame, for example in the previous configuration the dashed line shows the domain in the reference configuration and for the deformed configuration the dashed line shows the previous configuration. In the reference and previous configurations calculations take place on undeformed background meshes whereas in the updated Lagrangian approach equilibrium is satisfied on the current, deformed mesh. Table 1 highlights key quantities associated with the different formulations.
The proposed previously converged formulation shares some similarities with incremental non-linear elasticity where small deformations are superimposed on a given equilibrium state, see for example Bigoni [27]. However, in that work it is assumed that the deformation between the current and previous state is small, leading to a linear set of equations, whereas no such assumption is made here and we allow for finite (and non-linear) deformations within each step. There are also similarities with co-rotational formulations 5 that decompose the total deformation into stretch and rotation and map the equilibrium coordinates according to the material rotation. However, the formulation proposed in this paper uses the entire deformation (stretch and rotation) when defining the equilibrium state based on the deformation at the end of the previous load step.

Implementation aspects
This section presents key information regarding the implementation of the finite deformation formulations detailed in the previous sections. This includes the basis functions for the implemented material point methods, boundary conditions and material point position/domain updating.

Basis functions
The basis functions for the standard material point method are the same as conventional finite element analysis and the basis functions for the other methods are available from: [4,21] for the generalised interpolation method and [6] for CPDI2 with quadrilateral domains. 6

Boundary conditions
In this paper we ignore the external tractions as their general implementation within material point methods has yet to be fully clarified due to the uncertain definition of the boundary of the physical domain. It is noted that in domain-based material point methods (GIMP, CPDI1, CPDI2) it is possible to specify surface tractions on the outer boundaries of the material point domains and Bing et al. [29] and Remmerswaal [30] have gone some way to addressing this issue for the standard material point method.
In this work, Dirichlet boundary conditions are imposed directly on the background mesh in the same way as the standard finite element method. This places some restrictions on the form of the boundary conditions that can be applied to the material point method, which have been overcome by the recent work of Cortis et al. [31], Bing [32] and Bing et al. [29]. However, the focus on this paper is the use of different continuum mechanics formulations with the material point method and all of the problems analysed in Section 7 can be modelled by imposing the Dirichlet boundary conditions on the background grid.
The body force in (2) or (19) is approximated using where {b} is the body force associated with the material point. Point forces, if required, are held and convected with the material points. They are mapped to the background grid using where { f p } is the point force associated with the material point and { f P } are the equivalent nodal values.

Non-linear solution procedure
The nodal displacements within a load step, {∆d}, can be obtained by iteratively updating the nodal displacements until (2) or (19), for the updated and previous Lagrangian formulations respectively, are satisfied within a given tolerance using where k is the current iteration within the load step, [K ] is the global stiffness matrix and {δd k } is the iterative increment in the displacements from that iteration. The global stiffness matrix is obtained through assembling the individual material point contributions, for example for the previously converged formulation where A is the standard assembly operator acting over all of the material points in the problem and [k p ] is given by (21). { f R k−1 } is the residual out of balance force vector associated with the previous displacement value; the difference between the internal forces due to the stresses within the material and the applied boundary conditions, as given by (19) for the previously converged formulation. The global residual vector is assembled through where the first term is the internal force vector and the external force vector, { f ext }, which is constant over the load step, is given by The current displacement in a load step can be obtained by summing the iterative increments within the load step, that is where it is assumed that {δd 0 } = {0}.

Position and domain updating
At the end of each load step the material point positions, volumes and (if appropriate) domain half-lengths, l p , should be updated. The updated positions of the material points at the end of the load step are given by where (∆u p ) i is the displacement of the material point over the load step.
In the case of the standard material point method the volume at each material point is simply updated from the original volume using the determinant of the deformation gradient, (5). In other methods that explicitly define a finite sized domain for each material point the process of domain updating is more complex. For generalised interpolation material point methods, one approach is to apply uniform expansion/contraction of domains according to the volume change at the material point using where n D is the physical dimensionality of the analysed problem. A commonly adopted approach is to update the domain sizes according to the normal components of the deformation gradient through where F ii denotes one of the leading diagonal components of F i j . An alternative domain updating method, proposed by Charlton et al. [21], is to map the domain sizes according to the normal components of the material stretch tensor, that is The suitability of these domain updating methods for different deformation modes is summarised in Table 2, which shows that both not updating the material point domains and updating them according to (32) are not suitable choices when analysing problems involving simple stretch as they fail to correctly account for changes in the physical size of the global domain in the stretched/compressed direction. The only exception to this is when modelling problems involving hydrostatic compression/expansion where (32) is a suitable choice due to the imposed uniform volume change. (33) and (34) are suitable for these types of problems, however they have issues under pure rotation and simple shear, respectively. These issues are demonstrated in Fig. 4. Fig. 4 (a) shows the case of  pure rotation where the domain is updated according to (33), where the light and dark grey shaded regions are the original and updated domains, respectively. The mid lines of the domain are identified by the solid black lines with the black-shaded circles and the deformed location of these lines are shown by the black dashed lines with the white-shaded circles. Under continued rotational deformation the material point's domain shrinks but the physical volume of the particle should remain constant due to the isochoric nature of the deformation field. Fig. 4 (b) shows the case of pure shear where the domain is updated according to (34). In this case the domain expands in the y-direction whilst shrinking in the x-direction and the overall volume of the domain increases however, as in the case of pure rotation, the physical volume of the particle should remain constant due to the isochoric nature of the deformation field.
Here we propose a new domain updating method which is inspired by CPDI2 approaches in that the individual corners of a particle's domain are updated according to the nodal deformation but, in order to maintain a rectangular domain, the midpoints of the updated limits of the domain are used to define the new l p i . The proposed procedure is shown in Fig. 5 (a) for a single generalised interpolation material point domain. The original domain is shown by the grey-shaded region and the original and the deformed corners of the domain identified by the white and grey-filled circles, respectively. The corners of the domain are updated according to the basis functions associated with the background grid, that is whereX c i is the corner's position at the start of the load step. Once the updated positions of the domain corners have been determined, the midpoints of the edges can then be used to determine the new domain lengths through where x m i are the updated locations of the domain edge midpoints. An alternative is to directly map the midpoints using the nodal deformations, as shown in Fig. 5 (b). The final part of the updating procedure is to ensure that the domain updating is consistent with the change in volume at the material point. This is achieved using which uniformly scales the size of the domain about the material point's location and accounts for any spurious changes in volume caused by pure rotational deformation. Note that it is possible to apply a similar correction to the (33) domain updating method to ensure that the correct volume is maintained at the material point, such as However, there is still an issue with this approach when components of F ii are equal to zero, such as the case of 90 • rotation. An incremental form of (38), such as where l p n i is the domain length from the previously converged state, partially overcomes this issue, provided that excessive rotations are not experienced within a single load step. (32), (33) and (34) can also be cast in an incremental setting, if appropriate or convenient. Although the focus of the numerical examples presented in this paper is on two-dimensional analysis, it is straightforward to extend all of the above approaches to three dimensions and the above equations allow for this.

Computational procedure
This paper has described two approaches (updated Lagrangian and previous Lagrangian) to solve quasi-static large deformation elasto-plastic problems using implicit formulations of the standard, generalised interpolation and CPDI2 material point methods. Details of the computational procedure can be found in Coombs et al. [26].

Numerical examples
All of the results in this section use the previously converged formulation presented in Section 5. Adopting the updated Lagrangian formulation, as expected, gives identical results but is computationally less efficient due to the mapping of the basis derivatives between the previously converged state, where they are defined, and the current deformed state. Based on studies conducted by the authors, the computational saving of adopting the previously converged formulation is approximately 5% for two-dimensional analysis.
In all cases the global tolerance on the normalised out of balance force in the N-R process was set at 1 × 10 −9 , where the normalised out of balance force was defined as ∥(·)∥ denotes the L2 norm of (·), { f ext } includes the external forces applied to the problem (body forces, tractions, point loads) and { f react } are the reaction forces due to the prescribed displacement boundary conditions.

Compaction under self weight
The first example is an elastic column compressed under its own weight. Initially the column had a height of l 0 = 50 m and a width linked to the size of the background grid used to analyse the problem such that there was always one element in the horizontal direction. A plane strain assumption was set in the third direction. The base of the column was restrained vertically and both sides of the column were restrained in the horizontal direction. The column had a Young's modulus of 1MPa and Poisson's ratio of zero and a density of ϱ 0 = 80kg/m 3 . A body force of −800N/m 3 was applied in the vertical (Z , z) direction over 20 equal loadsteps.
The analytical solution for the vertical Cauchy stress is where g is the gravitational acceleration (taken to be 10 m/s 2 ), l 0 is the initial height of the column and Z is the initial vertical position within the column. If Poisson's ratio is zero the stresses in the horizontal directions are equal to zero when the behaviour is elastic. The convergence behaviour of the standard material point method with different numbers of material points per initially active background grid cell is shown in Fig. 6, where the dimensionless error is defined as In the above n p is the total number of material points, (σ p ) zz is the Cauchy stress in the vertical direction at a material point, Z p is the material point's original position, V 0 p is the original volume associated with the material point and σ a zz is the analytical Cauchy stress solution in the vertical direction, given by (41). The denominator of (42) is the product of the vertical Cauchy stress at the base of the column multiplied by the column's original volume, V 0 .  Fig. 6 presents the convergence data in two ways: (a) in terms of the normalised stress error against 1/ h z for different numbers of material points per background grid and (b) against 1/l 0 p with the data points being linked by the size of the background grid. 7 While it has been reported in the literature that the standard material point method does not converge for this problem due to cell crossing instabilities [21,26]. Fig. 6 shows this is not the case, but in order to achieve continued convergence with grid refinement a large number of material points is required. For this problem the required minimum number of material points per element in the vertical direction appears to be linked to For a given background grid size, it can be seen if sufficient material points have been used by the lines that link the discrete data points in Fig. 6 (b). These lines link together common background grid sizes and sufficient material points have been used if the error remains constant with increasing 1/l 0 p . It is clear that for h z = 0.3906 m and h z = 0.1953 m insufficient material points have been used for the background grid size and for these analyses it is predicted that 1,024 and 4,096 material points per element are required to minimise the error. Fig. 7 shows the convergence of the CPDI2 and generalised interpolation material point methods with two material points in the vertical direction per initially populated background grid cell. For the generalised interpolation method different domain updating methods are included. When the domains are not updated or when the domains are updated according to the determinant of the deformation gradient (32), the method fails to converge beyond the third and sixth analyses, respectively. However, when the domains are updated according to (33), (34) or (37) the method continues to converge with background grid refinement at the same rate as the CPDI2 method.
The method fails to converge when the domains are not updated due to the physical change in volume due to compaction under self weight not being accounted for at the material points, resulting in an over-stiff response. When the domains are updated according to (32) they contract in the horizontal direction due to the vertical compression. This causes a spurious reduction in the overlap between the background grid elements and the material point's domain and insufficient compression of the domain in the vertical direction. The result of this is that the method fails to converge towards the correct analytical solution beyond the sixth refinement step. 7 As the standard material point method is used in this example the length of the material point domains are obtained by assuming

Simple shear
The second example is a plane strain unit square subjected to simple shear with elastic and elasto-plastic material behaviour. The updated coordinates and deformation gradient for this problem are where det([F]) = 1, l 0 the original height of the domain and ∆l the top face shearing distance (as shown in Fig. 9). For the elasto-plastic analysis yielding of the material was governed by a von Mises yield function of the form where ρ = √ 2J 2 , J 2 = s i j s ji /2, s i j = τ i j − τ kk δ i j /3 and ρ y is the deviatoric yield stress of the material. The unit square had a Young's modulus of 10 kPa and Poisson's ratio of 0.3 and a density of ϱ 0 = 0. The top surface of the 1m square was sheared by ∆l = 6 m over 30 equal displacement controlled loadsteps. The shearing was imposed on all of the background cells that contained the outer boundary of material points, as shown by the black-filled circles in Fig. 9.  The normalised stress versus deformation response for an analysis using the standard material point method with h = 0.125 m and 3 2 material points per initially active background grid cell is shown in Fig. 10 with (a) elastic and (b) elasto-plastic material behaviour with ρ y = 5 kPa. This is a constant stress problem and all of the material points have the same stress response (to machine precision). The analytical solution for elastic behaviour is shown by the solid grey line in Fig. 10 (a) and the material point solution agrees with this solution to a normalised error of the order of 1 × 10 −5 . As this is a constant stress problem, the standard material point method does not suffer from the cell crossing problems in the previous example.
The standard, generalised interpolation and CPDI2 material point methods predict the correct response for the constant stress problem. The only exception was when the domains were updated according to the stretch tensor,  U i j , through (34), as explained in Section 6.4. This domain updating method introduces spurious volume changes and domain expansion/contraction under pure shear which should be an isochoric deformation field. The expansion of the domains in the vertical direction causes the areas associated with the top and bottom layer of material points to go beyond the boundary of the physical problem, losing partition of unity of the basis functions and causing a spurious variation in the stress field.

Elastic cantilever beam
The third example is a plane strain analysis of an elastic beam subjected to a vertical end load of 100kN applied over 50 equal loadsteps. The beam had an initial length l 0 = 10 m and depth d 0 = 1 m and the material had a Young's modulus of 12MPa and a Poisson's ratio of 0.2. The geometry and boundary conditions with h = 0.5 m and 2 2 material points per initially active background elements is shown in Fig. 11. The force was split across the two material points closest to the neutral axis at the end of the beam (the black-filled circles in Fig. 11). The neutral axis of the beam was pinned at its left hand edge and roller boundary conditions were imposed on the other background grid nodes at the root of the beam, as shown in Fig. 11.
The  The normalised horizontal and vertical displacement versus normalised force response for different background grid sizes and numbers of material points are shown in Fig. 13. The displacements have been normalised by the original length of the beam, l 0 , and the force by the maximum applied force of 100kN. Fig. 13 (a) shows the force-displacement response for h = 0.5 m and 3 2 , 6 2 , 9 2 and 12 2 material points per initially active background grid cells. 8 The analytical solution of [33] is shown by the discrete points with the horizontal and vertical solutions being shown by the black-and white-filled circles, respectively. Increasing the number of material points reduces the oscillations in the force-displacement response and the global response curves converge with increasing numbers of material points. The oscillations are caused by the sudden transfer of force and stiffness between background grid Fig. 13. Elastic beam: normalised force versus horizontal and vertical displacement response for the standard material point method. The analytical horizontal (black-filled circles) and vertical (white-filled circles) solution of [33] is also shown. cells as material points cross between background grid cells. Increasing the number of material points increases the instances of cell crossing but each of the material points has a lower volume resulting in a more gradual transfer of force and stiffness. Similar behaviour is seen in Fig. 13 (b) and (c) with background grid sizes of h = 0.25 m and h = 0.125 m, respectively. Fig. 13 (d) shows the response of mesh sizes h = 0.5 m, h = 0.25 m and h = 0.125 m with 6 2 material points per initially active background grid element. There is little difference between the forcedisplacement responses and the number of material points appears to be more important than the background grid size for this analysis.   (32). In this case, due to the large rotation of the beam, the material point domains artificially shrink and the force-displacement response deviates from the correct solution after approximately 15% of the load has been applied and fails to converge beyond load step 19 (highlighted by the grey-filled circles in Fig. 15 (b)). Similar behaviour of the (32) updating method was shown in Charlton et al. [21]. Fig. 16 shows the deformed generalised interpolation material point positions coloured according to σ yy for 6 2 material points per initial grid cell and h = 0.5, 0.25 and 0.125m background grid sizes. In these analyses the material point domains were updated according to (37) (volume-corrected corner update with mid-side determination). The detail around the loaded end of the cantilever beam is shown for h = 0.125 m which clearly shows the stress concentration caused by loading the two mid-depth material points (enclosed by the dashed-circular line, A). This stress concentration causes the over estimation of the vertical displacement shown in Figs. 13 and 15. Fig. 17 shows the normalised force versus displacement response of the CPDI2 method with 3 2 material points per initially populated background grid cell. The response is very similar to that of the generalised interpolation material point method provided that an appropriate domain update procedure is adopted. Table 3 shows the equilibrium convergence of the generalised interpolation material point method with h = 0.125 and 6 2 material points per element. The out of balance force residual (40) is shown for loadsteps 10, 20, 40, 49 and 50 with the other loadsteps having a similar response. In all cases asymptotic quadratic convergence is achieved as the algorithm approaches equilibrium satisfaction. This can be seen by the negative exponent on the residual  doubling with between iterations. This verifies the derivation of the previously converged consistent tangent (A.3) and its implementation within the generalised interpolation material point methods for large strain elasticity. Table 4 presents the normalised run times and total number of N-R iterations for the standard, generalised interpolation and CPDI2 material point methods where the run times have been normalised with respect to the finite element run time with h = 0.5 m. All of the material point analyses used 6 2 material points per initially populated element and the material point method implementation was based on [34]. The reference finite element code was taken from Coombs et al. [35], with four-noded bi-linear quadrilateral elements integrated using 6 2 integration   points per element such that the total number of integration points matched the total number of material points. 9 All analyses were conducted using MATLAB R2017b within a macOS V.10.14.5 environment on a 2.3 GHz Intel Core i7 with 16GB of RAM. For the coarsest analyses (h = 0.5 m) the material point run times are between 8% and 16% longer than that of the equivalent finite element analysis. The variation in the different material point methods is mainly due to the different number of global N-R iterations required to find convergence, as shown by the bracketed numbers in Table 4. As the problem is refined the run times for the finite element analysis and the material point methods increase at different rates. This is due to the structure of the implemented material point code being different from the finite element method, specifically that in the material point code the primary loop is over material points rather than elements in the finite element code. It is widely known that MATLAB is inefficient at handling large loops that operate on relatively small amounts of data within each loop execution. This is due to the overheads associated with RAM to cache communication (see [36][37][38], amongst others, for a detailed discussion on this issue, as well as methods to overcome it). Other differences include that for material point methods, the basis functions (and derivatives) need to be determined at each material point at the start of each load step whereas for the finite element method the basis functions can be determined once at the start of the analysis as the local positions of the integration points do not change. Table 5 provides a breakdown of the run time for the standard material point method with h = 0.5 m and 6 2 material points per element. The majority of the calculation time (almost 80%) is spent in determining the internal force and stiffness of the background grid nodes, (27) and (28), and within this the evaluation of the consistent tangent at each material point,Ã i jkl , takes over 30% of the total run time for the analysis. This can be attributed to the fact that these calculations must be performed for each material point (1,440 in total) for every N-R iteration (211 in this case) -303,840 material point evaluations in total. The other major contributor to the total run time is determining the interaction between the background grid and the material points in terms of element-material point searching, basis function calculation and determining a unique list of nodes for each material point. Updating the material points and assembling the external force are minor actions in terms of the run time, whereas the time spent in the linear solver is insignificant.

Elastic annulus under torsional deformation
The next example is of an elastic annulus with an inner radius of R i = 0.75 m and an outer radius of R o = 1.25 which was subjected to a torsional deformation field where the rotational deformation is dependent on the radial position, R, through where α controls the magnitude of the deformation. The problem is analysed by applying the appropriate body force 10 to cause this rotational deformation, determined via the Method of Manufactured Solutions (MMS). The inner and outer surfaces of the annulus are constrained by homogeneous Dirichlet boundary conditions. The problem is similar to that analysed by Kamojjala et al. [39] and Wang et al. [19]. The material had a Young's modulus of 1Pa and a Poisson's ratio of 0 and was modelled using the standard material point method using an unstructured background grid of linear three-noded triangles. The displacement field and associated body forces are shown in Fig. 18. Although the displacement field is relatively simple, the associated body forces required to cause this displacement field are complex. The problem was analysed with five different background grid sizes varying between 0.4363m and 0.0273m and 1 2 , 2 2 , 2 3 and 2 4 material points per element. A rotation of α = 0.1 radians was applied over 10 equal load steps. For the finest mesh this results in the material points at the midpoint of the annulus displacing 4 times the mesh size.
The convergence of the absolute displacement error with mesh refinement when analysed using the standard material point method is shown in Fig. 19. The convergence behaviour is similar to that shown for the compaction under self weight problem in that, provided that sufficient material points are used, the method approaches the correct convergence rate for the underlying basis. It is also clear that one material point per element is insufficient due to the complex variation in the body force over the domain. The figure also shows three of the background meshes used and the associated displacement magnitude of the material points at the end of the analysis with 2 2 material points per element.

Elasto-plastic collapse
The final example in this paper is the plane strain collapse of a 16m by 8m block of material with a Young's modulus of 1MPa and a Poisson's ratio of 0.3. Yielding of the material was governed by the von Mises yield criterion, (43), with ρ y = 20 kPa. Due to symmetry only half of the problem was modelled and roller boundary conditions were imposed on the base and the left hand edge of the background grid. The problem geometry and boundary conditions with h = 1 m and 2 2 material points per initially active background is shown in Fig. 20, where l 0 = 8 m. A body force of −10, 000N/m 3 was applied in the vertical, y, direction over 40 loadsteps. Fig. 21 shows the deformed material point positions coloured according to the vertical stress, σ yy , for different background grid sizes (h = 2, 1, 0.5 and 0.25m) and numbers of material points per initially populated grid cell (3 2 , 6 2 and 9 2 points per element). With coarse meshes, especially with low numbers of material points per initially populated element, there are severe stress oscillations through the body. These reduce with background and the material point refinement, however oscillations between adjacent elements still exist due to the combination of cell crossing and discontinuous spatial derivatives of the basis functions. Some authors, e.g. Wang et al. [40],  have combined the material point method with a plastic strain-driven softening constitutive relationship to model progressive slope failure. Although the distribution of the progressive failure segments appear to be physically realistic it is likely that they are triggered by the grid-aligned stress oscillations and, as with other localisation problems, will be highly mesh dependent. Table 6 gives the final extent and maximum height of the collapse using the standard material point method with different background grid sizes and numbers of material points per initially populated background grid cell. The extent and height converge with both grid and material point refinement. However, unlike with the elastic beam analysis the dominant factor is the grid spacing, h, and for finer meshes changing the number of material points from 6 2 to 9 2 has little influence on the final result.    (32) causes the domains to shrink artificially in the horizontal direction due to the vertical compression of the material. They also do not contract sufficiently in the vertical direction resulting in domains at the bottom of the problem moving beyond the background mesh. This causes a loss of partition of unity for those points and elevated the stresses at the base of the domain. As with the beam problem, if the domains are updated according to the normal components of the deformation gradient, (33), material points that experience large rotations spuriously shrink (as shown by the points in the dashed circle in Fig. 22 (b)). This results in a softer response compared to the other domain updating methods, as shown by the maximum extent and height values given in Table 7. This issue is removed by the volume-corrected deformation gradient update (38) or (39) which gives similar extent and  (34) and (37) domain updating methods as well as the CPDI2 method. Despite using generalised interpolation and CPDI2 material point methods, significant stress oscillations are still evident. Based on the findings of [26], it is unlikely that the oscillations are due to volumetric locking, instead they are a consequence of the bi-linear background grid. In the material point method the stress field is composed of the incremental summation of element contributions over a number of load (or time) steps. The variations in stress seen in Fig. 22 are due to the underlying basis of the finite element grid being unable to represent the stress field variation, rather than a manifestation of grid crossing. This is particularly evident on the right hand edge of the domain where the material points have moved, and crossed elements, horizontally. The stress variation within these elements clearly tracks to the background cell boundaries, with almost no variation in the vertical stress profile within the base elements which is consistent with the spatial derivatives of the background mesh.
The equilibrium convergence of the generalised interpolation material point method with h = 1 m and 6 2 material points per element is shown in Table 8 for loadsteps 10, 20, 30 and 40. Although in some cases the global iterations take a number of iterations to find the correct descent path, once approaching the correct solution asymptotic quadratic convergence is achieved. This verifies the consistent derivation and implementation of the consistent tangent for the previously converged Lagrangian formulation with large strain elasto-plasticity.

Observations
This section has presented the analysis of five large deformation problems using the previously converged Lagrangian formulation of the material point method. The examples have verified the formulation and its implementation through convergence towards analytical solutions and through obtaining the correct equilibrium rates of convergence, respectively. The following observations can be made from the analyses presented in this section: 1. Compaction under self weight: this simple large deformation problem highlighted issues associated with different generalised interpolation domain updating methods and, in this case, only (33), (34) and (37) are appropriate choices to achieve continued convergence. It has also shown that excessive numbers of material points are required to minimise the errors associated with cell crossing when using the standard material point method. 2. Simple shear: all of the material point methods predict the correct response for this problem. As the stress, strain, and material stiffness do not vary spatially within the problem domain there are no cellcrossing instabilities. The only method that did not predict the correct response was when the domains of the generalised interpolation material point method were updated according to the stretch tensor, U i j , through (34). 3. Elastic cantilever beam: results of this analysis are not highly dependent on the grid size provided that sufficient material points were used (above 6 2 material points per initially populated element). Low numbers of material points cause spurious stress oscillations through the beam. Applying the end load as a concentrated force on material points immediately above/below the neutral axis leads to an overestimation of the vertical displacement due to a local stress concentration. If the generalised interpolation material point domains are updated according to the normal components of the deformation gradient, F ii , through (32) the domains artificially shrink, due to large rotation components in the deformation gradient, and the method fails to complete the analysis. 4. Elastic annulus under torsional deformation: demonstrated the ability of the proposed formulation to be used on unstructured grids and reinforced the observation that the standard material point method will converge at the correct rate but only if sufficient material points are used.

5.
Elasto-plastic collapse: increasing the number of material points per background grid element reduces the magnitude of the stress oscillations, however local stress oscillations still exist. For the generalised interpolation material point method, updating the domains according to the determinant of the deformation gradient causes the domains to exceed the mesh due to insufficient contraction in the vertical direction, whereas updating the domains according to the normal components of the deformation gradient, F ii , results in an over-soft response due to spurious contraction of material point domains due to large rotation components in the deformation gradient.
The numerical examples presented in this paper have highlighted the importance of using sufficient material points per initially populated background element. However, we acknowledge that the use of 9 2 + material points per background grid cell is excessive for routine engineering analysis and we only include these results to demonstrate that errors associated with cell crossing can be mitigated via using a large number of standard material points. A more efficient approach is to adopt a domain-based material point method, such as the generalised interpolation or CPDI2 approaches and use fewer material points per element. The examples have also shown that only the volume-corrected deformation gradient domain update, (38)/ (39), and the volume-corrected corner update with mid-side determination, (37), lead to stable results for a range of physical problems that agree well with the more advanced CPDI2 method of Sadeghirad et al. [6]. Although the numerical examples in this paper have focused on two dimensional analysis, the formulations in this paper are general and suitable for three dimensional problems. There have been a number of papers that have analysed three dimensional problems with the material point method (see [7,26,41], amongst others). It is with these problems, where removing the need for remeshing whilst being able to model very large deformations, offers the most promise for the method.

Conclusion
Material point methods are unusual in that they are not ideally suited to traditional total and updated Lagrangian formulations of continuum mechanics. This is due to the fact that their basis functions are normally based on material point positions at the start of a load step and assume that calculations take place on the undeformed grid. In the case of updated Lagrangian formulations, this issue can be accounted for by mapping the basis function derivatives into the deformed configuration using the increment in the deformation gradient over the load step [21]. However this introduces an additional computational expense, and algorithm complexity, and the basis function derivatives must be individually mapped for each material point. In this paper we have proposed an alternative material point formulation where the equilibrium equations are satisfied at the nodal positions at the start of the load step. This previously converged Lagrangian approach predicts identical results to that of a conventional updated Lagrangian formulation but removes the need to map the derivatives of the basis functions. However, changing the equilibrium frame requires the definition of a new stress measure and also the derivation of an appropriate consistent tangent for large strain elasto-plasticity in order to achieve optimum convergence of the global equilibrium iterations. The new formulation is more computationally efficient than the conventional updated Lagrangian approach for the material point method and can be applied to all existing material point methods available in the literature. However, it should be noted that the primary focus of this paper is on implicit material point methods and the advantages of the proposed formulation are most evident within these methods. Explicit material point methods advance an analysis by forward-predicting from the previous state and for these methods the previously converged formulation reduces to the standard updated Lagrangian method.
This previously converged Lagrangian approach has been applied to the standard, generalised interpolation and CPDI2 material point methods and its performance demonstrated on a number of large strain problems with and without analytical solutions. In the case of the generalised interpolation material point method, different domain updating methods have been investigated and a new method proposed which avoids the degenerative properties of existing domain updates which can break down under pure rotation or shearing.
Although the focus of this paper has been on the material point method, the proposed previously converged Lagrangian formulation can also be applied to other numerical methods based on Lagrangian mechanics. and D alg i jmn is the small-strain algorithmic tangent obtained from the constitutive model. The derivative of the logarithm of the elastic Cauchy-Green strain tensor with respect to its argument is a special case of the more general formulation given by Miehe [42]. Note that all of the components of a i jkl should be evaluated in the spatial frame -the current deformed configuration.

A.2. Previously converged Lagrangian formulation
The consistent tangent modulus for the previously converged Lagrangian formulation can be expressed as The derivative of the Kirchhoff stress with respect to the increment in the deformation gradient is