An investigation of stress inaccuracies and proposed solution in the material point method

Stress inaccuracies (oscillations) are one of the main problems in the material point method (MPM), especially when advanced constitutive models are used. The origins of such oscillations are a combination of poor force and stiffness integration, stress recovery inaccuracies, and cell crossing problems. These are caused mainly by the use of shape function gradients and the use of material points for integration in MPM. The most common techniques developed to reduce stress oscillations consider adapting the shape function gradients so that they are continuous at the nodes. These techniques improve MPM, but problems remain, particularly in two and three dimensional cases. In this paper, the stress inaccuracies are investigated in detail, with particular reference to an implicit time integration scheme. Three modifications to MPM are implemented, and together these are able to remove almost all of the observed oscillations.


Introduction
For many years, the finite element method (FEM) has been the most used numerical technique to analyse and design structures, but it is well known that it is often unable to handle large deformations. In FEM, the geometry of a problem is attached to a mesh, and if the mesh suffers large distortions, the analysis is unable to continue. The material point method (MPM) is a numerical technique that overcomes this limitation [28,30], allowing for problems involving large deformations and multiple bodies to be analysed [5,29,33]. In MPM, the mechanical properties and the geometry of a problem are attached to a group of (material) points that move through an FEM mesh used to calculate the equation of motion in each time or load step. To enable this, the state variables are continuously mapped between the material points and the mesh. Many researchers have shown that MPM can be used to analyse some of the most common geotechnical problems, such as slope stability [1,4,9,15,31], foundation installation [17,21,25] and soil anchors [12]. However, the accuracy of MPM, in particular relating to the stress fields, is still far from the desired level. Indeed, it is noted that many publications do not display full results of the stresses, either presenting only deformations or limited data, and that the majority of work presented in the literature so far uses simple constitutive models. In some work, the stress oscillations and inaccuracies are acknowledged, and mainly attributed to the use of discontinuous finite element (FE) shape function (SF) gradients (e.g. [3,26,27,34]). Hence, a number of techniques have been developed to keep the SF gradients continuous between element boundaries, i.e. C1-continuous, for example: • GIMP [6], which distributes the influence of each material point over a characteristic or support domain, possibly extending the influence to multiple elements at a time. Both the SF and the SF gradients are modified.
• CPDI [24], which is an extension of GIMP in which the material point support domain can deform, maintaining the interaction between particles even after large extension. There are a number of CPDI variants, with different orientations and behaviour of the support domain. • B-spline MPM [26], which replaces the linear SFs by functions with higher-order B-spline basis functions, which are at least C1-continuous and positive definite. • DDMP [34], which preserves the linear SFs and replaces the SF gradients by smooth continuous functions, thereby allowing the usage of a local integration procedure rather than having a material point support domain.
These techniques have been proven to reduce the impact of cell crossing. Meanwhile, other techniques use material point integration together with Gauss point integration to reduce numerical inaccuracies [2,16]. However, a complete investigation of the causes of the stress inaccuracies has not been presented. Moreover, these techniques typically involve explicit MPM schemes, thereby ignoring the errors the proposed solutions can cause in the integration of the stiffness matrix in implicit schemes and not exploiting the advantages of implicit time integration. Finally, examples have often been investigated only for 1D cases (usually with 2D elements), so that oscillations caused by other deformations, e.g. material rotation or distortion, have not been examined.
This paper first presents the theoretical background of implicit MPM. Then, two benchmark problems are introduced to illustrate the oscillation problem. In Sect. 4, the main causes of stress oscillations are investigated. Then, a series of existing and novel solutions are presented and investigated. Finally, a comparison of regular MPM and the new proposed oscillationfree MPM is given for the simulation of a vertical cut failure, in order to demonstrate the relative performance for a problem involving both 2D geometry and elasto-plasticity.

Theoretical formulation
MPM shares the same continuum mechanics background as FEM. The equation of conservation of momentum is given as where σ is the symmetric Cauchy stress tensor, a is the acceleration, b are the body forces, and ρ is the mass density. In MPM, because of the partition of unity of the SFs, mass is automatically conserved. The weak form of Eq. 1 including traction as a boundary condition is where u denotes the virtual displacement, τ is the traction at the surface Γ (i.e. the boundary condition), and V is the body volume. Following standard FEM discretisation, Eq. 2 can be expressed in the matrix form [8,32] where M is the mass matrix, a is the vector of nodal accelerations, K is the stiffness matrix, u is the vector of nodal displacements, and F ext and F int are the external and internal force vectors, respectively. A quasi-static formulation is obtained by removing the Ma term from Eq. 3. Using the Gauss-Legendre quadrature rule and discretising the continuum body using a finite set of material points, the element (nodal) mass matrix can be expressed as where nmp is the number of material points in the element, ρ ρ is the material point density, N is the matrix of SFs evaluated at the material point position x p , J is the Jacobian matrix, and W is the material point integration weight (which is dimensionless and equal to the volume of the material point in local coordinates).
The element stiffness matrix K can be expressed in terms of a small or large strain formulation, but for simplicity it is expressed here in the small strain formulation (for details of the large strain formulation, see [32]), as where B is the strain-displacement matrix and D p is the elastic matrix at the sampling point. The element (nodal) external forces F ext considering gravity and boundary tractions are where g is the gravity vector. The element (nodal) internal forces F int are where σ p is the vector of material point stresses. Details of the axisymmetric form of the previous equations are presented in "Appendix A1". Using Newmark's [20] time integration scheme, where Δt is the time step, t+Δt , t+Δt and t+Δt are the respective vectors of displacements, velocities and accelerations at time t +Δt, and α and δ are time stepping parameters that are chosen to be α = 0.25 and δ = 0.5 to give a constant-average-acceleration. Substituting Eq. 9 into Eq. 3 and rearranging leads to where Δ = t+Δt − t is the vector of incremental displacements. In "Appendix B", a study of the conservation of mass and momentum of the implicit MPM is analysed. The above equation governs the behaviour of the body and it is therefore important to accurately evaluate each of the terms in order to ensure realistic behaviour. Following the solution of the updated displacements, the trial incremental stresses at the material points can be computed using the strain-displacement matrix as For an elasto-plastic material, stresses which are found to exceed the yield surface are redistributed using a consistent plastic return algorithm such that a new body force is calculated, and Eq. 10 is again solved to give plastic deformations. This is iteratively performed until no stresses exceed the yield surface. For more details see, for example, Bathe [8].

Material point method
MPM discretises the material into a series of (material) points which carry the information of the material (density, mass, deformation, velocity, acceleration and stresses). A mapping phase occurs at the start and at the end of each time or load step. At the beginning of each step the values required at the nodes in Eq. 10 (velocity, acceleration, etc.) are mapped via the SFs (see Wang et al. [32] for details). The matrices required for the calculation are calculated via element integration, as shown in Sect. 2, typically using material points as the sampling points. Afterwards, element assembly results in a set of global matrices representing nodal equations. A finite element calculation is then performed, with the state variables calculated at the nodes, in order to compute the deformation of the domain. Finally, another mapping step is undertaken to update the position and state variables of the material points. In Fig. 1, a sketch of the steps followed in MPM is shown. The SFs used to carry out the mapping and the integration are usually first order (e.g. bi-linear in 2D) to avoid negative values which cause instability.

Benchmarks
Two benchmarks are introduced to demonstrate and investigate the inaccuracies which occur in MPM. The first benchmark consists of an elastic quasi-static axisymmetric problem. The second benchmark is a 2D dynamic, elastoplastic, vertical cut problem.

Axisymmetric benchmark
The first benchmark is similar to that presented by Naylor [19] and Mar and Hicks [18] to explore stress recovery. It consists of a hollow cylinder which deforms due to an incremental pressure (Δp s ) applied on the internal boundary (s). The main benefit of this benchmark is that, unlike a 1D plane strain problem, the stresses inside the elements are not constant; moreover, they deviate from the real solution and, depending on the material point position, the deviation may be large or small. Figure 2a, b shows the initial conditions of the benchmark; that is, the top view of the cylinder and the finite element discretization of the cylinder wall, respectively. In both figures, the position of the boundary material point is shown (i.e. the material point nearest to the cylinder axis), which is used to determine the position of the boundary (s). Figure 2c and d illustrate that, during the loading, the distance r i to the inner wall (s) changes, and is equal to the distance between the cylinder axis and the nearest active node (this implies that r i remains constant until the boundary material point jumps to the next element). To enable the numerical (large strain) solution to be interpreted in terms of the analytical (small strain) solution, the methodology includes the following three features: (1) the applied pressure Δp s on the (2) due to the new location of the inner wall, Δp s is re-evaluated as Δp s (r i ) = A/ r i 2 + 2C, where A and C are constants associated with the initial geometry and boundary conditions of the benchmark, as shown in Fig. 3 (a description of the analytical solution and the constants A and C are presented in "Appendix A2"); (3) instead of accumulated stresses, incremental stresses at the material points are used throughout the analysis. These three features ensure that the incremental stress at the material points, for an arbitrary position of the cylinder wall, can be compared to the analytical stress related to the original geometry of the cylinder.
The inner (initial) and outer cylinder boundaries are located at r i = 0.5 m and r e = 1.5 m, respectively. The cylinder domain is discretised by elements of dimension Δr = Δy = 0.20 m, and each element initially contains four material points equally spaced. The elastic properties are Young's modulus, E = 1000 kPa, and Poisson's ratio, ν = 0.30. The initial applied pressure increment is Δp s = 100 kPa, and A and C are 19.56 kN and 10.87 kPa, respectively.
In Fig. 4, the incremental stress invariants (deviatoric stress Δq and mean stress Δσ m ) at material point mp 1 are plotted and compared to the analytical solution over 25 Δp s increments. It is evident that the stress invariants can deviate strongly from the analytical solution.

Vertical cut benchmark
A 2D elasto-plastic vertical cut problem has been simulated using the Von Mises constitutive model incorporating postpeak softening as described in Wang et al. [32]. Figure 5 shows the domain, boundary conditions and discretisation. The height H of the cut and length L of the domain are 3.0 m and 6.0 m, respectively; the element size is Δx = Δy = 0.10 m and each element contains initially four equally distributed material points. The elastic parameters are E = 1000 kPa and ν = 0.35, whereas the peak cohesion is c p = 12 kPa, the residual cohesion is c r = 3 kPa, and the softening modulus is H s = − 30 kPa. At the left boundary, the nodes are partly fixed to avoid displacement in the horizontal direction, whereas the nodes are fully fixed at the bottom boundary. The initial stresses in the domain are generated by fixing the locations of the material points and applying gravity loads until the internal and external forces are in equilibrium. After equilibrium is reached, the material points are released and deformation takes place. Figure 6a, b shows contours of the deviatoric and mean stresses, respectively. It is seen that during the movement of material points, both deviatoric and mean stress oscillations occur, although the overall failure mechanism is as expected. For Fig. 6b, the shown range was fixed between 10 and − 30 kPa because the oscillations are enormous in and around the shear band.

Oscillations in MPM
The MPM technique can be seen as an FE stepwise procedure, in which the integration points (now called material points) move together with the mesh, but keep their new positions while the mesh returns to its original position. This allows the simulation of large deformations since extreme distortion of the mesh is avoided, although the process is found to cause stress oscillations. There are a number of contributing factors causing these oscillations, which are investigated below.

Stress recovery
As is typical in many implicit FEM formulations, displacements have been used as the primary variable and stresses are back-calculated using the strain-displacement matrix and the elastic matrix (Eq. 11). During the back calculation of stresses, an oscillation occurs because the stresses inside the elements, interpolated using the element SF gradients, do not agree with the analytical stresses except at the superconvergent positions [7,19,35]. This problem is not observed in problems where the analytical stress is uniform across the element, e.g. as in a 1D bar. Figure 7 illustrates the radial stress inside a linear or quadratic axisymmetric element. It is seen that the computed stress distribution across the linear element (σ L ) is different from that across the quadratic element (σ Q ), and that both are different from the analytical stress (σ A ). However, the linear and quadratic stresses (σ L and σ Q , respectively) match the analytical solution exactly at the Gauss point locations. This means that, depending on the position of the material point, the recovered stresses can be either higher or lower than the analytical stresses, as illustrated in Fig. 4. Figure 8 shows the analytical radial stress distribution and the stress recovered using MPM (or FEM) at any stress recovery position for the first load step in the axisymmetric benchmark. It is evident that the exact solution is near the centre of the elements, and recovering stresses at any other position will cause oscillations. It can also be seen that there will be a large oscillation whenever a material point crosses an element boundary, since the radial stress is discontinuous across inter-element boundaries.

Nodal integration using SF gradients
The nodal integrations of F int and K are performed using SF gradients and the material point positions. However, considering that the SF gradients used in MPM are bi-linear (linear elements) and discontinuous, and that the material point positions change each time step, the resulting nodal values are inaccurate, especially if material points cross element boundaries. Next, a description of the SF gradients in MPM and the consequences of using them are presented.   Figure 9 shows the SF (Fig. 9b) and the horizontal and vertical SF gradients (Fig. 9c, d) of node 1 of a 4-node square element (Fig. 9a). It is noticed that the SF gradient is a maximum at the node, constant in the direction associated with the SF gradient, and decreases down to zero in the orthogonal direction. When a material point crosses an element boundary, the combination of the two element SFs must be considered. In Fig. 10, two elements are shown: E 1 and E 2 (Fig. 10a). The SFs and SF gradients in both directions of node 5 are shown in Fig. 10b-d, respectively. Figure 10b shows that the SFs are continuous between elements, while Fig. 10d shows that the vertical SF gradient is continuous between elements in the horizontal direction and constant in the vertical direction. On the other hand, Fig. 10c shows that the horizontal SF gradients at the inter-element boundary are discontinuous, and that they decrease in the vertical direction.

Integration of the internal forces F int and stiffness K
Using SF gradients in the integration of any variable (i.e. F int and K) results in an inadequate nodal distribution, whereas, if regular SFs are used, the nodal distribution is smoother (M and F ext ). Moreover, two differences should be noticed between the integration of F int and K. The first is that, to integrate F int , the strain-displacement matrix (B) is used once (Eq. 7), whereas the element stiffness is computed using both B and its transpose B T (Eq. 5). The second is that to integrate F int , the stresses of the material points are used, whereas to integrate K the elastic properties of the material points are used. The significance of this is that the elastic properties are constant throughout the analysis, whereas the material point stresses change during the analysis, causing possible accumulation of errors.
As an example of the inaccuracies caused by using SF gradients, the vertical and horizontal nodal internal force distributions ( F int x and F int y ) and the diagonal entries of the stiffness matrix (Eq. 5) corresponding to the vertical and horizontal degrees of freedom (K x and K y ) using two different material point distributions, are computed for nodes 1-5 of the plane strain finite element mesh shown in Fig. 11. In both cases the material points are equally distributed inside the elements; in the first case (Fig. 11a) the material points are located inside each element, whereas in the second case (Fig. 11b) the material points have moved and some are located at the inter-element boundaries. After the movement, the material points are still located inside their original element, except for material points a-d which have crossed the boundary by an infinitesimal distance. Stress components of σ x = σ y = − 1.0 MPa and σ xy = 0, a Young's modulus of E = 1.0 kPa and a Poisson's ratio of ν = 0 for each material point have been considered, while the distance between the nodes is 1 m and the material point weights are equal to 1.
In Fig. 11c, d, the vertical internal force is equal to zero in both cases. The force is unchanged because the horizontal displacement of the material points does not affect the values of the vertical SF gradients, and equals zero because the internal vertical forces on both sides of the nodes are the same but with an opposite sign. However, the distribution of the horizontal internal force is highly inaccurate due to the material point crossing the element boundary and the discontinuity of the horizontal SF gradients (Fig. 11d). When integrating the nodal stiffness, the horizontal and vertical stiffnesses are initially similar (Fig. 11e). However, as the material points cross an element boundary (Fig. 11f), the inaccuracies are evident again, although they are smaller than those of the internal forces. This is because the product BB T returns positive nodal values, so avoiding the change in sign of the SF gradients.

Nodal integration of the mass M and external forces F ext using SFs
The integration of M and F ext is performed using SFs rather than SF gradients, so that discontinuities between elements do not occur. In this example, only the external forces caused by gravity are considered. Since a lumped form of the mass matrix is used, and also because of the partition of unity of SFs, any initial distribution of material points inside the elements results in the same nodal mass (or external force), as long as the distribution is symmetrical. As an example, Fig. 12 shows two different material  Figure 13 shows the distribution of M for the same problem as in Fig. 11. It is clear that the movement of material points and the crossing of nodes does not cause any trouble for the nodal integration because of the continuity of the SFs. Also, since the integration of F ext is performed in a similar manner to M, the distribution would be similar to the one in Fig. 13.

Plastic stress redistribution
The stress oscillation caused by the plastic stress redistribution is an extension of the oscillations explained in the previous sections. As the stresses exceeding the yield surface are integrated as a new external force computed with SF gradients, additional oscillations comparable to the F int oscillations are introduced. Moreover, oscillating stresses could cause some points to yield spuriously, leading to an unrealistic system behaviour.

GIMP
The generalised interpolation material point (GIMP) method [6] was proposed to reduce oscillations derived from material points crossing element boundaries. In GIMP, FE SFs are replaced by functions constructed based on the linear FE SF and a material point support domain (SD). This means that each material point has a domain over which its influence is distributed. The GIMP SF (S ip ) and its gradient (∇S ip ) in one dimension are computed as where V p is the material point volume, Ω is the problem domain, Ω p is the material point support domain, i is the node, and χ p is the characteristic function delimiting the area of influence of the material point and is given as The support domain is often assumed to be square, with a size of 2lp (lp = half of the material point support domain), which is obtained by dividing the element size by the number of material points. In Fig. 14, a 1D comparison between an FE SF and a GIMP SF is plotted, considering a distribution of two equally-distributed material points per element. It is seen that the GIMP SF and GIMP SF gradients are no longer exclusive to a single element and that the GIMP SF gradients are continuous between elements.
The GIMP SFs in 2D and 3D are computed as products of the 1D GIMP SF in each direction; that is, where S k ip is the 1D GIMP SF in the k-direction. An additional advantage of including a support domain is that the material boundary is explicitly defined, and can be used to apply boundary conditions.

Modified integration weights
To reduce the problems caused by an irregular number of material points inside an element, it is here proposed to modify the material point integration weight to (13)   where W* is the modified material point weight (dimensionless), cmp is the current number of material points in the element, and omp is the original number of material points in the element. This modified weight is used considering only structured meshes, i.e. a mesh composed of equal-sized square elements, and equal mass material points, and its use with unstructured meshes or unequal mass material points is not part of this work. This modified weight technique differs from the approach of other researchers who have modified the weights based on volumetric strain (e.g. [11]), which, while compensating for 1D deformations of the material points (compression or extension), does not reduce the problems caused by the rotation or advection of the material points. Finally, it should be noted that for four noded elements this modified weight value reduces to 4.0/cmp.

Double mapping (DM)
Integration using SF gradients is seen to work only at Gauss point locations, whereas material point integration is stable when based on SFs. Therefore, mapping to the Gauss point locations using shape functions (via the nodes) is proposed.
As an example, the stiffness matrix is used. The elastic matrix is mapped to the nodes from the material points and then to the Gauss points, prior to the integration. Using FE SFs, the material point elastic matrix is mapped to the element nodes as were D i is the elastic matrix at node i, and D p is the elastic matrix of material point p. At this point, the total stiffness contribution of the material points is accumulated at the nodes, and this contribution is then redistributed to the original Gauss positions as were D g is the elastic matrix at the Gauss point, N i (x g ) is the nodal SF evaluated at the Gauss points, and nn is the number of nodes of the element. By substituting Eq. 16 into Eq. 17, D g is obtained as Finally, combining Eqs. 18 and 5 (in FEM form) results in the nodal stiffness: where ngauss is the number of Gauss points in the element and W FE is the weight associated with Gauss point g (as in FEM).

DM-GIMP (DM-G)
As mentioned in Sect. 5.1, the GIMP method was created to avoid problems caused by the use of discontinuous FE SF gradients. However, a simple example in calculating the stiffness reveals a key problem. Figure 15 shows the same problem as in Fig. 11, but in this case the stiffness is computed using regular SFs and GIMP SF gradients.
As shown in Fig. 15a, for the initial configuration of material points, the nodal stiffness distributions remain the same for both techniques, because at this position the MPM and GIMP SFs and SF gradients are the same. With the movement of the material points (Fig. 15b), the nodal stiffness computed with GIMP decreases, as the GIMP SF gradients drop to zero at the inter-element boundaries (as shown 14 a GIMP shape function (S ip ) and regular FE shape function (N i ) of node i, and b GIMP shape function gradient (∇S ip ) and regular FE shape function gradient (∇N i ) for node i in Fig. 14). In addition, the contribution of material points in neighbouring elements is not capable of compensating for this drop. This would be the case for other methods, including DDMP and CDPI, that have this same characteristic.
To overcome the problems of using GIMP to integrate nodal stiffness, it has been proposed that the double mapping approach be used alongside the local GIMP SFs [10]. The local GIMP SFs (S ip* ) are similarly created as regular GIMP SFs, but the influence of the material point support domain affects only the nodal FE SF in a single element rather than contributing to all contiguous elements. In Fig. 16, an illustration of the development of regular and local GIMP shape functions of a node is shown.
In a similar manner to the double mapping technique using regular SFs, by using local GIMP SFs it is possible to distribute the elastic matrix to the nodes of an element and afterwards to the Gauss positions. The element stiffness matrix is then constructed as where S ip* is the local GIMP SF of node i evaluated at the material point position, and smp is the number of material points with a support domain inside the element. The algorithm to compute the stiffness matrix using DM and DM-G is given in "Appendix C", together with a study of the computational performance.

Composite material point method (CMPM)
The composite material point method (CMPM) [14] is a modification of the composite finite element method (CFEM), proposed by Sadeghirad & Astaneh [23], in which the support domain used to recover the stresses is extended, i.e. a patch, improving the accuracy of the stresses computed. New shape (20)  functions enveloping all neighbouring elements of the element containing the material point are developed using Lagrange interpolation. In Fig. 17, the C 2 shape functions are shown in 1D, in which each shape function N 2 envelopes the local element plus the neighbouring elements.
Using Lagrange interpolation, each of the N 2 shape functions is computed as where ξ is the nodal local coordinate in the extended domain, n is the number of nodes, ξ j is the local coordinate of the N i 2 (21) shape function, and ξ m is the local coordinate of the remaining nodes. Solving Eq. 21 for each node, the CMPM shape functions for an element with two neighbours are If the material point is located at the boundary, as in Fig. 18, the CMPM shape functions are then It is important to mention that although the CMPM SFs extend beyond the limits of an element, the range of the functions remains between − 1 ≤ ξ ≤ 1. Also, this solution can only be used with a structured mesh. To extend the solution to a 2D domain, the new SFs are the product of the SFs in each direction. Finally, trial stresses using CMPM are computed as where ∇ 2 is the matrix of the CMPM SF gradients, and ext is the vector of nodal displacements in the extended domain.

Testing of proposed techniques
Since the novel techniques presented in this paper are designed for the integration of the stiffness, the testing performed in this section is focused on the stiffness matrix. To compare the stiffness using each technique, the stiffness magnitude is used, and this is computed as The test consists of computing the stiffness of an infinite space made up of square elements that are full of equally spaced material points, four per element, as shown in Fig. 19a. The infinite domain is then rotated 20° degrees around its centre (C), as in Fig. 19b. The elastic properties of the material are E = 1000 kN/m 2 and ν = 0.30. Figure 20 presents the stiffness computed using regular MPM and DM and the results are compared with the FEM stiffness, computed using four Gauss integration points (K mag = 3263.57 kN/ m 2 ). In addition, the stiffness using the modified integration weights (W*) and Gauss mapping (GM) separately (the two components of DM) are shown to highlight their comparative effects. Since the material points remain equally distributed after rotation, the stiffness of the domain should not change (i.e. be mesh independent). Finally, a further test is performed using two materials, by considering the properties of material points below line A-A′ to be E = 1500 kN/ m 2 and ν = 0.25.
Theoretically, the stiffness of the domain should be independent of the rotation of the field of material points, and should be equal to the FEM stiffness before rotation (for the case with one material). As can be observed in Fig. 20, the stiffness obtained using regular MPM is not accurate and improvements are needed. After including the modified integration weight (W*), which accounts for a varying amount of material points per element, the stiffness distribution oscillates, although with a different spatial pattern than in regular MPM. Using GM the oscillation also persists, as the number of material points per cell is still incorrect, but it is less than in regular MPM because it helps to reduce errors due to material point position. It is noted that including W* or GM separately is unable to fix the stiffness oscillation, and that the spatial distribution is almost opposite in pattern, i.e. where high values occur in GM, low values occur in W*, and vice versa. Using DM, i.e. combining GM and W*, the stiffness oscillation is reduced significantly, as it accounts for both the material point position and the number of material points per element. Moreover, the transition is smooth over the elements when two materials are used. In Fig. 21, the tests from Fig. 19 have been performed using GIMP and DM-G. As can be observed, the stiffness obtained using GIMP integration is significantly more inaccurate when compared to MPM integration, as it both oscillates and reduces in magnitude. Note that the results for GIMP are shown using a different contour range; this is because using GIMP SF gradients the stiffness reduces significantly, and it is necessary to change the contour range to visualize the stiffness distribution. On the other hand, using DM-G the stiffness oscillation is reduced further than using DM. This is because the W* approach, which only allows the impact of a discrete number of points in each element to be considered, is not being used. Utilising DM-G allows a gradual transition of mass from one element to another. Moreover, using DM-G, the transitions between materials appears sharper than in regular FEM due to an increase in the accuracy of the material stiffness distribution between the interface nodes.
In Table 1, the difference between the stiffness obtained using each technique is shown relative to the nodal stiffness magnitude of the real FEM stiffness. In this comparison, only the homogenous material is considered. As can be observed, regular MPM and GIMP give large stiffness oscillations relative to the FEM stiffness, but in the case of GIMP the stiffness only decreases (as observed also in Fig. 15b). Using only the modified integration weight the stiffness oscillation increases, whereas using the GM stiffness the oscillation decreases (compared to regular MPM), but not significantly. Using DM and DM-G, the dependence between the mesh and the position of the material points is reduced, and the nodal stiffness oscillations reduce significantly, especially using DM-G where the variation is smaller than 1%.

Benchmark problems including improvements
The benchmark problems introduced in Sect. 3 are now reanalysed using the improvements described in Sect. 5.   Figure 8 showed the stress oscillation caused by using regular SFs to recover stresses in the cylinder wall. In Fig. 22, GIMP and CMPM are compared against regular MPM for a single (i.e. the first) load step. As can be seen, the GIMP oscillation is the same as MPM close to the centre of the element, because there the SF gradients are the same for both techniques. However, stresses are continuous between the elements, due to the continuous gradients of GIMP. On the other hand, using CMPM the stresses remain discontinuous between elements, but the reduction of oscillation when compared to the analytical solution is significant. In Fig. 23, the evolution of the incremental deviatoric and mean stresses of material point mp 1 (over 25 load steps) are shown, comparing the stresses obtained using normal MPM (as shown in Fig. 4), DM and DM-CMPM (DM-C). As can be seen, there is a significant increase in the accuracy of the stresses recovered using the DM technique, due to the improved stress recovery and stiffness integration. Moreover, if CMPM is included in the analysis, the stress oscillation reduces still further to give stresses close to the analytical solution. Next, the same example using DM-G and DM-GIMP-CMPM (DM-GC) is studied. Using DM-G, the stiffness is computed with the DM-G method and the stresses are recovered using GIMP SFs. Using DM-GC, the DM-G method is again used to compute the stiffness, but the stresses are now recovered using CMPM rather than GIMP. In addition, since the inner wall boundary can be determined accurately using the material point support domain (as mentioned in Sect. 5.1), the distance between the cylinder axis and the inner boundary (s) is r i = r mp1 − lp as in Fig. 24. Then, the applied pressure Δp s is distributed linearly to the nodes of the boundary element based on proximity.

Axisymmetric benchmark
In Fig. 25 it can be seen that, using DM-G and DM-GC, the results approximate the analytical solution even better than DM and DM-C, respectively. This is because the stiffness computed using DM-G is closer to the FEM stiffness and also due to the accurate distribution of the external pressure considering the accurate location of the internal boundary.  Figure 26 shows the elastic stiffness magnitude in the vertical cut benchmark problem, using regular MPM and DM-GC. As can be observed in Fig. 26a-d, using regular MPM large stiffness oscillations occur, from the beginning (small deformations) up until the end (large deformations) of the analysis. In contrast, using DM-GC (Fig. 26e-h) the stiffness oscillation reduces significantly, although some small oscillation can be observed in the shear band and along the edge of the domain. In Fig. 27 the nodal F int magnitude is shown, once again comparing regular MPM and DM-GC. Analogous to Eq. 25, the magnitude of the nodal internal force is computed as It is seen that if GIMP and CMPM are included in the solution, a large reduction in the oscillations of F int is obtained. Using GIMP, the oscillation caused by the material points crossing cell boundaries is reduced. Furthermore, by including CMPM, the recovered stresses are improved, reducing the oscillation caused by the stress recovery position. Figure 28 shows the deviatoric stress contours from both analyses. It is evident that, after reducing the oscillation in the stiffness and the internal nodal forces by using DM-GC, the deviatoric stress distribution in the domain is significantly smoother. Similarly, Fig. 29 shows the comparison of mean stresses during the analyses, demonstrating that the In this case, as in the axisymmetric benchmark, some oscillation of the mean stresses still occurs, but this is thought to be due to incompressibility during plastic yielding. For methods to reduce locking behaviour in MPM using low order shape functions the reader is referred to Coombs et al. [13]. As can be seen from previous figures, the oscillation of material point stresses, nodal stiffness and internal nodal forces are reduced significantly using DM-GC. Plots for the nodal mass and external nodal forces are not included in the results, since the oscillation for both MPM and DM-GC is small.

Vertical cut benchmark
Finally, p-q curves have been plotted for 3 material points at key positions in the soil body. Figure 30 shows the location of the points chosen; material point A is located at the toe of the cutting, material point B is found in the middle of the soil layer in the shear band, and material point C is in the centre of the sliding block. Figure 31 shows the p-q stress paths at the 3 points, as computed using both techniques, as well as the initial position of the yield surface for a Von Mises material (F VM ). It is seen that, for material point A, both techniques give reasonable results; this is because the bottom of the domain is fully fixed, so that the material point does not move much throughout the analysis. For material points B and  (Fig. 31b, c), the oscillations are extreme. It is evident that were a constitutive model different from Von Mises to be used, in which plasticity does not depend only on the deviatoric stress, regular MPM would not perform well. On the other hand, using DM-GC, the stress path appears to be well-behaved (Fig. 31e, f), with only some small oscillations. Based on the results obtained with the benchmarks, Table 2 summarises the advantages and disadvantages of each of the methods studied in this paper.

Conclusion
MPM is a technique that is able to handle problems involving large deformations, since material properties and the body geometry are no longer attached to a mesh. Unfortunately, the use of regular bi-linear finite element shape functions causes significant oscillations when integrating internal forces and stiffness, decreasing the accuracy of the simulations. Moreover, the grid crossing of a material point between elements and poor stress recovery create additional oscillations. A series of improvements, both novel and building upon the work of others, have been studied and combined to obtain an almost oscillation free version of MPM. It has been shown that GIMP reduces the errors caused by grid crossing, but integration using SF gradients, shown via an example using the stiffness matrix, is inaccurate due to the use of SF gradients that drop to zero at the inter-element boundaries. Using GIMP together with a double mapping integration procedure significantly reduces the stiffness matrix oscillation. Also, it has been proven that CMPM increases the accuracy of the stresses computed at the material points compared to typical MPM and GIMP. These techniques combined (termed DM-GC) increases considerably the accuracy of the MPM simulations. Moreover, since it has been observed that DM performs well when using typical finite element shape functions, and better still when using GIMP shape functions, the combination of DM with other C1-continuous methods, such as CPDI, B-spline MPM or DDMP, is a possibility which can be studied in the future. The DM and DM-G methods have the benefit of being  where ax refers to the axisymmetric strain-displacement matrix, D ax refers to the stress-strain matrix, and r p is the radial distance between the sampling point and the axisymmetric axis.
reduces to zero within a specified tolerance), the following is true: Considering an isolated system, i.e. where momentum would not be altered by external forces, it can be stated that At the beginning of the time step for an isolated system, there is no net rate of change of momentum: where m ij is equivalent to M t . Moreover, acknowledging that Finally, substituting Eq. 53 into Eq. 46 leads to the conservation of momentum for the isolated system as Note that the previous elaboration also holds for DM, DM-G and DM-GC, because (1) the modified stiffness matrix (K t ) does not affect the conservation of momentum, and (2) the GIMP and CMPM SFs satisfy the partition of unity and Eq. 52. points with (part of) their support domain in the element (smp), for DM and DM-G respectively, (ii) the nodal shape function i represents N i or S ip* , for DM and DM-G respectively, and (iii) the material point weight W represents W* or W, for DM and DM-G respectively. Summarizing, DM uses cmp, N i and W*, whereas DM-G uses smp, S ip* and W. In Fig. 33, the stress recovery steps are shown considering CMPM. It is highlighted that all matrices are related to the patch used, so are larger than in regular MPM. For example, for 4 noded elements the patch is 4 × 4 nodes, instead of 2 × 2.
To demonstrate the computational performance of the DM-G algorithm (which has a higher computational cost when compared to DM), a series of tests were conducted on square meshes, with the results presented in Fig. 34. The tests consisted of computing the stiffnesses of meshes made up of 50, 75, 100, 150, and 200 elements per side, with each element filled with 4 equally distributed MPs. Each stiffness cycle was computed 200 times to obtain an accurate mean value of the time taken. In Fig. 34, the relationship between computational time using DM-G and regular MPM is plotted as a function of the number of equations. It is observed that the relationship is almost constant, with DM-G taking about 50% longer than regular MPM. However, the overall increase of computational time for the problems studied is close to 5%, although this is dependent on the solver and characteristics of the problem solved.