Imposition of essential boundary conditions in the material point method

There is increasing interest in the material point method (MPM) as a means of modelling solid mechanics problems in which very large deformations occur, e.g. in the study of landslides and metal forming; however, some aspects vital to wider use of the method have to date been ignored, in particular methods for imposing essential boundary conditions in the case where the problem domain boundary does not coincide with the background grid element edges. In this paper, we develop a simple procedure originally devised for standard finite elements for the imposition of essential boundary conditions, for the MPM, expanding its capabilities to model boundaries of any inclination. To the authors' knowledge, this is the first time that a method has been proposed that allows arbitrary Dirichlet boundary conditions (zero and nonzero values at any inclination) to be imposed in the MPM. The method presented in this paper is different from other MPM boundary approximation approaches, in that (1) the boundaries are independent of the background mesh, (2) artificially stiff regions of material points are avoided, and (3) the method does not rely on mirroring of the problem domain to impose symmetry. The main contribution of this work is equally applicable to standard finite elements and the MPM.


Summary
There is increasing interest in the material point method (MPM) as a means of modelling solid mechanics problems in which very large deformations occur, e.g.
in the study of landslides and metal forming; however, some aspects vital to wider use of the method have to date been ignored, in particular methods for imposing essential boundary conditions in the case where the problem domain boundary does not coincide with the background grid element edges. In this paper, we develop a simple procedure originally devised for standard finite elements for the imposition of essential boundary conditions, for the MPM, expanding its capabilities to model boundaries of any inclination. To the authors' knowledge, this is the first time that a method has been proposed that allows arbitrary Dirichlet boundary conditions (zero and nonzero values at any inclination) to be imposed in the MPM. The method presented in this paper is different from other MPM boundary approximation approaches, in that (1) the boundaries are independent of the background mesh, (2) artificially stiff regions of material points are avoided, and (3) the method does not rely on mirroring of the problem domain to impose symmetry. The main contribution of this work is equally applicable to standard finite elements and the MPM.

The MPM
The novelty of the MPM is its combination of Lagrangian and Eulerian approaches into a single procedure. Since its first presentation in the 1990s 7 many papers have appeared, both developing the MPM and applying it to various problems, as indicated above. In most papers, the MPM is presented in explicit form, however, there are advantages in using an implicit approach 8 ; for example, implicit approaches allow for much larger loadsteps to be applied, have improved stability, and contain error control within each loadstep. In particular, Guilkey and Weiss 8 demonstrated the advantage of adopting an implicit approach for quasi-static analysis in terms of accuracy, and we follow an implicit approach in this paper. The focus of this paper is on incorporating (and extending) an implicit boundary (IB) method within the MPM; we therefore restrict ourselves to infinitesimal strains and isotropic linear elastic material behaviour. It should however be noted that the IB method presented in this paper does not change if we remove these restrictions. The majority of this section uses index notation to detail the formulation; only the discrete nodal equations are given in matrix-vector form for convenience.
The weak form of equilibrium of a body subject to an admissible displacement field, u i , when inertia is neglected, can be expressed as where the domain, Ω, is subject to body forces, b i , and tractions, t i , on the boundary of the material domain, Ω, causing a set of admissible virtual displacements through the body, i . Within this statement of equilibrium, the conventional infinitesimal strain is the fundamental variable that characterises the deformation at an MP where x i are the coordinates of the MP and u i are the MP's displacements. In this paper, a linear relationship is assumed between engineering strain, ij , and the Cauchy stress, ij , that is, where D e ijkl is the conventional linear elastic isotropic material stiffness tensor.

Discrete implementation
We start by introducing the element approximation for the displacements at a point within an FE where {d e } and {d e } are the physical and virtual element nodal displacements, respectively, and [N] is the shape function matrix. The Galerkin form of the weak statement of equilibrium over an element, E, can be obtained from Equations 1 and 4 as where is the strain-displacement matrix containing the spatial derivatives of the shape functions with respect to the nodal coordinates, [L] is the gradient operator matrix, and the superscript (·) R denotes a residual. The first term in Equation 5 is the internal force within an element, and the combination of the second and third terms is the external force vector. Note that the integration of the traction vector remains over the boundary of the domain as we do not assume that the boundaries of the physical domain and the element edges coincide. The Galerkin form of the weak statement of equilibrium, Equation 5, can be discretised by MPs into the following form where the sub or superscript p denotes a quantity associated with an MP, n mp is the number of MPs within the element under consideration, and v p is the MP volume. In Equation 6, the MP volume, v p , replaces the quadrature weights and Jacobian volume mapping seen in conventional FE methods. The volume of an MP can be specified in a number of ways. If the MPs are equally spaced across the physical domain, the overall volume can be split equally between the MPs. An alternative is to specify the MPs' initial positions within each element based on a Gauss quadrature scheme. The volumes of these points would then be the appropriate Gauss quadrature weights multiplied by the determinant of the Jacobian of the parent element (in the case of a uniform grid, this would be a constant factor for all elements). Details of the approach for specifying MP volumes in this work are given in Section 4. In Equation 6, the shape function, [N], and strain-displacement, [B], matrices are determined by evaluating the conventional FE shape functions (or their derivatives) at the MP position within the element. Linearising Equations 5 or 6 with respect to the unknown nodal displacements and assuming that the applied body forces and surface tractions are independent of the nodal displacements gives the element contribution to the global stiffness matrix as In MP methods, Equation 7 is evaluated through the summation of the MP stiffness contributions within each element, where the nodal stiffness components of a single MP are given by and an MP's contribution to the internal force vector can be determined from The current stress and strain states at the MP are obtained through where {Δd e } are the incremental element nodal displacements associated with the nth loadstep, { p n−1 } is the MP stress from the previous loadstep (or the original state on the first loadstep), and {Δ } is the strain increment at the MP associated with the loadstep.
The MPM presented in this paper is an incrementally linear method. Within each loadstep, the material behaviour is assumed to be linear and the MP contributions to the internal force vector, Equation 9, are linearly dependent on the nodal displacements. However, the MP positions are updated at the end of each loadstep through where {x p n−1 } is the MP position from the previous loadstep and [N] contains the shape functions associated with the MP's parent element. This means that the stiffness contribution of an MP to its parent element will change between each loadstep. Therefore, the global stiffness matrix and the internal force vector must be recalculated at the start of each loadstep.

Algorithmic steps
In this paper, we adopt an implicit quasi-static formulation of the MPM and the steps in the implemented algorithm are concisely summarised below.
The applied body forces and/or tractions are split into a number of loadsteps, and for each of these steps, the following process is adopted: where {Δd} are the nodal displacements within the loadsteps and {f R } is the out of balance force residual between the stresses at the MPs and the external tractions/body forces. The MP stresses, { p }, used in Equation 6 are obtained from the end of the previous loadstep; 4. the MP positions and stresses can then be updated through interpolation from the incremental nodal displacements using Equations 11 and 10, respectively; 5. reset or replace the background grid.
Note that in this paper, we reset the background grid after each loadstep to the original regular background grid. However, at step 5, the grid can be changed or a completely new grid introduced as required. The implicit MPM described above has recently been extended to include generalised basis functions, the GIMP approach of Bardenhagen and Kober, 10 and full descriptions and validations are given in Charlton et al. 23,24

Imposition of essential boundary conditions
Because calculations in the MPM take place using the background grid, it is there that boundary conditions should be applied. Imposition of an essential boundary condition where a displacement component is set to zero can be achieved by aligning the problem domain boundary with a grid boundary and maintaining that grid boundary in all the grids used. However, this is often inconvenient and is not possible for an arbitrary geometry or for where a nonzero essential boundary condition is applied. Figure 1A illustrates this using a rigid footing foundation problem, where the footing is modelled by imposing nonzero displacements on a part of the problem domain boundary. Figure 1B shows the problem domain discretised with MPs, surrounded by a grid of 4-noded quadrilateral elements. While one could align the grid to the two sides and base of the problem domain to impose zero essential boundaries there, this is impossible for the nonzero essential boundary at the footing. Many demonstrations of the MPM to date have avoided tackling this issue either by modelling problems where the boundaries are aligned to the grid (e.g. most slope stability modelling) or other workarounds (such as introducing stiff regions of MPs, for example, see Sołowski and Sloan 25 ). An example can be found in Sadeghirad et al 11 where the zero essential boundary at the end of a cantilever is simulated by adding a whole "ghost" cantilever behind the support, subject to the same loading and with the same discretisation, thus ensuring symmetry at the fixed end of the cantilever. This type of workaround can only be applied in special cases and is computationally tiresome, at least doubling the cost of any numerical simulation. It is clearly an impediment to the modelling of more challenging problems, such as a moving plough share in 3D or even just a curved boundary. A recent advance in this area has been the development of a "moving mesh concept" in the MPM, e.g. Phuong et al, 26 in which part of the background grid is updated to move with a moving essential boundary condition; however, this approach loses one of the advantages of the MPM, that a regular structured grid can be reused without change between load or time steps and it is also incapable of modelling a boundary which changes shape. Although it is possible to model some of the problems in this paper with an MPM-based on a triangular background grid (see Wieckowski 27 and Beuth et al, 28 among others), The use of a regular quadrilateral mesh has two important advantages, namely, (1) a regular quadrilateral mesh makes MP-element searches trivial, (2) some extensions to the MPM, for example, the GIMP method, 10 can only be applied to rectangular background grids.
The issue discussed above for the MPM is linked to the wider problem of imposing constraints on nonmatching meshes or the use of so-called embedded domains, for which there is a growing body of literature, e.g. Freytag et al, 29 Schillinger et al, 30 Hautefeuille et al, 31 and Ramos et al, 32 that has links to attempts to impose essential boundary conditions in meshless methods (Fernández-Méndez and Huerta 19 ), i.e., the choice of penalty, Lagrange multiplier, or Nitsche methods. Here, however, we take a simpler approach using the work of Kumar and coworkers 21,[33][34][35] in which the same issue is tackled (in part) for the case of standard FEs; i.e., a problem domain is discretised with FEs and essential boundary conditions must be applied on surfaces/edges that are not aligned with element edges. Kumar and coworkers start from an idea originally developed in the 1950s by Kantorovich and Krylov, 36 which is also discussed in Freytag et al. 29 The approach in some related studies 21,33-35 is termed an IB method and was presented for the case of fixed essential boundaries parallel to one of the coordinate axes for the 2D case.
Here we extend the general IB method to allow for fixed, roller, and prescribed essential boundary conditions at any inclination with respect to the coordinate axes and then demonstrate its implementation in the MPM framework. The development is given for the 2D case for clarity and to maintain links with some related studies. 21,[33][34][35] The main contribution of this work is equally applicable to standard FEs and the MPM. The IB approach is conceptually similar to the earlier work of Höllig and coworkers [37][38][39] in that they both directly modify interpolation functions to impose the essential boundary conditions. 40 The weighted extended B-spline method of Höllig uses an R-function 41 to weight the interpolation functions to impose homogeneous Dirichlet boundary conditions. However, a clear distinction between the work of Höllig and coworkers [37][38][39] and the IB method is that the latter allows inhomogeneous Dirichlet boundary conditions to be imposed independent of the discretisation. The very recent work of Zhang and Zhao 40 provides an extension to the IB method through generalising the boundary value function. They combine these advances with the finite cell method, and the level set method to represent the boundaries, to provide a general method for the imposing of inhomogeneous Dirichlet boundary conditions in two dimensions for thermo-elasticity. In this paper, we adopt the earlier approach of some related studies, 21,33-35 extend it to allow boundaries at any inclination, and combine the approach with the MPM to provide, for the first time, a general method to impose inhomogeneous Dirichlet boundary conditions. At this point, it is important to note that other methods also exist, which remove the link between a mesh and a problem domain feature, notably the eXtended finite element method, 42,43 where the shape functions are enriched to allow the modelling of a displacement discontinuity representing a crack within an element, but this is a different problem to the one tackled here.

INCLINED IMPLICIT ESSENTIAL BOUNDARY CONDITIONS
In this section, we first present the IB method given in previous studies, 21,33-35 which applies to boundaries parallel to one of the coordinate axes. We then demonstrate how it can be amended for boundaries at any inclination. The key idea in the IB method is to define trial and test function spaces, which will implicitly lead to the enforcement of the essential boundary conditions, i.e., for the trial and test functions, respectively, where have been introduced into Equations 13 and 14 to distinguish them from the standard FE displacement approximation. In Equation 13, {u} is the standard approximation for displacement in an FE (i.e., as given by Equation (4)) while {u a } is the essential boundary condition (i.e., zero for a fixed degree of freedom or nonzero for a prescribed displacement). At the implicitly defined boundary, the first term in Equation 13 will be suppressed via the matrix [D], enforcing the essential boundary condition in the second term. Substituting the trial and test functions into the weak form (Equation 1) leads to expressions for the stiffness matrices of elements containing essential boundary conditions (given in detail below).
The implicit essential boundaries are defined by signed distance functions j in  2 for the 2D case, where j < 0 indicates the exterior and j > 0 the interior of the problem domain for the jth boundary. On the boundary, j = 0. Essential boundary (or Dirichlet) functions d j ( j ) are then defined for each boundary, which are equal to 0 on the boundary and outside the domain; inside the domain, the function rises to unity at a small distance away from the boundary. d j are often simple C 0 continuous quadratic functions in j . At points where more than one essential boundary is active, i.e., at a domain corner, and then the product of the d( j ) forms the essential boundary function and in general, we write the net essential boundary at a point as where k refers to the component of displacement defined at that boundary. D k are then the components of the diagonal matrix [D] in Equation 13. Here, following the work of Kumar et al, 21 we adopt the following Dirichlet function where is the dimension over which the Dirichlet function d j (for a boundary condition j) varies from 0 to 1. A key omission from the original development of this idea, and one of the motivations for the current paper, is the means of modelling "roller" boundary conditions (zero displacement in the normal direction combined with zero traction in the tangential direction); however, inspection of Equation 13 reveals this to amount to setting D k to unity for the coordinate direction that is free. Figure 2 shows the 2D essential boundary matrices for the case of roller boundary conditions along two boundaries.

Stiffness matrices
Substituting the trial and test functions into the weak form and following the algorithm given above leads to the element stiffness matrix in Equation 7 where the effect of the IB conditions is to alter the form of the [B] matrix. 21,34 For a bilinear quadrilateral The latter is built from blocks for each element node j as [k E ] in Equation 7 can be expanded as follows: or more concisely as or more compactly as and and again, this can be compactly expressed as The calculation of these additional stiffness matrices can be simplified by recognising that along a boundary the Dirichlet functions and their derivatives are constant along a tangential coordinate and that outside the band of width contributions to [K 2 ] and [K 3 ] will be 0. If we term the coordinate directions as (n, t) as indicated in Figure 3A, then we can integrate within the band as For the case of a roller boundary permitting translation in the t-direction (eg, Figure 3A) then

Inclined essential boundaries
In Kumar et al, 21,34 IBs are confined to those parallel to one of the coordinate axes and no guidance is provided on the case of inclined implicit boundaries. In Figure 3B, the boundary is aligned so that the n-axis is rotated by from the x-direction. The variation of the Dirichlet function, d, with distance from the essential boundary is also shown. Note that, as defined by Equation 16, the function varies from 0 at the boundary to unity at a distance into the solid material. If full fixity was required along this boundary, then one would simply define D 1 and D 2 in the x and y directions and use the procedure given above, because full fixity implies no effect locally of boundary inclination. However, with the motivation mentioned above, to develop a model of ploughing, we have determined the amendments necessary to the IB method given above. Matrix [D 1 ] contains values of the Dirichlet function and therefore is invariant with respect to the coordinate system used for its calculation in the band of width . However, the components of matrix [D 2 ] are derivatives of the Dirichlet functions with respect to the (x, y) system and, therefore, require a transformation if the boundary is inclined. This can be shown to be a simple transformation matrix so that for any IB, inclined or not and It is worth noting that in most cases, components of [K 2 ] are much smaller than components of [K 3 ], because the nonderivative terms are much smaller than the derivative terms. It should also be highlighted that the distance function matrix, [D 1 ], does not need to be mapped between coordinate systems as it contains distances that remain unchanged by a change in coordinate, rather than direction components as in [D 2 ].

NUMERICAL IMPLEMENTATION
In this section, the numerical implementation for imposition of implicit boundaries in MPM is discussed. Firstly, we discuss the preprocessing stages required to set up an MPM problem, and secondly, we explain how local stiffness matrices are evaluated that cover both the material stiffness and the implicit boundaries. Consider a problem domain ( Figure 4A) subjected to prescribed, roller, and fixed essential boundaries. The domain is discretized into a set of MPs, each with an associated volume covered by the MP. A regular background grid of elements on which the governing equation will be solved is then set up (here linear quadrilateral elements are used). Material points are then inserted in every grid element, as shown in Figure 4B. The MPs can conveniently be positioned at standard Gauss quadrature locations and volumes calculated accordingly, but there is no requirement for this. In the next step, all MPs with volumes entirely outside the domain are removed, as shown in Figure 4C. Volumes are adjusted for any MPs with volumes partly outside the domain ( Figure 4D), and the MP position is relocated to the centroid of the remaining solid region (see Figure 4E). This process ensures that the total volume of the domain is captured by the MPs.
After discretizing the domain into the set of MPs, the local stiffness matrices [k E ] of every grid element are evaluated using Equation 18. This contains both the material stiffness matrix [K 1 ] and the implicit penalty stiffness matrices [K 2 ] and [K 3 ] that ensure the imposition of any essential boundary conditions on the element ( Figure 4F). If no essential boundary conditions coincide with the element, then [K 2 ] and [K 3 ] are zero. To assemble the global stiffness matrix, each element [k E ] is evaluated for elements containing MPs and any empty elements are ignored. For elements intersected by essential boundaries, [K 2 ] and [K 3 ] are evaluated in a 2-stage process. Firstly, the inner integrals of Equations 26 and 27 across the band are evaluated using a 3-point Gauss quadrature scheme. Once the inner integrals are evaluated, the outer integrals are evaluated along the boundary segment within the element using a 10-point scheme (see Figure 5). If a nonzero essential boundary condition is required, i.e., {u a } ≠ 0, the equivalent reaction force to cause this displacement {f E R } is determined as and is included in the linear system to be solved. Equations 26, 27, and 28 are valid for both curved and straight boundaries, and the inner integrals in 26 and 27 are evaluated in the same way in both cases. For curved boundaries, the outer integral is approximated by a series of quadrature points within each element where the normal and tangential directions will vary along the segment. In this way, the IB method provides a piecewise linear approximation to the edge of the domain.

NUMERICAL EXAMPLES
{u p } and {u e } are the MP and analytical (or exact) displacement vectors, n mp is the number of MPs in the simulation, V is the total volume, and ||{·}|| denotes the L2 norm of a vector. The stress error over the domain, r ij , is computed in a similar way from the errors at each MP, r p ij , that is, where p ij and e ij are the MP and analytical stresses. In the sections that follow, Equations 29 and 30 are used to report both the error at an MP and the total error across the domain.

Compression of a contained block
The first example investigates the ability of inclined roller boundaries, modelled as described above, to deliver accurate results and also examines the effect of the choice of bandwidth . The problem is the compression of a 100 mm by 100 mm block of material contained between three roller boundaries, where a prescribed displacement of 1 mm is applied to the fourth boundary in a single loadstep, a problem introduced in Cortis et al. 44 Figure 6 shows the problem at three inclinations, to test the transformation for inclined roller boundaries. In each case, the background grid is aligned to the (x, y) axes of the 0 • model and single MPs are generated in each grid element, as described in Section 4. The width of the band over which the integrations in 21 and     Symmetry is maintained in all of the displacement error distributions, and the discontinuous nature of the stress error is due to the low order (bilinear) element basis noting that we do not apply smoothing to the error distributions. Table 1 gives overall error values for each model. The smallest errors occur in the 0 • and 90 • models with higher errors seen in the 30 • and 45 • models, likely indicating integration errors from partially filled elements and truncated MP domains. It should also be noted that the results from the 0 • and 90 • models are identical, demonstrating the objectivity of the method when integration accuracy variations are removed. Convergence plots for displacement errors with respect to the bandwidth and the degree of refinement are shown in Figures 10, 11, and 12, where h is the size of the background grid elements. Based on conventional FE theory, the displacement error should converge at a rate of p + 1, where p is the polynomial order of the FE basis. In this case, using linear FEs (p = 1), an FE implementation using the same boundary approximation method would converge at a rate of 2 (as shown by the triangular convergence rate indicator on Figures 10, 11, and 12). However, in the MPM, additional errors are introduced  45 Therefore, we expect the convergence rate of the MPM to be bounded by that of the underlying FE basis. All models show similar uniform rates of convergence for displacement error showing errors decreasing with h-refinement, illustrating that the solution is unaffected by the presence of the IBs. As becomes smaller, the condition number of the global stiffness matrix should increase because of the increasing penalty contributions from [K 2 ] and [K 3 ], however, these plots indicate that this does not seem to affect convergence for the investigated range of . Small variations in the convergence rates can be observed, particularly at the coarser discretisations, however, these reduce as the simulations are refined and can be attributed to inaccuracies in the integration of the material stiffness due to partially filled background elements. In particular, note that these variations are not observed in the 90 • model; in this case the elements are fully populated by MP domains. Figure 13 shows the convergence for the 90 • problem as a function of the ratio of the bandwidth to the element size, ∕h. As expected, decreasing ∕h decreases the error in the simulation, and for ∕h ≤ 1 × 10 −6 (the grey region in Figure 13), the maximum observed error is around 1×10 −8 . The remainder of the numerical examples presented in this paper adopt a bandwidth size of = 1 × 10 −6 h.

A deep beam under self-weight
The second example is a deep beam acting as a cantilever loaded with body force (see Figure 14) and demonstrates the influence on integration errors in the MPM when using the IB method. The problem domain is 100 mm square with full fixity at the left-hand support implemented using the IB method. A vertically downwards body force of 1 N/mm 2 is applied in a single step. Two convergence studies are conducted: firstly, with respect to the number of MPs in each grid element (keeping the element number constant) and, secondly, with respect to the number of grid elements, with MPs per element kept constant. The

FIGURE 14
Deep beam with self-weight: problem layout background grid in the first study comprises a mesh of 10 × 10 mm elements, and because no analytical solution is available (as the support is full fixity throughout), the MPM solution is compared with a FE solution using the same mesh, a suitable Gauss quadrature rule, and standard imposition of essential boundary conditions, the goal being to demonstrate that the MPM approach with IBs converges to the standard FE solution for the same discretisation. The IB bandwidth used throughout this study = 1 × 10 −6 h. The accuracy of the MPM model is assessed by a pointwise r u on the y-displacement measured at point A in Figure 14, and results are presented in Table 2, and a convergence plot in Figure 15A. As discussed in the previous section, the convergence rate of the MPM is governed by the underlying FE basis. However, the convergence rate of the method can be reduced through inaccurate integration and projection between the MPs and the background grid. In this example, the MPM results are being compared with the finite element method on the same computational grid; we are therefore investigating the error caused by inaccurate integration due to nonoptimum MP locations.
In the second study, the number of MPs per element is kept constant at 16, the grid is refined, and the results are compared with those from standard FE analysis with the same mesh. Table 3 gives the same pointwise r u on the y-displacement measured at point A in Figure 14, and the convergence is plotted in Figure 15B. Both of these studies demonstrate that the convergence of the MPM (reducing the errors associated with nonoptimum quadrature locations) is not significantly influenced by the method of imposing the essential boundary conditions. Increasing the number of MPs per element or reducing the size of the background grid cells, both reduces the discrepancy between the MPM and FE results.

Simple shear
The third problem presented also investigates integration errors in the MPM with IBs but concentrating on partially filled elements. In this example, we consider a 100×100 mm under simple shear conditions. The problem is shown in Figure 16 where two 100 × 100 mm grid elements contain the domain discretised using 20 × 20 MPs. The domain has a fully fixed essential boundary (using the IB method) on the left-hand side and is shifted from x = 0 mm to x = 100 mm using the coordinate system shown. A shearing displacement of 10 mm is imposed on the free end of the cantilever in a single step, while the top and bottom surfaces are fixed in the horizontal direction (u x = 0) and allowed to move in the vertical direction only to achieve simple shear conditions. In a simple shear problem, only uniform shear stresses ( xy ) exist, and the analytical solution can be expressed as xy = E 1+ xy , where xy is the prescribed shear strain (the sheared distance over the length of the body). The shear stress error r u is plotted with respect to the IB position ( Figure 16). From the plot, it is clear that the IB method leads to successful enforcement of an Only when the IB moves very close to the edge of the grid element does accuracy drop; however, it should be acknowledged that in all cases the error is very small, being less than 2 × 10 −8 , and any changes are relative to these very small numbers.
As expected, the lowest error is observed when the domain is contained within a single element (x A = 0 or 100 mm). As the domain moves between the 2 elements, the error initially increases because of the poor integration of the element containing a few MPs clustered at one edge (the right-hand element in the figure). The error reduces with increasing translation and reaches a local minimum at x A = 50 mm (half the domain in the two elements) before increasing again because of the problem's symmetry. The reduction in error between x A = 10 mm and x A = 50 mm is due to the improved integration of the right-hand element as MPs fill more of the element. The rate of improvement reduces with translation because of the increased level of error in the left-hand element as MPs move across the boundary between the two elements.

A slender cantilever beam
We now consider a slender cantilever beam of dimensions 20 × 100 mm. The cantilever is oriented at three different angles, 0 • , 45 • and 90 • as shown in Figure 17, but the loading and geometry are effectively identical; only the discretisation changes. All models are fully fixed at one end using the IB method and are subjected to an incremental body force applied in the y-directions shown in Figure 17. A background MPM grid, having 2×2 mm elements is used, while the cantilever domain itself is discretised using 100 MPs in every grid element. Note that we adopt a large number of MPs per element to reduce the cell-crossing instability inherent in the MPM.
Displacement at the origin (the free end) is plotted against total body force for the 3 models in Figure 18. While the difference between the 0 • and 90 • cantilever models is negligible, the difference is more significant between 0 • or 90 • and 45 • cantilever models. This is likely due to poor numerical approximation where grid elements lack sufficient MPs or where MPs are

FIGURE 19
Cantilever: background grid h-refinement convergence between the 45 • and 90 • models concentrated in a region of the element. As before, h-refinement of the background grid can be used to reduce this error. The percentage pointwise error between the 45 • and 90 • models is computed and plotted against the cantilever depth d normalised by the grid size h (Figure 19). The respective deflected shapes of both models using the finest background grid are shown in Figure 20 when the full body force is applied. To further demonstrate the quality of results obtained, stress contour plots for the 90 • case at full body force are shown in Figure 21 on the undistorted configuration close to the support. The stresses are reported for a body force of 0.2 N/mm 2 with a grid size of 1 × 1 mm and 25 MPs per background grid cell in the original configuration. The problem is one where the only major gradients in stresses occur at the fully fixed end, at the top and bottom of the section, and it can be seen that these concentrations are produced faithfully using the method.

Rigid footing penetration
The penultimate numerical example is a demonstration problem of the penetration of a rough rigid footing into an elastic domain, as shown in Figure 22. The elastic domain is 3×3 m, and the background grid has an element size of h = 0.2 m. Each element initially contains 4 MPs, as shown in Figure 22B. The bandwidth for the IBM method for all of the boundaries is set to = 1 × 10 −6 h. A displacement of 1 m is applied over 20 loadsteps, and Figure 22C-F shows the displaced MP profiles for loadsteps 2, 10, 15, and 20. The MPs are coloured according to their vertical displacement relative to their original position. The problem demonstrates the ability of the method to model a moving displacement boundary condition that does not coincide with the background grid.

Funnel penetration
The final numerical example is that of an elastic domain being forced through a curved funnel and is intended to demonstrate the use of the method for problems that cannot be rotated such that the boundaries coincide with the background grid. The problem geometry is shown in Figure 23A where the roller boundary conditions on both the top and bottom boundaries and the prescribed displacement on the right edge are enforced using the IB method. A 4mm × 4mm background grid was used  Figure 23B where the size of the points are determined by the MP volume. The elastic body, with the same material parameters as in the previous sections, is forced through the funnel with an incremental displacement of 4 mm per loadstep. The bandwidth for the IB method for all of the boundaries is set to = 1 × 10 −6 h, and the curved boundary is approximated by piecewise linear segments over each background grid cell.
The deformation of the MPs for the 6 loadsteps is shown in Figure 23. The MPs have been coloured according to their x-displacement measured from their initial positions, where the blue points have moved the greatest distance. Despite the very large displacement increment (3 times the grid spacing) in each loadstep, as new MPs move towards the inclined boundary condition, their displacement trajectories align with the boundary. All of the MPs remain on the correct side of the inclined boundary condition throughout the simulation. The reduction in the domain height causes a region, coloured in red in Figure 23, of reduced x-displacement while a region towards the top of the problem, coloured in blue, increases in displacement. The right-hand boundary remains vertical with a uniform value of displacement irrespective of its position within the background grid.

CONCLUSIONS
The MPM provides a means of analysing problems in solid mechanics for which conventional FEs, and other numerical techniques, struggle, notably those involving very large deformations. Because problem domain boundaries do not necessarily coincide with the boundary of the grid on which MPM calculations are performed, an issue arises when wishing to model nonzero essential boundary conditions and roller boundary type conditions on a regular background grid. In this paper, we have presented a simple technique, similar to a more limited version previously applied to FEs, and shown how it can be adapted to model inclined boundaries and then demonstrated its use on problems analysed using the implicit MPM. To the authors' knowledge, this is the first time that a general method has been proposed that allows arbitrary inhomogeneous Dirichlet boundary conditions (zero and nonzero values at any inclination) to be imposed in the MPM. In particular, the method presented in this paper is different from other approaches in the MPM in that (1) the boundaries are independent of the background mesh, 26 (2) it does not require artificially stiff regions of MPs to be introduced, 25 and (3) the method does not rely on mirroring of the problem domain to impose symmetry as seen in other MPM approaches. 11 The method significantly extends the applicability of the MPM to examine problems that were previously unassailable using the method, the work itself is simple to extend to 3D, and the inclined boundary method can equally be applied to standard FEs.
An extension to this work would be to combine the IB method with traction boundary conditions and thereby allow any type of boundary condition to be modelled within the MPM. The work could also be extended to allow for frictional boundary conditions (limiting shear traction dependent on the normal pressure) that would provide a more realistic interface for soil-structure interaction problems.