Advances in Engineering Software

The Material Point Method (MPM) is a computational tool ideally suited to modelling solid mechanics problems involving large deformations where conventional mesh-based methods struggle. Explicit and implicit formulations are available, but for both the learning curve for understanding the method and arriving at a useful implementation is severe. Researchers must understand and implement finite element analysis, non-linear material behaviour, finite deformation mechanics and non-linear solution methods before they can even verify their formulations. This issue represents a significant barrier for post-doctoral researchers, graduate students and undergraduate students to start working with (and understanding) the method. This paper presents A Material Point Learning Environment (AMPLE) based around implicit variants of the method, with the aim of softening this steep learning curve via MATLAB-based, accessible and compact scripts. The code is freely available from github.com/wmcoombs/AMPLE.


Introduction
Mesh-based Lagrangian approaches such as the finite element method dominate the analysis of solid mechanics problems. However there are issues with these methods for problems involving very large deformation in that the discretisation used (i.e. the mesh) becomes distorted, leading to inaccurate results and in extreme cases, eventual breakdown of the numerics due to element inversion. Various attempts have been made to address this issue within mesh-based methods themselves, notably with Arbitrary Lagrangian Eulerian (ALE) methods [2] where a standard Lagrangian step determines material displacements and is followed by an Eulerian step to generate a new mesh for the deformed problem domain. Apart from the complexity of the approach, which is off-putting, the major error source in using ALE methods comes in remapping of variables between the two meshes [27] The Material Point Method (MPM) is an alternative to pure Lagrangian approaches and is well suited to problems involving very large deformations. The method was developed in the 1990s by Sulsky et al. [31] as a solid mechanics extension to the FLuid Implicit Particle (FLIP) method [4] which itself was developed from the Particle-In-Cell (PIC) method [15]. The basic idea of the method is to discretise a problem domain with material points (sometimes called particles) which carry information (mass, volume, stress, state variables, etc.). The information is then mapped to the nodes of a regular background finite element grid where calculations are carried out. The results are mapped back to the material points and the mesh is then discarded.
This cycle is then repeated, each step using a regular mesh to conduct the calculations thus avoiding any issues with mesh distortion, as translation is recorded only at the material points. The MPM has been used for a variety of solid mechanics problems such as fracture [16,23], contact [21], machining [13] and impact [34]. Recently the method has caught the particular interest of the geotechnical community for the modelling of landslides, where coupled displacement-pressure formulations of the MPM have been used to model saturated and unsaturated soils, examples include [29,35]. Outside solid mechanics it is interesting to note that the MPM has also been used to create animations for films [18,30].
In the standard MPM, material points contribute stiffness to the nodes of the grid in which they are currently located. When a material point translates from one grid cell to another (often termed "cellcrossing") the jump in stiffness can cause problems such as oscillations in the stress resultants and consequential instability in the solution. For this reason, MPMs have been developed where each material point "owns" a finite subvolume of the problem domain, which in some methods can itself change shape. These so-called advanced domain methods include the Generalised Interpolation material point method (GIMPM), [1] and the more recent Convected Particle Domain Interpolation (CPDI) methods [25,26].
The material point method was initially developed in an explicit format for purely dynamic problems where material stiffness is ignored and a simple forward difference solution obtained for the particle velocities over a number of time steps. Providing the mass matrix is diagonal, an explicit approach avoids the need to solve a linear system (as required by an implicit approach) however, there are well-known restrictions on the step size possible without loss of stability, and a lack of error control. The first implicit MPM appeared in [14], and a semiimplicit approach for dynamics was presented in [32]. In the latter the cost of forming the precise material tangent matrix was partially avoided with an approximation, and an iterative method was then used to solve the linear system of equations. A glance at citations for the landmark papers in explicit and implicit MPM (i.e. [14] and [31]) indicates a continuing preference for explicit MPM, however the advantages of an implicit approach when working in some areas of modelling, such as materials with sensitive material nonlinearity, has not been fully highlighted.
Guilkey and Weiss [14] highlight the clear links between implicit MPM and standard finite elements, and indeed, one way to summarise the implicit MPM is a finite element method where the integration points (material points) are allowed to move independently of the mesh.
The purpose of this paper is to introduce a Matlab-based "learning environment" for those wishing to explore the MPM, based on the authors' own experiences working with students and other researchers over the past five years.

Material point learning curve
Despite the clear links between the material point method and the finite element method, the learning curve for researchers to arrive at a useful material point method code is far more severe. The authors' interpretation of the learning curve for those attempting to understand, and implement, the material point method is shown in Fig. 1. The starting point for any researcher looking to use/investigate the material point method is to understand finite elements. However, even if a small strain assumption and linear material behaviour are adopted, the material point method is still a globally non-linear method; each time step will be linear but the overall response will be non-linear due to changes in stiffness as material points move through the background mesh. This means that, in order to obtain physically meaningful results, researchers are forced to adopt a geometrically non-linear analysis framework and understand both finite deformation mechanics and nonlinear solution methods (implicit or explicit). In addition to this, the types of problems that are normally investigated using the material point method (geotechnical failures, impact, fragmentation, etc.) also require non-linear material (or constitutive) models to be incorporated into the code. Once all of these ingredients have been assembled one can undertake useful material point analysis and start to add additional features, such as advanced boundary conditions, etc.

Development principles
AMPLE has been developed as an environment through which researchers, especially PhD/MSc students, can understand the material point method. The development was guided by the following principles: • MATLAB-based -as discussed in the previous section, the material point method is scientifically complex and if users/developers also have to understand thousands of lines of Fortran/C/C++ code the hurdle to its use may become insurmountable. AMPLE has been developed in MATLAB to remove, or at least significantly lessen, the syntax learning curve and allow researchers to concentrate on understanding the key elements of the material point method.
• Compact -AMPLE's implementation was inspired by Trefethen's [33] philosophy of ten digit algorithms "... a little gem of a program to compute something numerical: Ten digits, Five seconds, And just one page" The basic idea is that if a program is compact enough to be viewed on a single page it allows the structure of the whole algorithm to be understood and visualised, making mistakes and misunderstandings far less likely. Each of AMPLE's scripts/functions have been written so that they are readable on a single A4 page (or computer screen), this allows users to understand the structure of each of the code segments.
• Modular & expandable -AMPLE has been written as a series of compact functions called by a core analysis script. For example, the link between the material points and the background mesh is contained within one function and all of the other functions remain unchanged if the background mesh is changed. The continuum mechanics formulation is contained in another function and the material model (stress-strain relationship) in another. Although the initial AMPLE release is focused on two-dimensional analysis, all of the continuum mechanics is implemented for three-dimensional analysis making it straightforward to reduce or increase the code to one and three-dimensions, respectively. This allows users/  developers to quickly understand the purpose of each of the functions and adapt/modify/replace functions as required.
• Rigorous -all aspects of AMPLE have been verified and are based on a rigorous updated Lagrangian continuum mechanics framework [9]. The implemented algorithms are based on published material [5,6,11], including convergence analysis to ensure that the fundamentals of the code are sound. This should provide confidence when using/expanding AMPLE.
• Proof rather than performance -it is well documented that MATLAB is not as computationally efficient as compiled languages (such as Fortran or C/C++). However, the focus of AMPLE is on proof of new concepts and ideas, not high performance computations. It provides a environment for researchers to understand the material point method and test out new ideas and explore the impact of these changes on the performance of the method.

Paper layout
The layout of the remainder of the paper is as follows: Section 2 details, as concisely as possible, the material point formulation implemented within AMPLE; Section 3 describes some of the implementation aspects associated with this formulation, such as the basis functions of the implemented MPMs, amongst other details; Section 4 details AMPLE's code structure including function formats and data structures; Section 5 provides three demonstration cases that explore key aspects of the AMPLE environment; finally, Section 6 briefly concludes the paper. A mix of index and matrix-vector notation is used for the continuum formulation and numerical implementation aspects of the paper, respectively.

Material point formulation
AMPLE is an implementation of a quasi-static implicit finite deformation elasto-plastic material point method based on an updated Lagrangian formulation 1 defined by the following weak statement of equilibrium φ t is the motion of the material body with domain, Ω, which is subjected to tractions, t i , on the boundary of the domain (with surface, s), ∂Ω, and body forces, b i , acting over the volume, v of the domain, which lead to a Cauchy stress field, σ ij , through the body. The weak form is derived in the current frame assuming a field of admissible virtual displacements, η i . Note that AMPLE does not include tractions so the surface integral term in (1) is neglected in what follows. Unlike finite element methods, the imposition of tractions in the standard material point method is not trivial and requires a representation of the physical boundary to be constructed; see the work of Bing et al. [3] and Remmerswaal [24] for methods of imposing traction boundary conditions in the material point method.

Finite deformation mechanics
In finite deformation mechanics, the deformation gradient, F ij , provides the fundamental link between the original and deformed configurations where X i are the original (reference) coordinates and = x X t ( , ) i i are the updated coordinates in the current (deformed) body, where φ is the motion of the body. It is assumed that the deformation gradient can be multiplicatively decomposed into elastic and plastic components [19,20] where the superscripts e and p denote the elastic and plastic components. In this paper we adopt logarithmic strains and Kirchhoff stresses and combine these measures with an exponential map of the plastic flow rule to allow the use of conventional small-strain stress integration algorithms with a finite deformation framework. This is a powerful combination as it allows existing constitutive formulations to be used directly rather than reformulating them for the particular choice of stress and strain measures used in the large deformation mechanics [17,28]. Within this formulation, the elastic logarithmic strain is defined as is the left elastic Cauchy-Green strain and the Kirchhoff stress, τ ij , can be obtained using where is the linear elastic stiffness matrix. The Cauchy stress can be obtained from the Kirchhoff stress through is the volume ratio between the deformed and reference configurations. In order to advance the non-linear solution, the finite deformation equations are discretised in pseudo-time by imposing the deformation over a number of load (or pseudo-time) steps. This allows the current deformation gradient to be defined using where ΔF ij is the increment in the deformation gradient between the previously converged state, denoted using a subscript n, and the current state. In order to obtain the updated Kirchhoff stress state for the current deformation gradient, a constitutive model requires an initial estimate (or trial) of the elastic strain (or stress) state. In this approach the trial elastic Cauchy-Green strain tensor is given by where the subscript t denotes a quantity defined in the trial state. The previous elastic Cauchy-Green strain tensor, b ( ) , n ij e can be obtained from the previous elastic strain state through and the trial elastic strain state follows as The adopted constitutive algorithm can then be used to return the updated elastic strain, , ij e and Kirchhoff stress, τ ij , states.

Constitutive formulations & stress updating
AMPLE includes two different constitutive models: (i) Isotropic linear elasticity and (ii) isotropic linear elasticity with perfect plasticity and a von Mises yield surface with associated plastic flow. In the case of isotropic linear elasticity, the updated elastic strain, , ij e is equal to the elastic trial strain, ( ) , t ij e and the Kirchhoff stress is given by (5). The elasto-plastic von Mises model is integrated using an implicit elastic predictor, plastic corrector algorithm. The von Mises yield function can be defined as 1 AMPLE adopts the same formulation as implemented in the generalised interpolation approach of Charlton et al. [5] but has been included here for completeness.
and ρ y is the yield strength of the material. Within this algorithm the elastic trial strain, ( ) , t ij e acts as the initial estimate for the updated elastic strain state. If the corresponding trial Kirchhoff stress is inside the von Mises yield envelope (f < 0), then the material is undergoing elastic behaviour and the updated elastic strain is equal to the trial state. If the trial stress state is outside of the yield envelope then it must be corrected back onto the yield surface. AMPLE uses a backward Euler procedure to perform this correction, details of which can be found in [6], amongst others.

Discrete material point formulation
The Galerkin form of the weak statement of equilibrium over each background grid cell, E, can be obtained from (1) as where [∇ x S vp ] is the tensorial form of the strain-displacement matrix containing derivatives of the shape functions with respect to the updated coordinates. 2 The first term in (12) is the internal force within an grid cell and the second term is the external force vector, in this case the body forces. AMPLE adopts a standard implicit Newton-Raphson (N-R) procedure to solve (12) for the unknown nodal displacements. This requires (12) to be linearised with respect to the nodal displacements to arrive at the stiffness of each cell in the background mesh km lm is the spatial consistent tangent modulus for a point within the background grid cell (see Charlton et al. [5] for details).
In material point methods the physical domain is discretised by a number of material points. These points are used to numerically approximate the stiffness (13) of the grid cells in the background mesh, essentially replacing the conventional Gauss-Legendre points (or other integration method). The key difference between material point and finite element methods is that these integration points move relative to the background mesh rather than being directly coupled to the positions of the background grid nodes. The stiffness contribution of a single material point to the background mesh is where V p is the volume associated with the material point in the spatial (updated) frame where V p n and V p 0 are the volumes associated with the material point in the previously converged state and the initial configuration, respectively. Comparing (14) with conventional finite element literature, the volume of the material point, V p , replaces the Gauss-Legendre quadrature weight and the determinant of the Jacobian that maps between the local and global frames. The product of the determinant of the Jacobian and the Gauss point weight provides the physical volume associated with the Gauss point which is equivalent to V p in the material point method. The internal force contribution of a single material point to the background mesh is Following the work of [5,10], the increment in the deformation gradient is obtained from where Δu i is the displacement increment within the current load step, X j are the coordinates at the start of the load step and n is the number of nodes that influence the material point. This allows the increment in the deformation gradient to be obtained from derivatives of the basis functions based on the coordinates of the nodes at the start of the load step. The spatial derivatives of the basis functions can subsequently be calculated using the method proposed by Charlton et al. [5], that is It is essential that the spatial derivatives are used in the strain-displacement matrix, [∇ x S vp ], to both ensure the correct order of convergence in the N-R process and convergence towards the correct solution based on the internal force contribution, (16), of each material point [5,10]. It is not necessary to map the basis functions between the start of the load step and the current configuration because it is assumed that, during a load step, the displacement of a material point, (Δu p ) i , is linked to the nodal displacements through the basis functions, S vp , that is where n is the number of nodes that influence the material point and (Δu v ) i are the displacement of the nodes over the current load step. Therefore the basis functions evaluated at the start or the end of the load step are identical (see [9] for a detailed discussion on this point).

Implementation aspects
This section presents key information regarding the implementation of AMPLE. This includes the basis functions for both the standard and generalised interpolation material point methods, boundary conditions and material point position/domain updating. An outline of the computational procedure is also given at the end of the section.

Basis functions
The basis functions are one of the key areas where the material point method and the finite element method diverge. For the standard material point method the basis functions of a material point are simply the basis functions of the underlying finite element grid. In this case, one way to obtain the basis functions (and the spatial derivatives of these functions, if required) is to calculate the local position of the material point within its associated background grid cell and then evaluate the basis functions as in the finite element method. However, if the background grid comprises a regular grid with the cell boundaries aligned with the global coordinate system, it is straightforward to evaluate directly the basis functions of a material point. This is the approach taken in this paper as the initial release of AMPLE assumes that the background grid consists of regular two-dimensional bi-linear quadrilateral grid cells with their edges aligned with the global Cartesian coordinates. With this assumption, the basis functions for the both the standard and generalised interpolation material point methods can be expressed as the convolution of material point's domain with the basis functions of the underlying finite element grid, that is where Ω p is the influence domain associated with the material point, χ p is the material point's characteristic function which defines the zone of influence (or domain) of a material point, N v are the underlying shape functions of the finite element grid which are dependent of the position of the material point at the start of the load step, X p . The basis function is only given in one-dimension; the extension to higher dimensions is obtained through the Cartesian product of the shape functions in each direction. 3 The gradient of the basis function can be expressed as For multi-dimension problems the gradient of the basis functions are obtained from the product of the gradient in one direction with the shape function in the other direction, for example the derivative of the basis functions with respect to X in a two dimensional problem is The following sections provide the basis functions for the standard and generalised interpolation material point methods.

Standard interpolation
The basis functions for the standard material point method are obtained by replacing the characteristic function, χ p , with a Dirac delta function. With this substitution, (20) becomes where h is the size of the background grid (distance between the nodes in each direction) and X v is the position of the node (or vertex) associated with the basis function at the start of the load step. The gradients of the basis functions with respect to the material point position are

Generalised interpolation
The particular form of the generalised interpolation material point method adopted in this paper assumes a unity characteristic function (a hat function with the value of one inside the material point's domain and zero outside) of length 2l p centred on X p . 4 This characteristic yields the following basis functions 1˜: 1˜: As with the basis functions in (25), when a material point's domain lies entirely within a background grid the gradient of the basis function equals that of the standard finite element function.

Boundary conditions
In this paper we ignore the external tractions as their general implementation within material point methods is complex. Dirichlet (essential/displacement) boundary conditions are imposed directly on the background mesh in the same way as the standard finite element method. This imposes the restriction that essential boundary conditions must be imposed along parts of the problem domain that are contiguous with background grid cell boundaries, as with the majority of MPMs presetned to date. Deviating from this requires special treatments such as those developed in [12]. The nodal body forces in (12) are approximated using { } is the body force associated with the material point. For two dimensional analysis, 1}, where ϱ is the material's density, g is gravitational acceleration and = m V p p is the mass associated with the material point. Point forces, if required, are held and convected with the material points. They are mapped to the background grid using

Non-linear solution procedure
The nodal displacements within a load step, {Δd}, can be obtained by iteratively updating the nodal displacements until (12) is satisfied within a given tolerance using k is the current iteration within the load step, [K] is the global stiffness matrix and {δd k } is the iterative increment in the displacements from that iteration. The global stiffness matrix is obtained through assembling the individual material point contributions, that is  where A is the standard assembly operator acting over all of the material points in the problem and [k p ] is the stiffness contribution of a single material point. f { } k R 1 is the residual out of balance force vector associated with the previous displacement value, i.e. the difference between the internal forces due to the stresses within the material and the applied boundary conditions, as given by (12). The global residual vector is assembled through where the first term is the internal force vector and the external force vector, {f ext }, which is constant over the load step, is given by

Position and domain updating
At the end of each load step the material point positions, volumes and (if using the generalised interpolation material point method) domain half-lengths, l p , should be updated. The updated positions of the material points at the end of the load step are given by

Computational procedure
The steps in the implemented algorithm are concisely summarised below. The applied body force and/or point forces are split into a number of load steps and for each of these steps the following process is adopted: 1. calculate the stiffness contribution, [k p ], of all of the material points using (14) and assemble the individual contribution of each material point into the global stiffness matrix, [K]; 2. calculate the internal force contribution, {f p }, of all of the material points using (16) and assemble the contributions into the global internal force vector, {f R }, in (12); 3. solve for the nodal displacements within a load step, using the N-R process (29) until the out-of-balance force converges to a specified tolerance; 4. update the material point positions, volumes and domain lengths through interpolation from the incremental nodal displacements, deformation gradient and stretch tensor using (33), (15) and (34); 5. reset the background grid.
These steps are displayed graphically in Fig. 3. Fig. 4 and comprises three loops:

AMPLE's high-level structure is shown in
1. a loadstep for loop that first determines the external forces at the nodes for the current step and then calculates the nodal basis functions, S vp , and spatial derivatives, S ,

X vp
for each material point. The code then enters a while loop to find equilibrium between the internal and external forces. Once equilibrium has been obtained the material point positions and internal variables (deformation gradient, stress, elastic logarithmic strain, domain size, etc.) are updated. 2. a Newton-Raphson while loop to solve the global equilibrium equations on the background mesh. The loop first solves the nonlinear system of equations to determine the displacement and reaction force increments at the nodes before moving onto the determination of the current internal force and stiffness via a material point for loop. 3. a material point for loop that determines the internal force and stiffness contribution of each of the material points to the background mesh.
AMPLE's main script, ample.m, is shown in its entirety in Fig. 5. The format of the code aligns with the algorithm shown in Fig. 4. Comments are shown in green and have been truncated on the right hand side for clarity of the executed code. Only two of the three loops described above can be seen in the main ample.m file: (1) A loadstep for loop spanning lines 9 through 33 and (2) a Newton-Raphson while loop over lines 21 through 30. The material point loop is contained within the detMPs.m function on line 25.
The main ample.m script calls eight functions which are described in Table 1. The purpose of each of these functions is explained below: setupGrid: called on line 3 of ample.m, the function returns the analysis-specific information such as details of the background grid and any information held at material points (position, material properties, etc.). This information is held in two structured arrays, mpData and mesh, that are explained in Section 4.2. In addition to this the function returns the total number of load steps, lstps, and the gravitational acceleration, g, which are returned as variables independent of the material point and mesh data.   Table 3) mesh -mesh structured array (see Table 4)  It is important to note that the majority of AMPLE's functions do not depend on: (i) The form of the background mesh, (ii) the type of material point methods adopted or (iii) the number of physical dimensions (1D, 2D plane strain or 3D). The exceptions to this are: elemMPinfo.m which depends on the background mesh type and the material point variant and updateMPs.m which depends on the material point variant as different updating procedures are required for different material point methods. Obviously setupGrid.m depends on all three of the above points as it contains the analysis-specific information that will change depending on the user's requirements.
The functions are organised into a number of folders depending on their purpose, specifically: constitutive contains the constitutive functions, setup contains the functions that are required to provide the analysis-specific information, plotting contains the post-processing files and functions contains the remaining .m files. AMPLE's file structure also contains output and documentation that contain the generated VTK files and AMPLE's html documentation. Table 2 lists all of the variables used in the main ample.m script along with their mathematical symbol, size and description. The 24 variables are listed in terms of their order of appearance in the ample.m script. Note that mesh and mpData are structured arrays containing the mesh and material point data, respectively, and are

Function format
All of AMPLE's functions share a common file format. Fig. 7 shows an example AMPLE function, specifically the function to determine the external forces on the background mesh based on body and point forces at material points, detExtForce.m, which is called on line 12 of ample.m. The function starts with a common comment block that includes the following information: (i) brief description/title; (ii) author and date information with a more detailed description of the function's purpose; (iii) function call format; (iv) required input information; (v) the function's output; and (vi) a see also section detailing the functions that are called by the function (in the case of detExtForce.m, none).
All of AMPLE's functions support the MATLAB help command, for example typing help detExtForce will return the top comment block shown in Fig. 7.

Data structures
The majority of the analysis information required by AMPLE is stored in two structured arrays: mpData: this structured array contains material point information, such as the point's position, deformation gradient, Cauchy stress, etc. The 21 fields within mpData are detailed in Table 3, along with their mathematical symbols and dimensions per material point. A number of the quantities depend on the number of physical dimensions, nD. mesh: this structured array contains information about the background mesh, such as the positions of the nodes and the topology of each of the background mesh cells. The five fields of mesh are detailed in Table 4, where nels, nodes and nen are the number of background grid cells, the total number of nodes and the number of nodes per background grid cell, respectively.
MATLAB's structured arrays provide a convenient way to store the material point data as different material points will potentially influence different numbers of nodes and therefore will require different amounts of storage for, for example, the basis functions that influence the point and their spatial derivatives. 6 It is worth highlighting two of the fields that provide different options within AMPLE. mpType controls the material point type used in the analysis: mpType=1 results in a standard material point whereas when mpType=2 the point uses generalised interpolation basis functions. cmType controls the constitutive model used by the material point: cmType=1 specifies the isotropic linear elastic model contained within Hooke3d.m, whereas cmType=2 calls the linear elastic-perfectly plastic von Mises constitutive model, VMconst.m. Both of these could be extended to include other material point variants, such as CPDI methods, and different material models.

Demonstration cases
The analyses presented in this section explore changing a number of AMPLE's features/options. The physical problems modelled and the features modelled are as follows: 1. Compaction of a two-dimensional column under self weight: changing numbers of material points, material point types and mesh density; 2. A large deformation elastic beam subject to an end load: visualisation of material point data; and 3. Elasto-plastic collapse of a two-dimensional body: changing the constitutive model.
In all cases the global tolerance on the normalised out of balance force in the Newton-Raphson process was set at × 1 10 , 9 where the normalised out of balance force was defined as ||( · )|| denotes the L2 norm of ( · ), {f ext } includes the external forces applied to the problem (body forces, tractions, point loads) and {f react } are the reaction forces due to the prescribed displacement boundary conditions.   * depends on the problem analysed; number of displacement constraints. 6 It is noted that structured arrays are not the most computationally efficient way to store data in MATLAB, however AMPLE's focus is on proof (and clarity) rather than performance.

Convergence: Compaction under self weight
The first example is the one dimensional compaction of an elastic column with an initial height of = l 50 0 m under its own self weight. The material has a Young's modulus of 10kPa and a Poisson's ratio of zero. The background mesh is comprised of square background grid cells 7 with roller boundary conditions on the base and sides and the column is discretised by a 2 × 2 grid of equally spaced material points in each initially populated background grid cell. A body force of 800N/m 2 ( = g 10m/s 2 and an initial density of = 80 0 kg/m 3 ) is applied over 40 equal loadsteps.
The basic setup information for this analysis, as contained on lines 57 through 67 of setupGrid, is shown in Fig. 8. The figure shows the setup information for an analysis with 64 grid cells in the y-vertical direction (nelsy) and a single grid cell in the x-direction (nelsx). The physical body is discretised by generalised interpolation material points, mpType=2 (material point Type), with linear-elastic material behaviour, cmType=1 (constitutive model Type).
The analytical solution for the normal stress in the y-direction for this problem is     7 Note that the width of the problem is changed depending on the number of grid cells in the vertical direction to ensure that the cells remained square.
(nelsy) was varied between 2 2 and 2 13 , that is between 4 and 8,192, in powers of 2 while maintaining the same material point/cell ratio 8 It is clear that the standard material point method (mpType=1) does not converge with uniform h refinement due to cell-crossing errors whereas the generalised interpolation method (mpType=2) converges at an approximately linear rate. Fig. 9 also shows AMPLE's run time with progressive background mesh refinement 9 For most of the analyses, the run time scales approximately linearly with the number of material points ( × 4 10 2 seconds per material point). However, when the number of grid cells in the y-direction exceeds 4,096 the time spent in the linear solver starts to dominate, with a corresponding increase in the gradient of the run time. Interestingly the generalised interpolation analyses typically have a slightly lower run time due to the Newton-Raphson process taking fewer iterations to find convergence.

Visualisation: Large deformation elastic beam
The second demonstration problem is the large deformation bending of an elastic cantilever beam subjected to a point load at its free end using the generalised interpolation material point method   Fig. 8. 8 The number of material points in each grid cell was kept constant, therefore the total number of material points increased from 16, for the 2 2 mesh, to 32,768, for the 2 13 mesh. 9 All analyses were conducted using MATLAB R2017b within a macOS V.10.14.5 environment on a 2.3GHz Intel Core i7 with 16GB of RAM.
boundary conditions are applied at = X 0 and the beam is fully-fixed at its mid-depth (ly-d/2); lines 86-90: Storage of the background mesh information within the mesh data structure;   solution and there is very little difference in the force-displacement predictions of the two discretisations. Fig. 12 also shows the deformed material point positions at the end of the = h 0.125m analysis that have been coloured according to the normal stress in the y-direction, σ yy . The figure was produced from the VTK output files that are generated on line 32 of the main ample.m script (shown in Fig. 5) via postPro.m (shown in Fig. 13) and saved into the output folder. 10 The postPro.m script generates two different types of VTK files: mpData_#.xtk: On line 5 where # is the current load step number (lstp). A VTK file based on the converged state at the end of each load step, including the following material point data: current position, total displacement between the current and initial positions and the Cauchy stress (all components). The script could be extended to include other data as appropriate, such as the volume at each material point, the level of straining, the deformation gradient, etc. AMPLE generates one mpData_#.xtk file per load step. mesh.vtk: On line 7 which includes the background mesh information. AMPLE assumes that the background mesh is unchanged through the analysis and therefore only a single mesh.vtk file is generated and this happens on the first load step (lstp=1). If the mesh was different between loadsteps, a series of mesh.vtk files could be generated in a similar way to the material point data.
The number of material points per background cell is a point of debate in the material point method literature. In the authors' experience 2 6 material points per element for a two dimensional problem represents a reasonable balance between stability (in terms of conditioning of the global stiffness matrix), accuracy (in terms of reducing quadrature errors) and efficiency. In order to demonstrate this point, Fig. 14 presents the normalised global force versus displacement response for = h 0.25m with different numbers of material points per background grid cell, alongside the analytical solution shown by discrete white-filled circles and squares. The analysis with 4 material points per cell fails to converge in the 8th load step due to ill conditioning of the global stiffness matrix -the final stable load step for this analysis is highlighted by the grey-filled circles. The analysis with 3 2 and 4 2 material points per cell have very good agreement with the 6 2 analysis for the vertical displacement and converge to the 6 2 analysis for the horizontal displacement. The increase in accuracy is very minor between the analysis with 3 2 material points and that with 6 2 ; if efficiency is a user's primary concern then 3 2 material points per background grid cell is a reasonable choice.

Material model: Elasto-plastic collapse
The final demonstration example is the elasto-plastic collapse of a rectangular body of material subject to a gravitational body force using the generalised interpolation material point method (mpType=2). Due to symmetry only half of the body is modelled and the initial discretisation of the problem is shown in Fig. 15 Fig. 15. These Dirichlet boundary conditions are imposed directly on the background grid by eliminating the degrees of freedom associated with these constraints from the global system of equations when solving the linear system during each N-R iteration. The degrees of freedom associated with the nodes that do not have a contribution from any material points are eliminated in the same way.
The material is modelled using a linear-elastic, perfectly-plastic von Mises constitutive formulation with the yield function defined in (11). In order for AMPLE to adopt this material behaviour it is necessary to set cmType=2 for all of the material points in the analysis. For this analysis, the yield strength of the material is = 20 y kPa, the Young's modulus is 1MPa and Poisson's ratio is 0.3. These material properties are stored in mCst for each material point with the mpData structured array (see Table 3) in the following order [E, ν, ρ y ]. A body force of 10 kN/m 3 (density of 1000 kg/m 3 and gravitational acceleration of 10 m/s 2 ) is applied over 40 equal loadsteps. Fig. 16 shows the deformed material point positions at the end of the analysis for = h 0.5 m and = h 0.25 m. In both cases, 6 2 generalised interpolation material points (cmType=2) were included within each initially populated background grid cell; a total of 16,384 and 65,536 material points for the = h 0.5 m and = h 0.25 m analyses, respectively. The material points are coloured according to σ xx , σ yy and σ xy , where the blue and red regions show areas of low and high stress, respectively. There is very little deformation for the first 20 loadsteps (up to a body force of 5 kN/m 3 ) but beyond this value the material rapidly collapses and achieves the final deformed state shown in Fig. 16. Table 5 presents the horizontal extent and maximum height of the collapsed body at the end of the analysis for the generalised interpolation material point method with different numbers of material points and background mesh sizes. The results from an updated Lagrangian finite element analysis using bi-linear (q 1 ) and bi-quadratic (q 2 ) quadrilateral elements are shown for comparison 11 . The generalised interpolation material point method results are between the linear and quadratic finite element results which is consistent with the basis functions of the method being linear if the material point is contained within an element and quadratic if the material point contributes to multiple elements. In this case the analyses using 3 2 and 6 2 material points give similar results.
All of the examples presented in this paper have assumed quasistatic conditions -any inertia effects are ignored. However, in the case of elasto-plastic collapse it could be argued that dynamic effects would impact on the physical process. We acknowledge this point, however the focus of AMPLE is on quasi-static analysis. An interesting future extension to this framework would be to dynamic problems via an implicit or explicit time stepping procedure.

Conclusion
The MPM is a relatively new approach for the modelling of large deformation solid mechanics problems and is attracting interest from a diverse range of applications. However, those interested in exploring its capabilities have to navigate their way through a number of challenging  Table 5 Elasto-plastic collapse: Horizontal and vertical extents at the end of the analysis, where nPs denotes the number of integration or material points in each element/cell. For the finite element analyses the q 1 and q 2 denote bi-linear and bi-quadratic quadrilateral elements, respectively.  aspects in order to understand the method and then have to turn the understanding into code. AMPLE is software for quasi-static implicit MPMs and has been developed to take some of the pain out of this exploration process. It provides a robust framework in which to explore MPMs and to develop new features appropriate for specific applications. AMPLE also provides a useful teaching tool for the MPM suitable for graduate students and post doctoral researchers wishing to consider using MPMs. The code is available here [7] and more information is available on the AMPLE web pages [8].