iGIMP: An implicit generalised interpolation material point method for large deformations

The Material Point Method (MPM) uses a combined Eulerian-Lagrangian approach to solve problems involving large deformations. A problem domain is discretised as material points which are advected on a background grid. Problems are encountered with the original MPM when material points cross between grid cells, and this has been tackled by the development of the Generalised Interpolation MPM, where material points’ domains of influence extend beyond the currently occupied grid cell. In this paper, the Generalised Interpolation Material Point (GIMP) Method has been implemented implicitly in a manner that allows a global stiffness matrix to be constructed similar to that in the Finite Element Method (FEM) by combining contributions from individual elements on the background grid. An updated Lagrangian finite deformation framework has been used to ensure non-linear behaviour within each of the loadsteps. The weighting functions used for this which make the GIMP method different to standard MPM are presented and the implementation is explained. Specific details on computing the deformation gradient to be consistent with the updated Lagrangian framework and the updating of the material point influence domains are outlined, both of which are currently unclear in the published literature. It is then shown through numerical examples that for both small and large deformation elastic and elasto-plastic problems, the implicit GIMPmethod agrees well with analytical solutions and exhibits convergence properties between that of linear and quadratic FEM. 2017 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
When modelling continuum mechanics problems, specifically those involving large deformation, difficulties often occur when using mesh-based methods; such as displacement without remeshing being limited by mesh distortion. Due to this, there has been increased interest recently in particle-based methods, in particular the Material Point Method (MPM). The MPM models a problem domain as a collection of Material Points (or particles) which move through a fixed background mesh on which calculations are carried out. This offers an advantage over many other meshfree methods, which are also good for modelling large deformations and non-linear behaviour for example [1][2][3], because the existence of a background mesh removes the computational expense of undertaking nearest neighbour searches. In the MPM, properties are mapped between nodes on this background mesh and material points, during each load (or time) step. The majority of previous MPM research has looked at explicit formulations , with a few exceptions [47][48][49][50][51][52][53][54][55]. The advantages of adopting an implicit approach include allowance of much larger loadsteps, improved stability and error control, in comparison to explicit methods. An implicit formulation has also been shown in [49] to achieve superior accuracy. For static stress analysis problems, which are commonly tackled using an implicit Finite Element Method (FEM), there are benefits in an implicit MPM approach as there are many commonalities.
One of the main issues of the MPM is a well documented grid crossing instability. This occurs when material points move between elements in the background grid. Errors are introduced at the material points due to the non-continuous nature of the shape function derivatives between elements. These errors have been investigated in [56] and there are a selection of methods that have been proposed to address this issue and improve the MPM including CPDI [57,58], DDMP [59] and the Generalised Interpolation Material Point (GIMP) method [60]. Research into the GIMP method has also almost exclusively used an explicit approach, for example [60][61][62][63][64][65][66][67][68][69][70][71][72] with one notable exception [73] where the GIMP method is implemented implicitly using a matrix-free approach for dynamic problems using hyperelasticity. In an alternative approach [74] implicit MPM is used, however, rather than tackling the instabilities using GIMP, an additional non-physical stiffness term is introduced. Although this increases the numerical stability of the method, it destroys its ability to converge to analytical solutions.
In this paper, an implicit implementation of the basic GIMP method using an updated Lagrangian formulation is described using an approach that allows the local calculation of element stiffness matrices. The large deformation elasto-plastic implicit GIMP method described in this paper is implemented using an updated Lagrangian, geometrically non-linear formulation to accurately capture the behaviour of large strain problems. The formulation presented in this paper allows technology (such as constitutive models) to be easily shared between the implicit GIMP and implicit FEM and this is demonstrated by introducing von-Mises plasticity in some of the numerical examples. This is the first time that a fully implicit formulation has been proposed for an elasto-plastic GIMP method. The formulation includes the implementation of the spatial algorithmic consistent tangent to ensure optimum convergence of the global Newton process. This has required the calculation of the deformation gradient for implicit MPM methods to be clarified and a suitable domain updating procedure to be established. Additionally, we are able demonstrate the convergence properties of the GIMP method. The GIMP method is presented first before introducing the weighting functions (Section 2), outlining the finite deformation theory (Section 3) and describing the implementation (Section 4) before demonstrating the method using numerical examples (Section 5) and presenting conclusions (Section 6).

Overview of MPM methods
The MPM was first developed by Sulsky et al. [75,76] as an extension to solid mechanics of the Fluid Implicit Particle (FLIP) method [77,78], which itself was an extension to the Particle in Cell (PIC) method [79] used in fluid dynamics. The problem domain is divided into a number of background cells forming an Eulerian mesh through which the material points travel. Information from the material points is interpolated to this background mesh where the calculations are carried out. The new values are then mapped back to the material points and the material point positions are updated. This results in a method that combines advantages of both Lagrangian and Eulerian approaches, allowing boundaries and history dependent variables to be tracked while not encountering the problems associated with mesh distortion in large deformation problems in other methods such as the FEM. Since its inception, the MPM has been developed and improved as well as being used in a number of applications including granular materials [41,50,80,81], fracture [16,25,43,69,70,82] and geotechnical applications [20,21,37,36,42,46,54,55,63,68,72,[83][84][85][86].
In the standard MPM a problem can arise when a material point crosses the boundary between one background grid cell and another. This is due to the fact that shape function derivatives are not continuous between elements and this results in incorrect stresses being calculated. A significant advance in the MPM was made in [60] where the cell crossing instability was reduced with the GIMP method. To attempt to alleviate the problem, the GIMP method introduces a material point characteristic function describing the influence domain of each particular material point. This modification from the MPM means that it is possible that a material point can influence nodes other than those associated with the element it is inside. This occurs when the material point is close enough to the edge of an element that the domain overlaps adjacent elements. The introduction of the GIMP method is shown in [60] to give an improved stress response to the MPM. The errors introduced when material points cross element boundaries are reduced (although not completely removed [87]) because of the increased smoothness of the shape functions. Despite not being a complete remedy to grid crossing error the improvement shown by GIMP is significant.

Weighting functions
In the MPM, shape functions identical to those used in the FEM are used to relate nodal values to values at material points. For example, forces can be mapped to grid nodes from material points through where subscripts v and p refer to grid nodes (or vertices) and material points (or particles) respectively, n p is the number of material points in elements surrounding the node and N i are the standard shape functions as used in the FEM. A straightforward choice would be linear Lagrange shape functions given in 1D as where n is the local coordinate (in a domain À1,+1). These shape functions are also used to map properties from grid nodes to material points at the end of each step and their derivative are used to compute the stiffness matrix.
In the GIMP method [60], the standard FEM shape functions used in the MPM are replaced by weighting functions S vp which are constructed based not only on the FEM shape functions but also a material point characteristic function v p specifying the influence domain of the material point. This is the key difference between the GIMP method and the MPM, it can be thought of as giving each particle an associated domain rather than being a single point in a specific location. The weighting function ðS vp Þ can be calculated in a local coordinate system in one dimension (n) as where V p is the material point volume (or length in 1D), N v are the shape functions as shown in (2) with subscript v indicating values are associated with grid nodes (or vertices), X is the problem domain and X p is the influence domain of the material point. The gradient of the weighting functions ðrS vp Þ, can also be calculated using The original MPM can be recovered by setting the material point characteristic function equal to the Dirac delta function, that is In the GIMP method, the use of different functions for v p ðnÞ means that smoother weighting functions can be obtained. The simplest extension is to use a hat function with a value of unity within the material point's influence domain and zero elsewhere. This characteristic function, which is used in the development below can be expressed as It can be seen that a material point positioned at a node would not solely contribute to that node and would instead also have a small amount of influence on the surrounding nodes. Despite this the GIMP weight functions still possess partition of unity. To construct weighting functions in multiple dimensions, the tensor product of one dimensional functions is taken. This use of separate functions relating to material points has similarities to that of other meshless methods [88].
Using weighting functions that extend outside of an element presents a problem if one wishes to calculate element stiffness matrices in a manner similar to the FEM. In this case it must be ensured that although material points outside an element can have influence, the influence is only for the part of that material point's influence domain overlapping with the element. To address this problem, new weighting functions are constructed, where the integrations in (3) and (4) are only calculated over the area of each element. Fig. 3 shows (in a similar manner to the GIMP weighting functions) how the overlapping area between the material point characteristic function and the standard FEM shape functions within an element can produce new functions S vp a and S vp b that allow element stiffnesses to be calculated in a manner not previously possible. The weighting functions associated with the element a-b in Fig. 3 are where n 1 and n 2 are the integration limits for (3) in the local coordinates of the current element which can be given as where n p is the material point location and p is the material point domain size.
By summing these weighting functions at nodes from the contributions from different elements it is possible to recover the GIMP weighting functions as introduced earlier. This is illustrated in Fig. 4, where S vp b is reconstructed from contributions from elements a À b and b À c. It can be seen that these weighting functions mean that each node is not over-or under-accounted for.
The gradients of the weighting functions are calculated similarly using (4) and can be visualised as the overlap between the material point characteristic function and the gradients of the standard shape functions. Fig. 5 shows the gradient of the weighting functions within each element and the sum of these at node b. It can again be seen that these functions extend beyond the element but the gradients of the GIMP shape functions are recovered when contributions from both elements are considered. The area with a constant gradient is the section where the material point's influence domain is fully inside the element; at this point it is equal to the standard FEM shape functions. The gradient weighting functions for element a À b in Fig. 5 can be expressed as rS vp a ¼ n 1 À n 2 2l p and rS vp b ¼ where n 1 and n 2 are given by (8) and (9).

Geometrically non-linear GIMP
In this paper an updated Lagrangian finite deformation formulation is combined with logarithmic strain and Kirchhoff stress  measures that control the constitutive behaviour at each material point. This formulation is one of the most straightforward ways to implement large strain elasto-plasticity within a finite element framework [89]. In this framework all static and kinematic variables are referred to the previously converged state, rather than the original state in a total Lagrangian formulation. The majority of this section uses index notation to detail the geometric nonlinear formulation; only the discrete equations are expressed in matrix-vector form for convenience. The finite deformation framework adopted in this paper is based on implementations given in [90,91]. The framework has more generally been accepted and used by a number of authors including [92][93][94]. It is possible to extend this to allow plastic anisotropy using the formulation of Eterovic and Bathe [95] or to allow elastic and plastic anisotropy following the work of Caminero et al. [96] without modifying the overall framework. In this paper examples are restricted to isotropic elasto-plasticity for simplicity to ensure clarity of the GIMP method.
The weak form of equilibrium for an updated Lagrangian formulation can be expressed as where u t is the motion of a material body, B, subject to body forces, b i , and tractions, t i , on the boundary of the material domain, @B. The weak form is derived using a field of admissible virtual displacements g i . Within this statement of equilibrium, the deformation gradient is the fundamental variable that characterises the deformation at a material point where X j are the original reference coordinates and x i ¼ uðX i ; tÞ are the updated coordinates of the material point. It is assumed that the deformation gradient can be multiplicatively split into elastic and plastic components [97,98] where the superscripts e and p denote elastic and plastic terms, respectively. When implementing large strain elasto-plasticity there is a choice of which stress and strain measures to adopt. Here, we use the combination of logarithmic strains with Kirchhoff stresses as it allows the use of conventional small strain constitutive   equations within a finite deformation framework. In particular, the stress integration algorithm of plasticity models does not change provided that these stress and strain measures are combined with an exponential map of the plastic flow equation (see [99] or [100] for full details). All of the constitutive models used in this paper adopt a fully implicit stress integration algorithm based on backward Euler. This allows the updated stress state to be determined given a trial stress state (or trial elastic strain state) and the relevant constitutive parameters, again see [100] for full details.
Here we restrict this framework to the case where a linear relationship exists between the elastic logarithmic strains e e ij and the Kirchhoff stresses s ij that is where D e ijkl is the linear elastic isotropic material stiffness tensor. The elastic logarithmic strain is defined as where b e ij ¼ F e ik F e jk are the components of the elastic left Cauchy-Green strain tensor. The Cauchy stress can be obtained from the Kirchhoff stress using the relationship where J ¼ detðF ij Þ is the volume ratio. In order to obtain the current Kirchhoff stress state, s ij , the constitutive model requires a trial stress (or elastic strain state) to act as the initial estimate in the backward Euler stress integration algorithm. The trial elastic left Cauchy-Green strain tensor ðb e t Þ ij is obtained from where the subscripts t and n denote trial and previously converged states, respectively, rather than a physical index. DF ij is the increment in the deformation gradient for the current loadstep (see Section 3.2). The previous elastic left Cauchy-Green strain tensor ðb e n Þ ij can be obtained from the elastic strain state from the previously converged loadstep ðb e n Þ ij ¼ exp 2ðe e n Þ ij ð18Þ and the trial elastic strain is obtained using The updated Kirchhoff stress and the updated elastic strain states can then both be obtained from the constitutive algorithm.

Discrete implementation
Introducing the element approximation for the displacements at a material point where d i and d g i are the physical and virtual nodal displacements, respectively, a denotes the node number and n en is the number of nodes associated with an element. The Galerkin form of the weak statement of equilibrium over an element, E, can be obtained from (11) and (20) as where ½S vp the GIMP shape function matrix and ½G is the tensorial form of the strain-displacement matrix containing the derivatives of the GIMP shape functions with respect to the updated nodal coordinates. The first term in (21) is the internal force within an element and the combination of the second and third terms is the external force vector. Eq. (21) is non-linear in terms of the unknown nodal displacements and can be efficiently solved using the standard Newton-Raphson (NR) procedure. The nodal displacements within a load step, fDdg, are obtained by iteratively updating the displacements until (21) is satisfied within a given tolerance, that is where k þ 1 denotes the current iteration number, fdd kþ1 g are the iterative nodal displacements, ff R k g is the global residual out-ofbalance force vector (21) from the kth iteration and ½K is the global tangent stiffness matrix. The current displacement increment within a load step can be obtained through Linearising (21) with respect to the unknown nodal coordinates, and assuming that the applied body forces and surface tractions are independent of the nodal displacements, gives the element contribution to the global stiffness matrix as The non-symmetric spatial material tangent modulus of a material point is given by where D alg ijmn is the small strain algorithmically consistent tangent, that is, the tangent that is consistent with the adopted stress integration algorithm [101]. The use of this tangent allows for asymptotic quadratic convergence of the global residual (21). L ijkl can be determined as a particular case of the derivative of a general symmetric second order tensor function with respect to its argument; see Miehe [102] for details.
In material point methods, (24) is evaluated through the summation of the material point stiffness contributions where the nodal stiffness components of a single material point can be obtained through where V p is the material point volume in the spatial frame, that is V n p is the material point volume at the previously converged state, obtained from the product of the global influence domain lengths in the Cartesian directions. A material point's contribution to the internal force vector is given by Note that it is essential to use the volume in the spatial frame (28) in both (27) and (29) to ensure the correct order of convergence of the NR process. It may not initially be clear that detðDF ij Þ must be included to obtain the correct volume.

Deformation gradient calculation
One point of departure of implicit MP methods from conventional finite elements is the calculation of the deformation gradient. In an updated Lagrangian formulation, the deformation gradient is normally obtained through [99] F n ij is the deformation gradient from the previously converged state, Du i is the displacement increment in the current load step, x i are the updated coordinates (position in the spatial configuration determined from the nodal positions). However, in the MPM the concept of the current (or original) coordinates of the nodes does not exist.
The reason for this is that in MP methods the shape functions, and their derivatives, are defined assuming that the global coordinates of the background mesh remain in a regular grid. It is therefore not possible to use (30) 2 to determine the deformation gradient increment. Instead, the increment in the deformation gradient must be obtained using (see [99] amongst others) where e X i ¼ x i À Du i are the coordinates at the start of the load step.
The derivatives with respect to global coordinates at the start of a loadstep can then be obtained as where @ e X i @n j ¼ X nen a¼1 @ðNÞ a @n j ð e X i Þ a : Eq. (31) allows the determination of the increment in the deformation gradient based on a regular (undeformed) background grid. However, in order to form the stiffness matrix and internal force vector for an updated Lagrangian formulation we require the derivatives of the shape functions with respect to the current coordinates, x i . The mapping that links the current coordinate, x i , to the coordinate at the start of the load step, e that is, the inverse mapping of the increment in the deformation gradient (obtained from (31)). The derivatives of the shape functions with respect to the updated coordinates follow as where a is the node number. These derivatives are required for the construction of ½G, first seen in (21), as this matrix contains the derivatives of the basis functions with respect to the current nodal coordinates.

Domain updating
MP methods usually model a problem over a number of loadsteps and this presents an opportunity to update the influence domains of the material points at the end of each load step. Two ways of doing this labelled uGIMP and cpGIMP were presented by Wallstedt and Guilkey in [62]. uGIMP keeps the material point influence domains unchanged between loadsteps. This is the simplest approach to take, however it often results in domains overlapping each other or separating. cpGIMP addresses this problem by updating the size of the influence domain using the diagonal components of the deformation gradient. This rectifies the problem when deformation is aligned with the grid, however the method fails to improve matters when any rotation of the material occurs. Sadeghirad et al. [57] developed another approach known as the Convected Particle Domain Interpolation (CPDI) method which improves the MPM by updating the influence domains associated with material points. In the CPDI method, the initially rectangular material point domains are allowed to deform into parallelograms. Thus, the CPDI method is an extension to GIMP which removes a problem that exists when rotations occur for updating the material points. A further extension can be achieved by tracking the domain corners as shown in [58,73] which can ensure material points are contiguous. However, the CPDI method induces an additional approximation in the way that the basis functions of a material point are determined. Unlike in the GIMP method where the basis functions are determined analytically, in the CPDI method, to obtain the integration of the grid shape functions over the particle domain, a linear approximation between domain corner points is introduced. If all of the corners of a particle domain do not lie in the same element then errors are introduced into the basis function determination as the discontinuous nature of the grid shape functions is not captured by the linear approximation.
A simpler way of addressing the rotation problem is to use only the stretch part of the deformation gradient for updating domains in a cpGIMP fashion. In the cpGIMP method, it is necessary to update the material point domain lengths instead of updating the material point volumes as an independent variable In 1D this is the same, however in 2D or 3D it is important to take note of the changes in each direction. One option is to update the domain lengths, l p i , using components of the deformation gradient according to where l p 0 i are the original domain lengths. However, problems arise when the rotational component of the deformation gradient is nonzero [57]. Instead, here we propose a new approach where the domain lengths are updated according to the symmetric material stretch tensor where F ij ¼ R ik U kj and R ij is the rotational component of the deformation gradient. It should be clear from the above equation that the material stretch tensor is equivalent to the deformation gradient rotated back into the original reference frame. The material point domains can then be updated according to This updating is performed at the end of a load step once the NR process has converged. The important consequence of this rather minor modification to the theory is demonstrated numerically in Section 5.3.

Numerical implementation
When implementing the iGIMP algorithm, the first step is to discretise the problem into a set of material points spread over the domain, within a regular background grid which extends beyond the physical domain. A notable difference to the standard FEM is this requirement for the grid to extend to where material is expected to move into during a simulation. It is possible to keep track of how much of the grid is covered by material points and extend it if necessary. As in the standard MPM, the background grid is not restricted to any particular shape, however for convenience a regular mesh is usually selected. In this work, twonoded line elements are used in 1D and four-noded elements are used in 2D aligned with the coordinate axis. It is possible to use the same techniques in 3D. Elements initially containing material are populated with material points and a weight is assigned based on the volume of material represented by each material point. The influence domains are defined to initially cover the whole of the material with no gaps or overlaps. In this work, material points are arranged inside the elements in a uniform manner, however other initial positions can be chosen.
At the start of each load step the location of each material point with respect to the background grid must be determined. The local coordinates of the material point ðn; gÞ are calculated and, from these coordinates, the weighting functions (3) can be computed. Grid elements void of material points are also determined and not included in the calculation during the current load step. At this stage, any external forces on the material points should be incremented and then mapped to the grid nodes using The out-of-balance force is calculated as in (11) and the NR process is started. Displacements are calculated from the out-of-balance force and the stiffness calculated on the previous iteration (22). Afterwards, the strain displacement matrix and the derivatives of the displacement required for calculation of the deformation gradient can be calculated. It should be noted here that these quantities at the material points must be calculated as the sum of the different contributions from elements the material point overlaps, and particular attention should be made when calculating ½G to take into account the mapping outlined in (36). Due to material points potentially having influence over different numbers of nodes, the size of ½G can change between material points. The structure of ½G in 1D for a material point overlapping two elements is as follows where superscripts a and b refer to the derivatives of the weighting functions in elements a and b. The stress, internal force and stiffness can be calculated as shown in Section 3.1, and the out-of-balance force can be evaluated to determine whether the NR process has converged. Fig. 6. Implicit GIMP algorithm.
At the end of each load step, once the NR process has converged to within the designated tolerance, the material point positions and domains are updated and the background grid is reset. The algorithm is outlined in more detail in Figs. 6 and 7 showing the implicit Generalised Interpolation Material Point (iGIMP) algorithm in a similar way to [100]. It is possible to replace the material model in Fig. 7 with other models such as those presented in [96] to allow both elastic and plastic anisotropy, however for clarity the elasto-plastic model outlined in this paper has been included.

Numerical examples
In this section, three numerical examples are presented to demonstrate the iGIMP method. For each example, the geometrically non-linear updated Lagrangian GIMP method described in Section 3 is used.

One dimensional compression under self weight
The response of a column to the application of a body force due to increasing gravity, as shown using an explicit GIMP method in [60]. The column has an initial length (L 0 ) of 50 and is restrained at the bottom with uðz ¼ 0Þ ¼ 0. A total force of w ¼ 40; 000 is applied by assigning a density of q ¼ 80 and incrementing gravity (up to g ¼ 10). The Young's modulus is E ¼ 1 Â 10 6 , in compatible units. The analytical solution for the vertical Cauchy stress r in this 1D problem is now derived. The initial vertical position within the column is Z, therefore where q 0 is the initial density of the material and b is the body force. The Cauchy and Kirchhoff stresses are linked through (16). In one dimension, the logarithmic strain is defined as and we assume that the Kirchhoff stress is linked to the logarithmic strain through By combining the above equations, the Cauchy stress can be expressed as From (42) we can obtain r for any point in the problem domain and the deformation gradient can be found using a Newton process to Fig. 7. Implicit GIMP algorithm.
solve for F in (45). This analytical solution differs from the incrementally linear solutions given in [59,60] to be consistent with the fact that the method described in this paper is geometrically non-linear within a load step. For the case shown in Fig. 8, the domain is discretised into 50 elements with each element initially containing two material points positioned so that the influence domain of each material point consists of half the element, or V p ¼ 0:5. The stresses at the end of the simulation using the iGIMP method are compared against stresses calculated using the standard MPM using the same discretisation, and the analytical solution given below. The MPM and iGIMP simulations were both run using 20 loadsteps. It can be seen that the MPM simulation experiences an oscillation in its response, deviating significantly from the analytical solution, due to material points crossing element boundaries. In the iGIMP method this problem is alleviated as movements of material points between elements happen more gradually giving a smoother change in stiffness as opposed to a sudden jump. The effect of changing the element size in the background grid is now demonstrated. Both the number of background grid elements and the number of material points per element are changed and the error relative to the exact solution plotted both for the problem outlined above with a total load of w ¼ 40; 000 shown in Fig. 9(b), and with a load of w ¼ 10; 000 as shown in Fig. 9(a).
To aid comparison and to study convergence with refinement, a dimensionless error measure is used It can be seen in Fig. 9 that varying the number of material points does not have a large influence on the solution to this problem. The convergence rate varies between 1 and 2 for most numbers of elements with a degradation towards higher numbers of elements. It is possible that this can be attributed to the fact that there will be more material points crossing boundaries contributing additional error which cancels out the benefit of additional elements. The same problem is also modelled with a weight of w ¼ 400; 000, ten times larger than the initial problem to show the large deformation capabilities of the method. Fig. 10 shows the stress against position and the corresponding analytical solution, and Fig. 11 shows the convergence of the error with an increasing number of elements. Here, the convergence is also compared against linear and quadratic finite element solutions. It can be seen that for a given number of elements, the error for the iGIMP code is less than the linear FEM simulation with 2 Gauss points per element. The convergence rate for the linear FEM simulation is constant at 1 whereas the convergence rate for the iGIMP  simulations varies between 1.8 and 0.6. The convergence rate for the quadratic FEM code with 3 Gauss points per element is 2, again as expected. It should be noted that if 2 Gauss points per element are used then the FEM code achieves machine precision for any number of elements. This is because the two sampling points are correctly positioned to approximate the solution exactly. For a linear finite element the same applies to a single Gauss point in the centre of the element. In iGIMP the material points are not located at Gauss quadrature positions so the same does not apply.
From the analytical solution given in Eq. (45), the deformation gradient at the base of the column is calculated to be 0.74292 and the displacement at the top of the column to be À7.3347. Using two material points per element, and taking the top displacement from the top material point and bottom deformation gradient from the bottom most material point, it can be seen in Table 1 that the displacement is accurate to 5 significant figures for all numbers of grid elements shown, and the error in deformation gradient decreases with increasing elements with a linear rate of convergence, where the deformation gradient error is given as The error and rate of convergence are comparable to that of linear finite elements. The same problem of a column under self weight is now investigated but this time using a Von-Mises constitutive model with a deviatoric yield stress of q y ¼ 3 Â 10 4 . The yield surface is defined as Fig. 10. Stress against position for w ¼ 400; 000.
The material will yield when s zz ¼ q y which should occur at a posi- Below this point, the material will experience elasto-plastic behaviour and despite zero deformation being enforced in the out of plane directions, stresses s xx and s yy will be introduced. Due to fact that the boundary conditions are the same these two stresses will be equal, because of this from here on variables in the y direction will not be discussed. Using this it is possible to write the deviatoric stress in this situation as Here the stress in the vertical direction should remain the same, following the analytical solution (42) however in the section near the base of the column where this stress is reached, stresses will appear in the out of plane direction. This is shown in Fig. 12 where it can be seen that in the elasto-plastic region there are stresses in both the vertical and horizontal directions. The analytical solution for these out of plane stresses is not immediately obvious, using knowledge of the deformation due to the boundary conditions a relationship can be found between elastic parts of the deformation gradient It can further be shown that in the vertical direction a relationship between plastic and elastic components of the deformation gradient exists as Using this the Cauchy stress in the vertical direction can be written as This result allows us to calculate F e zz using a Newton process which yields F zz using (52). F e xx can also be calculated using (51) which allows the calculation of The derivation for this can be found in Appendix A.
The problem was repeated with a body force increased by a factor of 10. The solutions are shown in Fig. 13 to also show close agreement with the analytical solution. Fig. 14 shows the convergence for these simulations when refining the mesh which exhibits the same behaviour as for the elastic case. The convergence plot for quadratic finite elements is not shown in this case as it experienced locking with a 3 Â 3 quadrature scheme and reached machine precision in one step using reduced integration.

2D compaction under self weight
The second problem presented is the behaviour of a material compacting under self weight. The response of the material with increasing gravity is modelled. The problem domain at the beginning of the simulation has a height of 8 units and a width of 8 units, Young's modulus of E ¼ 1 Â 10 5 and Poisson's ratio of m ¼ 0:3, as shown in Fig. 15. Vertical movement along the bottom edge is prevented and due to symmetry, only half the problem is modelled. Self weight is applied incrementally over 20 loadsteps with a total weight of w ¼ 4 Â 10 5 . A 10 by 10 background grid is used to allow for material movement during the simulation. The initial position of the material is modelled using 4 material points per element (shown in a single shade of grey) along with the final (nonexaggerated) displaced shape of the material. The shading of the material point domains corresponds to the vertical stress at each material point at the end of the simulation. Fig. 16 shows that the convergence within a load step is near asymptotically quadratic. This convergence rate of the NR process indicates correct implementation of the method. The values shown in Table 2 correspond to the norm of the out-of-balance force Fig. 12. Vertical (rzz) and horizontal (rxx ¼ ryy) stress using iGIMP plotted against analytical solution for a load of 40,000 and yield stress of q y ¼ 3 Â 10 4 .
(ff oob g) at the end of each NR iteration for the final 5 loadsteps. This is calculated as outlined in Section 4. In Fig. 17, the maximum horizontal displacement of the material is compared against results obtained from an updated lagrangian Finite Element analysis of the same problem. The material properties are as described above. The FEM analysis used linear finite elements with 2 x 2 Gauss quadrature with 1,000,000 elements. By changing the number of grid elements, it can be seen that the displacement converges towards a constant value. Displacements have been normalised to the displacement given by the finest FEM mesh. When calculating the displacement using the iGIMP method, the value used is the displacement of the bottom right material point plus half of any extension to its influence domain. The convergence is investigated using different numbers of material points per element. With an element size of 1 and smaller, the solutions using different numbers of material points per element all come within 1% of the converged solution. This suggests that the number of elements in the background mesh has more influence on the accuracy of the solution than the number of material points per element; this can be seen clearly in Fig. 17 by the fact that for finer meshes the displacements for varying numbers of material points are all very similar.
The same problem is analysed with a Von-Mises constitutive model with a deviatoric yield stress of q y ¼ 1:2 Â 10 4 ; the yield function being defined in the same way as in (48). This leads to significantly larger displacements as can be seen in Fig. 18. The convergence for the final 5 steps is given in Table 3 where it can be seen that more steps are needed for the NR algorithm to find the correct path now that the material behaviour is elasto-plastic but then reaches near asymptotic quadratic convergence as before, until running into the precision of the machine for lower errors.

2D cantilever beam
The final example is an elastic cantilever beam with a point load of 100kN applied at the vertical mid point on its free end. To achieve this loading using iGIMP, this load is split between the  two end most material points above and below the neutral axis, as highlighted in Fig. 19. The beam has a length of 10 m and depth of 1 m, and the material has a Young's modulus of 12 MPa and Poisson's ratio of 0.2. The load is applied over 50 loadsteps with the problem initially split into 40 elements each containing 3 Â 3 material points. The beam is fixed at the left hand end in both directions at the neutral axis with roller boundary conditions applied to other nodes along the boundary. Fig. 19 shows the original and final (unexaggerated) configurations. Here, it is important to note that the material point influence domains are updated using the stretch rather than the full deformation gradient. The reasons for this are as discussed in Section 3 and it can be seen in Fig. 20(a) how the analysis collapses when using the deformation gradient for these updates (highlighted by the circled region on the right hand figure). Due to the material point rotations (a rotation of 90 degrees would cause the leading diagonal of the deformation gradient to go to zero) the size of the material points gets very small leading to an non-physically small stiffness in those elements.
In Fig. 21, the normalised horizontal and vertical displacement at the loading point are plotted against the analytical solution and results from a finite element analysis. For the iGIMP solution, this is the average between the two loading points above and below the neutral axis. The analytical solution is provided in [103]. The FEM analysis uses 8 noded quadratic elements with 3 Â 3 Gauss quadrature with the same 20 by 2 element discretisation as initially used in the iGIMP analysis. The load is applied at the neutral axis and the boundary conditions applied at the left hand end of the beam the same as the iGIMP analysis. Good agreement can be seen between the iGIMP displacements and both the FEM and the analytical solutions.

Conclusions
In this paper a fully implicit formulation of an elasto-plastic GIMP method is presented for the first time. The construction of the weighting functions and implementation of the method have been explained in detail. The element stiffnesses are shown to be calculated based on contributions from overlapping parts of material point influence domains allowing the global stiffness matrix to be assembled from the element stiffness matrices as in the FEM. In this paper, a new way of computing the deformation gradient and updated derivatives, which are consistent with the updated Lagrangian approach has been developed. Using the implicit GIMP method it has been shown that for both small and large deforma-tion elasto-plastic problems, accurate results can be achieved. It has been shown that by using these consistent values of deformation gradient and derivatives of shape functions when forming the consistent global stiffness matrix it is possible to maintain the correct convergence rate of the global equilibrium equations, and that within a load step the Newton Raphson process converges asymptotically quadratically as expected indicating a correct implementation. By increasing the number of background elements, the iGIMP method shows convergence properties between that of linear and quadratic FEM. Additionally, the novel use of updating material point influence domain lengths using the stretch tensor has been shown in Fig. 20 to allow the iGIMP method to model problems including rotation of material which previously has been problematic for the standard GIMP method when updating material point influence domains using the deformation gradient.

Appendix A. Analytical solution to column under self weight with plasticity
The response of a column to the application of a body force due to increasing gravity is modelled. The column has an initial length (L 0 ) and is restrained at the bottom with uðz ¼ 0Þ ¼ 0. Displacement is only permitted in a vertical direction. A Young's modulus of E and a density of q 0 in compatible units are assigned to give a total force, once gravity (g) is applied, of w. This time a deviatoric yield stress of q y is introduced. The yield surface is defined as f ¼ q À q y ¼ 0 ðA:1Þ as outlined in Section 5. The material will yield when s zz ¼ q y which should occur at a position of Z ¼ l 0 À  material will experience elasto-plastic behaviour and despite zero deformation being enforced in the out of plane directions, stresses s xx and s yy will be introduced. Because the boundary conditions are the same these two stresses will be equal, because of this from here on only variables in the x direction will be discussed. Using this it is possible to write the deviatoric stress as q ¼ js xx À s zz j.
The Cauchy stress in the vertical direction is the same as that for the elastic case can be determined from the initial vertical position within the column, Z, through r zz ¼ q 0 bðl 0 À ZÞ;

ðA:2Þ
where q 0 is the initial density of the material and b is the body force. The Cauchy and Kirchhoff stresses (s) are linked through r ¼ s J where J ¼ detðFÞ which in this case is equal to F zz . Due to the boundary conditions it is known that When there are only normal components, the logarithmic strain is defined as e ð0Þ ¼ The elastic part of the deformation gradient can then be found using a Newton process to solve for F e zz in (A.16). Using the above relationships it is possible to calculate the remaining parts of the deformation gradient and find the out of plane stresses as