Development of an implicit contact technique for the material point method

An implicit contact algorithm for the material point method (MPM) has been developed to simulate contact. This allows recently developed implicit MPM codes to simulate large-scale deformations and interaction with external bodies. The performance of the method has been investigated and compared to an existing explicit method using benchmark and geotechnical examples. In particular, the proposed formulation has been shown to conserve energy in a similar way to the explicit formulation and reach similar results. The method typically converges to the analytical solution when an adequate time step and mesh size are used, with the time step generally being around ten times larger than the explicit method, although during the contact phase this does not always result in faster computation due to the iterative solution procedure.


Introduction
The material point method (MPM) is an improvement over the classical finite element method (FEM) for problems involving large deformations. In this technique, the material is no longer attached to elements, but to particles or points (i.e. material points) that move freely through a background mesh. MPM has proven to be reliable for simulating mechanical behaviour (Abe et al. 2014, Beuth et al. 2008, Coombs et al. 2018, González Acosta et al. 2020, Jassim et al. 2013, Ma et al. 2013, Nairn 2003, Wang et al. 2016b, York et al. 1999, as well as large deformations (Coetzee 2014, Huang et al. 2015, Llano-Serna et al. 2016, Wang et al. 2016c, Więckowski 2004, Woo & Salgado 2018. As the fundamental attributes of the method involve material moving substantial distances, for example, as in landslides with an extensive runout (e.g. Wang et al. 2016a,c), the impact of this displaced material on other structures is a logical aspect to be investigated and has great value in many branches of engineering. In a seminal MPM paper, Sulsky et al. (1994) showed that different bodies of material could be simulated together on the same background grid, with a single value velocity field, allowing contact to be automatically calculated. However, the simulated contact was non-slip, (background) mesh dependent and, due to the lack of a specific contact formulation, inaccurate when considering bodies with different properties or initial state conditions.
The first dedicated methodology in MPM to simulate interaction between different bodies, i.e. contact, was proposed by Bardenhagen et al. (2000). In this method, the nodal velocities of distinct bodies are compared with the combined nodal velocities (i.e. nodal velocities which have changed due to the influence of multiple bodies), and contact is defined when these two velocities are not equal. When this happens, constraints on the normal and tangential kinematics are applied, such that the bodies cannot penetrate each other and that the tangential behaviour is defined based on constitutive behaviour, e.g. Coulomb friction. In this way, Bardenhagen et al. (2000) were able to simulate frictional slip behaviour within the MPM framework. Since then, many case studies and improvements have been presented (Bardenhagen et al. 2001, González Acosta et al. 2018, Homel & Herbold 2017, Müller & Vargas 2019, Nairn 2013, Pantev 2016, Zhang et al. 2006). Most of this work has been implemented using an explicit scheme, in line with the method of solving the equations of motion in the respective publications, and this has advantages with respect to the simplicity of the implementation. However, versions of MPM which solve the equations of motion implicitly (e.g. Guilkey & Weiss 2003, Wang et al. 2016b) have been recently developed, which include the advantages typically associated with implicit solutions; in particular, a larger time step. Therefore, the development of an implicit contact scheme is desirable.
Some work regarding implicit contact in MPM has been reported, where, in line with contact methods used in implicit FEM, Lagrange Multipliers (Chen et al. 2017) or the Penalty Method (Liu & Sun 2020) were used. The procedures followed in these techniques to simulate contact are different from the traditional method used in the explicit MPM. In this paper, an alternative form of the implicit contact algorithm is presented. This algorithm adopts the classical equations used in the explicit approach (Bardenhagen et al. 2000) for the implicit scheme and considers the Newton-Raphson iterative procedure. Moreover, to demonstrate the accuracy and performance of the proposed solution, three benchmarks and two geotechnical problems are analysed. Finally, comparative remarks are presented relating to the computational performance, followed by some recommendations for increasing the accuracy of the simulations.

MPM Formulation
The formulation of the explicit and implicit MPM equations is briefly presented and, since MPM shares similarities with the FEM continuum mechanics framework, most of the following equations are the same as in FEM. It should be noted that this formulation refers to a single body. For a detailed elaboration of these equations, the reader is directed to Bathe (1996), Belytschko et al. (2013), Sulsky et al. (1995), and Wang et al. (2016b).

Explicit MPM
The respective equations of mass and momentum conservation are dρ dt (1) ρa = ∇⋅σ + ρb (2) where ρ is the mass density, v is the velocity vector, t is time, a is the acceleration vector, σ is the Cauchy stress tensor, and b is the body force.
Besides these conservation laws, a stress-strain relationship is also required. The strain increment computed from the velocities, and the resulting stress increment, are respectively given by where dε is the strain increment vector, and D is the material elastic matrix. Using the principle of virtual power, the weak form of eq. (2), considering traction, is found by taking the product of a test function with the equation of momentum conservation and integrating over the current configuration to give ∫ where δv is the test function, s s is the traction load at the body surface Γ (i.e. the boundary condition), and V is the current body volume. Then, since in MPM two discretisations are used (i.e. material points representing the computational body and elements representing the computational space), the integrals in eq. (5) can be expressed as nodal equations via the summation of material point values in the elements. This leads to the (element) equation where N is the element matrix of shape functions (SFs) (see Appendix A for the definition of the matrices used in the paper), B is the strain-displacement matrix, a is the vector of nodal accelerations, nmp is the number of material points in the element, bmp is the number of boundary material points having a traction load, x p is the position of each material point, σ p is the material point stress tensor, ρ p is the material point density, J is the Jacobian matrix, and W p is the material point integration weight. It should be noted that the shape functions are used to interpolate material point variables to the element nodes (i.e. Nδv(x p ) ≈ δv(x p ), where δv is the vector of nodal test function velocities), and that the volume of each material point is represented by the Jacobian determinant and the material point weight (i.e. V p = |J|W p , where V p is the material point volume). Finally, since the element test function (δv) in eq. (6) should hold for any arbitrary virtual velocity that satisfies all the boundary conditions, it vanishes from eq. (6) since δvA = δvC, in which A and C represent the left and right sides of eq. (6), respectively. Hence, eq. (6) reduces further and may be expressed for an element (| elem ) using matrix notation as where In eqs. (7) and (8), m represents the element consistent mass matrix, the use of which can lead to an increase in the solution time. Hence, the element lumped mass matrix is used here, and is computed using the "row-summation" technique as Furthermore, the lumped mass matrix can be replaced by an element mass vector (m i = m ii ). This replacement improves the computational time, since matrix operations (i.e. matrix inverse computations) are replaced by linear operations. Finally, note that the use of m is particular to the explicit scheme; in the implicit scheme m will be used.

Nodal mapping, solution and update phase (explicit)
At the beginning of each solution step, the state variables are mapped to the nodes. For example, the velocity at node i is computed as where N i is the nodal SF, and the summation is over all material points surrounding the node. After nodal mapping and integration to form the element equations, an assembly operation over the mesh is undertaken (as in FEM), yielding global vectors and matrices at time t, i.e. the global mass matrix m t , the global external force F ext,t and the global internal force F int,t , which can be solved to find the global acceleration vector. Since the vector of nodal masses is used, the solution of the equation of equilibrium (in terms of accelerations) is and, using the Euler time integration method, the material point variables are updated at time t + Δt as v t+Δt x t+Δt where σ p t is the previous material point stress. The vector of (element) nodal velocities at time t + Δt in eq. (16) is computed from the material point velocities at time t + Δt as where m p is the material point mass. Note that this vector of nodal velocities is computed from the material point velocities at time t + Δt rather than using the Euler time integration method (i.e. v t+Δt = v t + a t Δt) to ensure conservation of momentum.

Implicit MPM
To derive the equation for implicit equilibrium, the principle of virtual work is used. This states that the difference between the internal and external work should be equal to zero ( i.e. Π = W int − W ext = 0 ) , and is written as where u represents the continuous displacement field. In eq. (18), the inertia term (i.e. a) is included as part of the body forces using D'Alambert's principle and, if this term is excluded, a quasi-static scheme is obtained. Next, by evaluating δΠ = 0 with respect to the displacements and knowing that δε 2 = 2εδε (Davies 2011;Bathe 1996) where δu and δε represent the test function in terms of displacements and strains. Then, following similar procedures as in eq. (6) and splitting the total stresses into the sum of the previous stress and the incremental stress (i.e. σ = σ 0 + Δσ), eq. (19) is expressed as where B is used to interpolate the virtual strains using virtual nodal displacements, i.e. B(x p )δu ≈ δε(x p ), where δu are the nodal virtual displacements, σ 0 represents the vector of initial stresses, and u represents the vector of (element) nodal displacements which is obtained from the interpolation of strains, i.e. Bu = ε. As for eq. (7), the previous equations represent equilibrium and can be expressed in matrix form as where the element stiffness matrix is and F ext | elem and F int | elem are similar to eq. (9) and eq. (10), respectively. If the inertia term is ignored, the equation of equilibrium reduces to which represents quasi-static equilibrium. Again, an assembly operation over the mesh is undertaken (as in FEM), yielding global vectors and matrices. By using Newmark's (1959) time integration technique, which is recognised to be highly effective (Bathe 2007), eq. (21) can be solved. The velocities and displacements at time t + Δt are computed as where α and γ are time stepping parameters that are here chosen to be α = 0.25 and γ = 0.5 (giving a constant-average-acceleration approach). Then, isolating a t+Δt from eq. (25) leads to and substituting eq. (26) into eq. (24) results in Finally, by substituting eq. (26) into eq. (21) and adding the Newton-Raphson iteration procedure, the equation of equilibrium is written as or where Δu = u t+Δt -u t is the vector of incremental displacements, and In these equations, K represents a modified stiffness matrix, F kin are the kinetic forces, and the left superscript k refers to the Newton-Raphson iteration step. The iteration is stopped after the desired convergence criterion has been reached, which in this paper is computed as (Guilkey & Weiss 2003) where tol is the tolerance value. There are several alternative convergence criteria which could be used. However, this convergence criterion has been selected due to its simplicity, although it is noted to be susceptible to false convergence due to iteration stagnation.

Nodal integration and update phase (implicit)
The integration of nodal external loads, mass, velocity and stiffness has already been given in eqs. (9), (11), (12), and (22), respectively. In contrast to the explicit scheme, nodal accelerations are also needed and are computed as After assembling the stiffness and mass matrices (eq. (30)), the equilibrium in terms of displacements is computed (eq. (29)) and the vectors of internal forces, nodal accelerations and kinetic forces are updated every iteration step using eqs. (10), (26), and (31), respectively. Finally, after the desired iteration tolerance is reached, a material point's acceleration, velocity, position and stress are updated as

MPM mapping, integration and stress recovery enhancement
Since standard MPM suffers from oscillation problems, the procedures elaborated in González Acosta et al. (2020) are used to mitigate such inaccuracies. In this paper, the double mapping DM-GC technique is used, in which the (element) vectors of external and internal nodal forces are integrated using standard GIMP SFs (Bardenhagen & Kober 2004) as where S ip is the element matrix of GIMP SFs, and ∇S ip is the strain-displacement matrix for the GIMP SFs. Since the implicit solution scheme is used, the (element) stiffness matrix is computed using double mapping (DM) procedures as where N i is the nodal SF, S ip* is the local GIMP SF of node i, ngauss is the number of Gauss integration points in the element, smp is the number of material points with a support domain inside the element, and W FE is the finite element (FE) weight associated with the Gauss point g. In González Acosta et al. (2020), the DM implementation and its performance are presented. The incremental stress of a material point is then updated for the explicit and implicit schemes using CMPM (González Acosta et al. 2017) as and where B C is the (element) matrix of the CMPM SF gradients (i.e. the strain-displacement matrix for CMPM), and u ext and v ext are the vectors of nodal displacements and velocities, respectively, in the extended CMPM domain. Note that since GIMP and local GIMP SFs are used to integrate both the forces and stiffness matrix, the stencil (i.e. support domain of each material point) is the same in both cases, including when using double mapping. CMPM uses an extended domain (u ext or v ext ) to recover stresses, which is not necessarily equal to the GIMP stencil. An explanation of the stencil is given in Appendix B, and the accuracy of using DM-GC is demonstrated in Appendix C. Finally, it is important to remark that the state variables (i.e. mass, velocity, acceleration) are mapped using GIMP SFs.

Contact formulation
The contact formulation relies on the equations of motion presented in the previous section, which are solved for a series of separately defined bodies, with the interactions (contact) between them being included in the equations of equilibrium as additional external forces. The bodies may be pre-defined (as in the examples here), or identified by the program to allow for the separation or joining of bodies. The original contact procedures in MPM (Bardenhagen et al. 2000) use two steps to simulate contact: (i) the velocity field sharing step, and (ii) the contact force evaluation step. In this work, an extra (proximity detection) step to ensure mesh independency has been included, using a distance limit to activate the contact. In the velocity sharing step, the velocities are interpolated to the nodes to detect if any of the bodies may be interacting (which is possible if two or more bodies interpolate velocities to the same node). Using the individual nodal velocities of each body and the combined domain velocities (i.e. the nodal velocities accumulated from all material points in the domain), velocity sharing is detected at node where v i,C is the nodal combined velocity accounting for all bodies in the domain, v i,bod is the nodal velocity of each independent body, m i,C is the combined nodal mass, bod denotes the body, and nb is the number of bodies in the domain. If eq. (43) is true for any node i, the proximity detection rule is then evaluated. This proximity condition is where x i p, bod1 and x i p, bod2 are the coordinates of the closest material points p from each body (bod1 and bod2) to the possible contact node i, d is the distance between x i p, bod1 and x i p, bod2 , and d min is the minimum distance required to activate the contact. Note that d min should have a maximum value of the cell size to ensure mesh independence, but it has no minimum value. It should be noted that, in eq. (46), only the contact between two bodies (bod1 and bod2) is considered, and the formulation to account for contact at a single node by more bodies is not included. If the velocity and proximity conditions are satisfied (i.e. eq. (43) and eq. (46)), it is then necessary to evaluate if the bodies are approaching each other (contact evaluation step). This last condition is relevant because, if the bodies are separating (i.e. the first two conditions are true but the bodies are moving in opposite directions), interaction is considered not to be taking place; otherwise it would add energy into the system. This step is formulated as where n i,bod is the unit normal vector of the body at the contact node. The normal vector can be computed using different approaches (Nairn 2013), but in this paper the approach described in Huang et al. (2011) is used, which is where the normal n is the individual normal of the body bod around node i, in which the direction is governed by the mass m g concentrated at the (central) Gauss position, and ngp is the number of (central) Gauss positions with a concentrated mass m g around node i. Note that the normal n can be inconsistent (i.e. may not be perpendicular) to the surface of the body due to the irregular distribution of m g around node i. This inconsistency can be lessened by using n, the normal direction at node i in which the influence of the neighbouring body is considered. Furthermore, to compute n, bi-linear or GIMP SF gradients can be used, as the results are identical due to the centrally located Gauss point. Fig. 1 illustrates the different variables during contact between body A and body B, and highlights the difference between n and n. If eqs. (43), (46), and (47) are true, contact is deemed to occur, and contact conditions must be applied by assigning contact forces. To do so, first the corrected nodal velocities of each body are calculated to avoid interpenetration. Hence, where *v ic, bod is the corrected nodal velocity, and ic represents the contact node. Then, the normal contact force is computed as where F nc ic, bod is the nodal normal contact force. The tangential force required to ensure non-slip (i.e. stick) conditions is calculated as where F stick ic, bod is the frictional force required to ensure non-slip conditions. The frictional force, obeying Coulomb's theory, is then where F fric ic,bod is the final frictional force between the bodies, μ is the friction coefficient which depends on the material characteristics, the quotient ensures that the direction of the force opposes the relative motion, and the min function determines the magnitude of the force. Finally, if contact occurs, nodal accelerations are no longer computed using eq. (13), since the contact loads should be included as

Explicit contact algorithm
The implementation of contact using the explicit scheme is straightforward. Fig. 2 shows the algorithm of the explicit scheme including the contact steps. Several minor steps are omitted, such as the activation of the boundary conditions and the computation of the material point local coordinates, in order to maintain the algorithm simplicity. The square brackets located at the left side of the figure indicate a loop for each body; i.e. the combination of steps should be repeated for each body before continuing with the following steps. The plastic iteration loop is indicated, but the procedures are not elaborated; the reader is directed to Smith et al. (2013) and Sloan et al. (2001) where the plastic procedures are described in detail.
Note that the steps shown in Fig. 2 indicate that the modified update stress last (MUSL) method (Sulsky et al. 1995) is used. The standard update stress last (USL) method consists of computing the increment of stresses after the update of nodal velocities, but, since the nodal velocities are computed from the updated material point velocities, the MUSL is obtained. The MUSL algorithm has been shown to conserve energy better than the USL approach (Nairn 2003, Pantev 2016.

Implicit contact algorithm
To simulate contact using the implicit scheme, the same approach as in the explicit scheme is used. Newmark's time integration technique and the Newton-Raphson method are used, so the equations used to estimate contact must be updated in each iteration step. Moreover, if contact occurs, it is assumed that it persists throughout the Newton-Raphson iterative process, i.e. between times t and t + Δt. Considering these assumptions, the equations used to update the nodal velocities during the iterative procedure are then and The contact loads are computed as Note that these equations are similar to those used in the explicit scheme, with the exception that the nodal combined velocities (v t+Δt ic, bod ) are updated every iteration, which leads to a new set of corrected velocities and nodal forces. Finally, if contact occurs, the incremental displacements should be computed using Fig. 3 shows the steps followed to simulate contact using the implicit scheme. It should be noted that the vector of external forces (F ext ) and the stiffness matrix (K) are evaluated at time t + Δt outside the iteration loop, as both remain constant throughout the iterative procedure, i.e. the modified Newton-Raphson method is used. The small time steps used to capture the kinematics ensure that this results in a reasonable solution.

Benchmark problems
To evaluate the performance of the implicit contact method, three benchmarks are analysed and the solution of each compared against the explicit solution, as well as against the analytical and FEM solutions where available. The first two benchmarks (i.e. a collision and a 1D vibrating bar) are used to study energy conservation. The first benchmark considers a contact occurring instantaneously between two blocks, and the second considers contact which persists throughout the whole simulation. The third benchmark is the simulation of a block sliding on a rigid surface, and is used to analyse the interaction between bodies when considering frictional forces. Each benchmark was simulated using four equally spaced material points per background element, which were  initially positioned at the local coordinates ξ = ± 0.5 and η = ± 0.5. All the benchmarks and geotechnical simulations introduced later are simulated using plane strain conditions.

Collision benchmark
In Fig. 4, the initial configuration of the collision benchmark is shown. The benchmark consists of two square blocks (A and B) of size 0.6 m, moving freely towards each other through a background mesh with a mesh spacing of Δ x = Δ y = 0.20 m. Each block has an initial horizontal velocity of 0.5 m/s. The elastic properties of both blocks are Young's modulus, E = 500 kPa, and Poisson's ratio, ν = 0.30. The tolerance values used in the Newton-Raphson iterative procedure, and to establish contact, are tol = 1.0 × 10 -8 and d min = 0.2 m, respectively. Gravity forces have not been considered. Fig. 5 shows the results of the collision in terms of total energy (TE), for both the explicit and the implicit schemes. In Fig. 5a, for the explicit scheme, it is observed that, at the moment of the collision (t c ), the energy decreases significantly (by 17%) and, depending on the time step used, the total energy either continues decreasing or stays almost constant after the collision. Fig. 5b shows that, using the implicit scheme, the total energy behaves in a similar way to the explicit solution, except that some oscillation of the energy is observed after the collision. These oscillations are attributed to not capturing accurately the internal stress waves due to the large time step used, which is 10 times greater than that used with the explicit scheme for almost the same performance.
Additionally, Fig. 6 and Table 1 demonstrate the energy conservation and computational time required to complete each simulation with respect to the mesh (element) size for the implicit scheme. The study of energy conservation is performed using a time step of Δt = 1.0 × 10 -4 s, and mesh sizes of 0.20, 0.10 and 0.05 m. Also, an additional curve (solid line with markers) computed using a mesh size of 0.05 m and a time step of Δt = 1.0 × 10 -5 s is added to Fig. 6, to study the effects of reducing the time step and the mesh size simultaneously. The contact tolerance values used are the maximum values for the mesh, i.e. d min = 0.20, 0.10 and 0.05 m, respectively. It is observed that the drop of energy at the point of contact decreases with the reduction of the mesh size, improving the conservation of energy. Furthermore, by reducing the mesh size and time step together, the drop of energy and the oscillations after contact reduce drastically, reinforcing the theory that the oscillations in Fig. 5b were caused by the large time step. However, despite the improvement in the results, some reduction in energy is still observed.
An analysis of the computational time has been performed, considering the time steps Δt = 1.0 × 10 -5 s and Δt = 1.0 × 10 -4 s, for the explicit and the implicit schemes, respectively. It is observed that the computational time needed to complete a simulation grows substantially with the reduction of the mesh size. Nevertheless, the time relationship is almost constant, with the implicit solution being on average 30% faster than the explicit solution.
The reason for the unrealistic decrease of energy observed in Figs. 5 and 6 is the inconsistency between the contact loads and the internal loads developed at the contact (F nc ∕ = F int ). One method to prevent this inconsistency is to compute the contact loads as a function of the internal loads developed by each body at the contact (Xiao-Fei et al. 2008). At the contact interface, the normal acceleration of each body can be assumed to be equal (i.e. a t i, A ⋅n i, A = − a t i, B ⋅n i, B ). Hence, from Newton's second law, it is established that which reduces to where F nc i, A and F nc i, B are the contact forces between the bodies as a function of the internal forces developed. In Fig. 7, the energy conservation of the collision benchmark is plotted using eq. (63) to compute the contact loads, a step size of Δt = 1.0 × 10 -5 s, a mesh of size of Δ x = Δ y = 0.20 m, and the explicit scheme. It is seen that the energy conservation improves considerably because the contact loads and the internal loads are equivalent. However, this approach should be implemented with caution, since the contact forces are dependent on the material point stresses which are known to suffer from numerical oscillations, leading to oscillations of the contact loads. Since stress oscillation problems are typical in MPM, in order to implement eq. (63) in any problem more complicated than this example, further improvements are needed. Additionally, the implementation of this contact technique using the implicit scheme is not straightforward, due to the interdependence of the internal and contact loads and the iterative nature of the implicit scheme. Since the improvements needed to implement this technique in both solution schemes are beyond the scope of this paper, all remaining examples have been computed using the initial approach.

1D bar benchmark
This benchmark is similar to other 1D vibration benchmarks, except that this 1D bar is made up of two bars (bodies A and B). Initially, the bodies have zero internal stresses, but, due to the influence of gravity, they compress as soon as they are allowed to move. The equations of motion are calculated separately for the two bodies and the interaction between them is calculated via the contact forces. Figs. 9 and 10 show the sum of the total energy of each 1D bar and the displacement of the top material point over 3 s using various time steps, for the explicit and implicit solution schemes, respectively. An FEM result is also included to provide a reference solution. As can be seen in Fig. 9a, by reducing the time step the energy conservation is improved. The largest drop of energy is during compression, when most of the penetration occurs. Using the smallest time step of Δt = 1.0 × 10 -5 s, the loss of energy is the lowest (4.7% after 3 s). On the other hand, using the largest time step, Δt = 1.0 × 10 -4 s, the loss of energy is considerable (22.5% after 3 s). Consistent with the calculated loss of energy, Fig. 9b shows that the (a) (b) Fig. 9. a) Sum of the total energy of each 1D bar using the explicit scheme, and b) displacement of the top material point. top material point is unable to reach its original position when a large time step is used. In contrast, when the smallest time step is used, the result is nearly equal to the FEM result (i.e. accurate). This is because the computed contact loads can adjust satisfactorily to the change of internal forces and velocities of each body during contact when smaller time steps are used, restricting the penetration of the bodies.
When the implicit scheme is used (Fig. 10), the results for energy conservation and material point position are similar to those using the explicit scheme, based on each implicit simulation using a time step 10 times larger than in the explicit simulation.

Sliding benchmark
In this benchmark, a block under gravity loading is located on a flat surface and slides. Fig. 11 is a sketch of the problem, showing that part of the background mesh where each element was initially filled with 4 material points. Body A is a rigid surface and a body B slides over body A due to the gravity load (g) which is rotated by θ = 15 • to the vertical. The lengths of the two bodies are L 1 = 1.5 m and L 2 = 6.0 m, respectively, and the height of each body is H 1 = H 2 = 0.75 m. The mesh size is Δ x = Δ y = 0.25 m, and body A is fixed at the bottom boundary and at the vertical boundaries. The elastic properties of the bodies are Young's modulus, E A = 1000 kPa and E B = 5000 kPa, in which E A and E B refer to body A and B, respectively. The Poisson's ratio and unit weight for both bodies are ν = 0.35 and γ = 20 kN/m 3 , respectively. Two simulations using two different friction factors, μ = 0 and μ = 0.15, have been performed, and the results are compared with the analytical solution of Jiang & Yeung (2004). Time steps of Δt = 5.0 × 10 -5 s and Δt = 5.0 × 10 -4 s were used for the explicit and the implicit simulations, respectively. The tolerance values used in the Newton-Raphson iterative procedure and the contact distance are tol = 1.0 × 10 -8 and d min = 0.25 m, respectively. Fig. 12 shows the results of the benchmark over two seconds of sliding. It is seen that the simulated results using the explicit and the implicit schemes are close to the analytical solution. Nevertheless, small deviations can be seen, which are attributed to the non-smooth distribution of the contact loads at the interface between bodies A and B. Fig. 13 shows this distribution at the beginning of the simulation (i.e. no slip had yet occurred) and after some sliding considering μ = 0.15. Fig. 13a shows that the normal load F nc is constant at the centre of the block and drops close to the corners. In addition, it shows that the frictional load F fric is developing from right (F fric = -0.43 kN at 1 m) to left (F fric = 0.29 kN at 0.25 m) due to the weight of body B and the inclined gravity force. The frictional force has both positive and negative components. After slip occurs, Fig. 13b shows that the frictional loads are constant below body B, indicating that the frictional force is fully developed. Additionally, it is observed that the results obtained using the explicit and the implicit schemes are similar. The reader is directed to the work of Oden & Pires (1984) in which a similar simulation was performed, although under static conditions with loading applied via a traction on a single edge. The similarities in the (a) (b) Fig. 10. a) Sum of the total energy of each 1D bar using the implicit scheme, and b) displacement of the top material point.  qualitative distributions of the normal and frictional loads before the movement of the block give confidence in the analysis. Figs. 14 and 15 show the average of the normal and tangent loads, respectively, below body B during the sliding. It is observed that both loads are close and oscillate around the analytical solution (which is computed considering static conditions). These oscillations are caused mainly by the vibration of body B during sliding. This is more evident in the explicit results, which exhibit larger oscillations.

Geotechnical applications
Two geotechnical problems are simulated using both the explicit and the implicit contact schemes. The first problem consists of a shallow foundation that penetrates into the soil due to an increase in gravitational loading applied to the foundation. The second problem consists of a vertical cutting that, after failing, slides and collides against a wall. In these problems, the Young's modulus used for the rigid elements (i.e. the foundation and the wall) are lower than realistic values for such elements. The reason for this is to permit a limited amount of bending of the 'rigid' bodies, to demonstrate the capabilities of the method.

Foundation
This problem involves a foundation slab that penetrates into an incompressible, linear elastic, perfectly plastic von Mises soil, with a rough interface, i.e. no tangential slip is allowed. The initial stresses in the soil are computed as σ y = γ soil h y and σ x = σ z = K 0 σ y , where γ soil = 17.5 kN/m 3 is the soil unit weight, h y is the depth of the material point below the soil surface, and K 0 = 0.7 is the coefficient of earth pressure at rest. The elastic parameters of the soil are Young's modulus, E soil = 5.0×10 3 kPa, and Poisson's ratio, ν soil = 0.49, and its undrained shear strength is S u = 10 kPa. The slab Young's modulus and Poisson's ratio are E slab = 1.0 ×10 4 kPa and ν slab = 0.45, respectively. Based on the results of the benchmark analyses, the time steps selected for the explicit and implicit simulations were Δt expl = 1.0 × 10 -5 s and Δt impl = 1.0 × 10 -4 s, respectively. The pressure exerted by the slab on the soil has been applied by increasing the self-weight of the slab by gravitational loading; i.e. by assuming the density of the slab to be ρ slab = 3000 kg/m 3 and increasing the gravity from zero at a rate of 10 g/s. Finally, the tolerance value used in the Newton-Raphson iterative procedure and the contact distance are tol = 1.0 × 10 -8 and d min = 0.20 m, respectively. Fig. 17 shows the results for the foundation problem after a settlement of 0.20 m using the implicit scheme (with similar results also being obtained with the explicit scheme). In Fig. 17a, the plastic deviatoric strain contours indicate a failure mechanism resembling the failure described in Fig. 16. Foundation problem.  Fig. 15. Interface tangent force (F fric ).   Terzaghi (1943), while Fig. 17b shows the shear stresses developed. By comparing both parts of the figures, it is seen that the larger plastic deviatoric strains develop in zones where the shear stresses shift from positive values to negative values. It is also evident that, due to the soil yielding and the use of a large Poisson's ratio, some checker-board oscillation occurs, as is typical for von Mises incompressible materials using bi-linear elements. The reader is referred to Coombs et al. (2018) and González Acosta et al. (2019), where the volumetric locking problem has been studied. Fig. 18 shows the pressure (P L )-displacement responses computed using the explicit and implicit schemes, plotted against two analytical solutions of Terzaghi (1943). Both of these solutions for the bearing capacity consider rough contact between the soil and foundation. However, one solution does not include any overburden load (Q 1 = S u N c ), whereas the second solution does (Q 2 = S u N c + γ soil D f , where D f is the depth of the base of the foundation). The pressure load between the slab and the soil is computed from the contact forces (F nc ) and the size of the foundation. As can be seen, despite the stress oscillations observed in Fig. 17, the maximum simulated bearing capacity of the soil is close to both analytical solutions, especially the second solution in which the influence of the deformation has been included (via the overburden component).

Vertical cutting
This problem consists of a vertical cutting that fails due to self-weight and then hits a protection wall. Fig. 19 shows a schematic of the initial state of the problem, including the boundary conditions. As for the previous example, the initial stress states in the vertical cutting and the protection wall are computed using the depth of each material point from the surface of each structure (i.e. h y1 and h y2 , respectively), and the unit weights and earth pressure coefficients of the soil and the wall, γ soil = 18.0 kN/m 3 , K 0-soil = 0.7 and γ wall = 25.0 kN/m 3 , K 0-wall = 0.5, respectively. To generate the initial stresses in the vertical cut, a static step is applied in which the movement of the material points is prevented. This step causes large localised shear stresses to develop at the toe of the vertical cut, which cause the triggering of the landslide when the material points are released. The boundary conditions are that both bodies are fully fixed at the base, and the soil is also fixed in the horizontal direction at the lefthand vertical boundary. The height and width of the modelled soil domain and the protection wall are H 1 = 2.5 m and L 1 = 6.0 m, and H 2 = 1.0 m and L 2 = 0.8 m, respectively. The distance between the soil and the wall is L 3 = 1.0 m, and the element size is Δ x = Δ y = 0.10 m.
The elastic parameters of the soil and the wall are Young's modulus, E soil = 1.0 × 10 3 kPa and E wall = 3.0 × 10 3 kPa, and Poisson's ratio, ν soil = 0.38 and ν wall = 0.30, respectively. The soil is modelled as a linear elastic, linear strain-softening von Mises material, with a peak cohesion, c p = 14 kPa, a residual cohesion, c r = 5 kPa, and a softening modulus, H s = -18 kPa (see Wang et al. 2016b for details of the constitutive model). The wall is simulated as a linear elastic material. The time steps selected for the explicit and the implicit simulations are Δt expl = 5.0 × 10 -5 s and Δt impl = 5.0 × 10 -4 s, respectively. The tolerance value used in the Newton-Raphson iterative procedure and the contact distance are tol = 1.0 × 10 -8 and d min = 0.10 m, respectively. Fig. 20 shows the distribution of plastic deviatoric strains in the vertical cut, and the distribution of deviatoric stresses in the wall, after a) 0.8, b) 1.25 and c) 2.4 s in the implicit simulation. Fig. 20a shows a moment after the vertical cut has failed, with the soil mass moving downwards and to the right, and increasing its kinetic energy before the Fig. 19. Vertical cutting and protection wall (mesh indicative and not to scale). collision. Fig. 20b shows the instant at which the contact load and the internal deviatoric stresses in the wall are a maximum. Fig. 20c shows the end of the simulation, when the wall has partially recovered its position. It can be seen that, due to the adopted proximity detection rule, the two bodies are closely connected during contact. Fig. 21 shows the average pressure at the wall surface and the deviatoric stress at the wall base. The wall pressure in Fig. 21a has been computed by adding the contact force from each contact node and dividing it by the wall height H 2 . The deviatoric stress q in Fig. 21b has been computed, using the principal stresses, as The material points selected to compute the deviatoric stress are those inside four elements at the bottom left corner of the wall (Fig. 21b). These points were chosen because, in this section of the wall, the deviatoric stress is a maximum. It is observed that the results in Fig. 21 are consistent with each other, in that the maximum deviatoric stress occurs at the same time as the maximum contact pressure. Also, due to the elastic behaviour of the wall, an oscillatory behaviour is observed.

Computational time
The time steps used and the computational times recorded for the simulations from Sections 5.1 and 5.2 are presented, to further assess the capability of the implicit solution. Note that the time steps chosen for the geotechnical simulations were based on the results obtained from the benchmarks (i.e. the largest time steps which returned accurate results). The codes were compiled using Intel Parallel Studio XE, and run on a computer with an Intel Xeon E5-1620 processor and 16 GB RAM. However, it should be noted that no special effort was made to optimise the code (e.g. no parallel computing), which means that the installed RAM was not exploited to its maximum capacity. Table 2 summarises the implicit/explicit time step relationship, computation time (CT), and memory usage (MU) between the simulations. It is observed that, in most of the simulations (i.e. except the block collision simulation), the time step used for the implicit simulations is 10 times larger than the time step used for the explicit simulations. Due to the use of a larger time step, the computational time is smaller (in most of the simulations) when using the implicit solution. The reason for the larger computational time in the foundation simulation is the incompressibility of the soil, which leads to more iteration steps in order to reach equilibrium. However, in more realistic simulations, where contact only occurs for part of the simulation time, the implicit algorithm should lead to substantial computational time savings. Finally, it is observed that the memory usage is almost the same in both cases. This is because the implicit scheme stores the global matrix in the form of a vector using the skyline strategy, and the adopted solver preserves the skyline during the solution step. Nevertheless, when the problems simulated contain a large number of points (e.g. the foundation and the soil/wall collision problems), the differences in memory usage are more evident.

Conclusion
MPM is able to simulate large-strain mechanical behaviour and, with the development of the contact solution, the interaction between multiple bodies is feasible. In this paper, an implicit contact solution was developed for MPM, based on existing explicit procedures and using the Newton-Raphson iterative procedure and Newmark's time integration scheme. The relative performance of the formulation with respect to an explicit scheme was investigated through analysing three benchmark problems. It was shown that the results obtained using the explicit and implicit methods were almost identical. It was observed that both methods dissipate energy in a similar manner due to an inconsistency between the contact loads and the increment of material point stresses, which is proportional to both the time step and mesh size discretisation. Moreover, it was observed that by reducing the time step and mesh size to realistically small sizes, a reasonable energy conservation can be achieved. The ability of the method for simulating geotechnical problems was demonstrated through the analysis of foundation slab settlement in a cohesive soil, in which the failure mechanism and the bearing capacity obtained were similar to Terzaghi's analytical solution, and through the analysis of the collision of a mass of soil against a protection wall. Finally, it was shown that the simulation using the implicit approach is faster when the contact loads in the problem are moderate; otherwise, the computational time can be similar to the explicit approach.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  where H represents the matrix of FE bi-linear SFs (N) or GIMP SFs (S ip ), B represents the strain-displacement matrix B, ∇S ip or B C , the subscript indicates the node number, and n = 4 when bi-linear SFs are used or n is the number of nodes inside the material point support domain when GIMP is used. On the other hand, the relationship between stresses and strains (σ = Dε) considering a plane strain condition is where the vector of strains (ε = Bu) is where n is the number of nodes in the extended domain as shown in Fig. A1. Note that the extended domain can be constructed using square arrays of C 2 elements, where C > 1.

Appendix B
In this appendix, the combination of using GIMP SFs and the CMPM extended domain is explained. Fig. B1 shows a domain composed of 9 elements and 16 nodes (both numbered in the y direction, from left to right). The domain is filled with several material points, and a particular material point (indicated by the larger red cross) is used to demonstrate the GIMP and CMPM stencils. Since GIMP SFs and local GIMP SFs are used to integrate the nodal forces and stiffness, the material point support domain is used to select the elements. In this case, four elements are selected (4,5,7,8). On the other hand, CMPM has a different "extended" domain which, in this case, uses the full domain to enhance the accuracy of the stresses computed using Lagrange interpolation functions. It should be noted that the use of a different domain to perform nodal and stiffness integration and to recover stresses does not represent a problem, since the equilibrium (Eqs. (13), (21), (23), (28), (54) and (60)) is computed using only GIMP and, after equilibrium, the stresses are computed using CMPM. The reader is directed to González Acosta et al. (2020) for further information on the implementation of DM-GC, where the case of boundary material points is presented. Fig. A1. Extended domain for a material point using an array of 3 2 elements.

Appendix C
To demonstrate the accuracy of DM-GC, a 1D benchmark problem is analysed. This problem is similar to the benchmark analysed in Coombs et al. (2018), in which a 1D column of material points is subjected to an increasing gravity force, and the error as a function of the number of square background elements (nel) is investigated. In this benchmark, the initial gravity load is zero, and the gravity is increased in increments of Δg = 1.0 × 10 -4 g until maximum gravity loads of g SD = 0.1g and g LD = 20g, in which the first maximum gravity load (i.e. g SD ) is used to study small deformations, and the second (i.e. g LD ) is used to study large deformations. The 1D bar has an initial height and width of H 0 = 10 m and L = 1 m. The elastic parameters of the bar are Young's modulus, E = 1.0 × 10 3 kPa, and Poisson's ratio, ν = 0.0. The density of the material points is ρ = 1500 kg/m 3 . The bar is fully fixed at the bottom, and fixed in the horizontal direction along both vertical sides. Each background element is filled (initially) with four equally-spaced material points, and the simulations have been performed using the static scheme to avoid kinematic effects. The error during the calculations is measured as   in which σ p,z is the material point vertical stress, σ a,z (z) is the analytical vertical stress in the bar at the height z (i.e. σ a,z = ρg(H-z), in which z is measured from the bottom), ‖•‖ is the L 2 norm of (•), V 0 p is the initial material point volume, V 0 BAR is the initial bar volume, and ρgH 0 is the vertical stress at the bottom of the bar used to normalise the difference in the stresses. Fig. C1 shows the convergence of MPM, GIMP and DM-CG considering small deformations. It is observed that MPM and GIMP have the same convergence (α ≈ 1, in which α is the order of convergence). Nevertheless, it is observed that the MPM error cannot reduce linearly at a certain point, which is attributed to the element crossing problem. On the other hand, using DM-GC the order of convergence is 1 < α < 2. Fig. C2 shows the results considering large deformations. In this case, it is observed that the MPM error grows due to the large number of material points crossing elements, which causes highly inaccurate internal stresses. On the other hand, the GIMP and DM-GC errors reduce at a similar rate (with DM-GC having a smaller error).