A Discontinuous Galerkin Material Point Method for the solution of impact problems in solid dynamics

Abstract An extension of the Material Point Method [1] based on the Discontinuous Galerkin approximation (DG) [2] is presented here. A solid domain is represented by a collection of particles that can move and carry the fields of the problem inside an arbitrary computational grid in order to provide a Lagrangian description of the deformation without mesh tangling issues. The background mesh is then used as a support for the Discontinuous Galerkin approximation that leads to a weak form of conservation laws involving numerical fluxes defined at element faces. Those terms allow the introduction of the characteristic structure of hyperbolic problems within the numerical method by using an approximate Riemann solver [3] . The Discontinuous Galerkin Material Point Method, which can be viewed as a Discontinuous Galerkin Finite Element Method (DGFEM) with modified quadrature rule, aims at meeting advantages of both mesh-free and DG methods. The method is derived within the finite deformation framework for multidimensional problems by using a total Lagrangian formulation. A particular attention is paid to one specific discretization leading to a stability condition that allows to set the CFL number at one. The approach is illustrated and compared to existing or developed analytical solutions on one-dimensional problems and compared to the finite element method on two-dimensional simulations.


Introduction
A wide variety of engineering problems in solid dynamics, such as impact problems or high speed forming techniques are modeled with hyperbolic initial boundary value problems (IBVP). The solutions of this class of mathematical problems involve waves which are important to accurately follow in order to understand physical phenomena occurring within the medium during the deformation. Indeed, knowledge of waves paths provides information about the loading history undergone by each particle, which in turn allows to properly compute the residual strains and stresses, and to optimize the final shape of the body. Furthermore high speed forming processes such as electromagnetic forming [4] often involve finite deformation than can complicate the development of an efficient numerical scheme. A Lagrangian mesh approach may require re-meshing and projection steps in order to deal with mesh entanglement effects. Those issues can be avoided by using an Eulerian method. However interface tracking techniques and convection steps are required in order to follow the boundaries and transport internal variables which is less convenient for solid mechanics.
The numerical simulation of IBVP has been and is still widely performed in solid mechanics with the Finite Element Method (FEM) [5] due to, among other assets, its ability to easily deal with nonlinearities such as historydependent constitutive models [6]. However, despite its well-known advantages a Lagrangian finite element formulation becomes less efficient and accurate for problems involving very large strains due to difficulties mentioned above. The Material Point Method (MPM) developed in the 90's [1] can be viewed as an extension of the FEM in which quadrature points can move in an arbitrary computational grid. Those material points represent a continuum with a Lagrangian description and are used to store all the fields and internal variables of the problem. The underlying mesh is only used to compute gradients, apply boundary conditions and solve discrete equations. The method hence uses many transfer steps of fields between material points and nodes. Since the data are stored at material points, the grid can be discarded at each time step and re-constructed for computational convenience without additional projections, thus avoiding mesh entanglement. In addition to problems caused by finite deformations, classical explicit time integrators used in solid dynamics with FEM and MPM introduce high frequency noise in the vicinity of discontinuities. The removal of these spurious numerical oscillations with additional viscosities is difficult to achieve without loss of accuracy, and can be troublesome for the wave tracking. Hence, although the material point method has been successfully applied to a wide range of complex large strain engineering problems [7] this approach seems unable to accurately deal with discontinuities owing to numerical oscillations, diffusion introduced by mapping steps, and its low maximum admissible Courant number.
On the other hand, the Finite Volume Method (FVM) [8] initially developed for fluid dynamics has proven to be a very efficient tool in order to accurately follow the path of waves. Indeed, the formulation involves a conservative form leading to the same order of accuracy for velocity and gradients, and numerical fluxes allowing to account for the characteristic structure of hyperbolic equations. Moreover, several works on FVM provided high order schemes satisfying the Total Variation Diminishing (TVD) stability condition [9] and recent studies extended the approaches to solid mechanics for problems involving history-dependent models [10,11], or finite deformations with a Lagrangian formulation [12,13].
The discontinuous Galerkin (DG) approximation [2] enables to build numerical schemes that benefit from both finite element and finite volume methods. DG-methods are based on piece-wise polynomial shape functions and their combination with the finite element formulation (DGFEM) led to several development steps to reach Total Variation Diminishing in the Means (TVDM) and Total Variation Bounded (TVB) high order schemes that ensure convergence to entropy-satisfying solutions [14]. In addition, the same order of accuracy is reached for velocity and gradients within a finite element framework when the weak form is based on a conservation laws system. Although this approach can easily handle mesh-adaption strategies due to the relaxation of fields continuity, the mesh tangling problems do not vanish. Furthermore, the Courant condition that applies is restrictive [14] and prevents to accurately capture discontinuities.
The purpose of this work is to develop a numerical method gathering the strengths of both material point and DG methods and avoiding their respective limits by means of the discontinuous Galerkin approximation. The computational grid of the MPM is used to build a piece-wise linear approximation of the solution of conservation laws, which weak form is written on each element. Boundary integrals arising from integration by parts involve numerical fluxes that connect the elements together and are computed with an approximate Riemann solver [3]. Projection steps between the underlying mesh, used to solve the discrete equations, and material points that store the fields of the problem, provide a Lagrangian description of the deformation in an arbitrary computational grid. The resulting Discontinuous Galerkin Material Point Method (DGMPM) hence aims to enable an accurate wave paths tracking in solid media undergoing large strains. A similar approach based on the Particle-In-Cell (PIC) method has already been developed in [15] for electromagnetism equations. This discontinuous Galerkin version of PIC, unlike the DGMPM, uses classical Gauss points to integrate the weak form. Moreover, electric and magnetic fields are stored at nodes whereas charge and current densities are stored at particles which prevent the arbitrariness of the grid. In what follows, after a brief review of continuum equations (section 2), the DGMPM formulation is derived within the large strains framework with a total Lagrangian approach for hyperelasticity in section 3. Section 4 is dedicated to the computation of interface fluxes arising in the weak form. First, Riemann problems are written at interfaces with a quasi-linear form which spectral analysis is performed. Second, an approximate Riemann solver that can take into account transverse corrections for multi-dimensional problems is derived. Finally, the enforcement of boundary conditions is discussed. Next, after a summary of the numerical procedure in section 5, the stability analysis of the one-dimensional scheme is carried out in section 6. This analysis shows that a reduced integration of the weak form provides an improvement of the stability condition that limits the DGFEM. At last, the method is illustrated in section 7 on one-dimensional cases in order to compare it to analytical solutions, and two-dimensional problems for which comparisons are performed with the Finite Element Method.

Conservation equations
We consider a three-dimensional isothermal deformation of a homogeneous continuum elastic body from the reference configuration Ω 0 ∈ R 3 , to the configuration of the same body at a different time Ω t ∈ R 3 . The motion of the body is described by the smooth mapping ϕ(X, t) ≡ x from the reference coordinates system X to the spatial coordinates system x. The second-order two-point deformation gradient tensor and the material velocity vector are respectively defined as: where the superposed dot denotes the material time derivative, Greek indices stand for material coordinates and Latin ones for spatial coordinates. Combination of equations (1) and (2) leads to the Lagrangian kinematic equations: The solid must locally satisfy the conservation of linear momentum written in the reference configuration within the time interval of interest t ∈ τ = [0, T ]: where ρ 0 (X) = ρ (X, t = 0) is the reference mass density, Π is the first Piola-Kirchhoff stress tensor, ρ 0 b is the body forces vector and ∇ X · (•) is the divergence vector with respect to reference coordinates. For simplicity, body forces are neglected. The use of the Cartesian coordinates system leads to a homogeneous system of partial differential equations written with equations (3) and (4) [16]: where U is the vector of conserved quantities and F α is the flux vector in the α th direction: In what follows, vectors (resp. tensors) denoted by typographic symbols belong to R p (resp. R p × R p ) where p = 2, 3, whereas calligraphic symbols denote column arrays.

Constitutive model
We are interested here in modeling hyperelastic materials characterized by a local objective stored-energy function Ψ. In order to ensure mathematical well-posedness of the problem this function is assumed to be polyconvex [17] namely, Ψ is convex with respect to the deformation gradient, its determinant J and its co-factor tensor H = JF −T . Such a stored-energy can then be expressed as a function of the deformation gradient only, leading to the classical definition of hyperelastic materials: The fourth-order tangent modulus tensor C and the second-order acoustic tensor A are respectively defined as: where N α are the components of an arbitrary material normal vector N.

Space-time discretization
As for MPM, a continuum body Ω t is discretized into a set of P material points in an arbitrary Cartesian grid. The grid is here supposed to be made of E non-overlapping elements, or cells, for convenience. In this grid, the set of faces separating cells that contain particles and having empty neighboring elements defines the boundary of the mesh (see figure 1 for a two-dimensional case). In this section, the formulation of the new approach is developed within the large strains framework.
boundary of the domain

Weak form of the system of conservation laws
We seek an approximate solution to system (5) in a finite-dimensional function space using discontinuous Galerkin approach. The key idea of DG methods is to allow jump of fields across mesh elements faces by using broken polynomials spaces [2,18,19]. In what follows, we will only consider linear polynomials and the corresponding broken spaces: where H 1 (Ω e ) and P 1 (Ω e ) denote respectively Sobolev space and the linear (Lagrange) polynomials space, in Ω e .A weak form is written by multiplying equation (5) by a test function V ∈ V 1 h and integrating over each element e of the mesh: Integration by parts leads to: where N is the outward normal vector to the boundary, and F · N = F N denotes a normal flux at cell faces. Those numerical fluxes can be computed from the solution of an approximate Riemann solver [20] (see section 4). In particular, the stationary solution of Riemann problems defined at element boundaries Γ i such that ∪ i Γ i =Γ e gives the well-known Godunov's method. Specific quantities are introduced and delta Dirac distributions are used as particles characteristic for the mass density [1]: where m p is the mass of the p th material point. Such a discretization of the reference mass density combined with the writing of conservation laws (5) leads to a total Lagrangian formulation, and equation (12) thus reads: which, with Dirac integration properties transforms into:

Semi-discrete and discrete equations
The particles are viewed as interpolation points so that fields of the weak form (16) are evaluated from nodal values: where S ip denotes the shape function attached to the i th node evaluated at X = X p , andŪ i the specific vector of conserved quantities at node i. In equation (17) and in the following of the paper, index p stands for a quantity defined at material points whereas (i, j) stand for nodal quantities. Combination of equations (16)- (17) and arbitrariness of the test field yield: which is the semi-discrete system that must be solved on the grid. Note that the consistent mass matrix M ij can be singular when only one material point lies in an element (reduced integration) [21]. However the use of a diagonally lumped mass matrix permits to avoid the singularity problem, and in what follows the lumped matrix will be written: Furthermore this configuration is equivalent to the DGFEM [22] with a mass matrix integrated with only one Gauss point.
Finally, the time interval τ is discretized in N T sub-intervals of size ∆t n such that N T n=1 ∆t n = T and time integration is performed with an explicit forward Euler algorithm leading to discrete equations: where the superscripts denotes the time step at which a field is evaluated.

Computation of interface fluxes
Interface fluxes involved in boundary integrals of equation (18) can be used to introduce the characteristic structure of hyperbolic problems within the numerical scheme, providing that they satisfy some mathematical requirements, namely they must be E-fluxes [23]. One possibility is to solve a Riemann problem normal to each boundary of a given cell and to compute the Godunov flux of the solution. An approximate Riemann solver that computes a solution which can be used to calculate those fluxes is presented in this section.

The Riemann problem
For the multi-dimensional case the problem at a given face consists in finding the solution of the following auxiliary initial value problem (IVP): where the top equation is the projection of system (5) along the normal vector to an interface with the auxiliary variable ξ = X · N. Initial conditions (Ũ L ,Ũ R ) arise from the DG approximation which introduces multi-valued fields on interfaces and can be viewed as a duplication of nodes (figure 2). Due to the expressions of flux and conserved quantities vectors (6), the calculation of Godunov fluxes based on the solutionŨ to problem (21) requires the computation of constitutive equations. However, this procedure can be costly for nonlinear problems, especially for history-dependent models. Hence, the quasi-linear form is written introducing an auxiliary vector Q [24]: A wise choice of the auxiliary vector, namely Q is a rearrangement of the flux components, permits to compute the interface flux from the solution Q * of problem (22) with a simple assembly: In what follows, the auxiliary vector will be taken as Q T = Π T , v T . One single Riemann problem is solved at each interface center with initial conditions obtained by averaging the nodal values on each side as depicted in figure 2, so that the computational cost is reduced. Figure 2: Duplication of nodes at an interface and building of initial conditions of the Riemann problem (2D). IVP (22) admits self-similar solutions, constant along each ray traced from the origin of the (ξ, t) plane. A linearized Riemann solver allows to determine explicitly at a moderate cost an approximate value of the stationary solution [3] providing that the characteristic structure of the solution is known.

Spectral analysis
The characteristic structure of the solution of problem (22) is given by the spectral analysis of the Jacobian matrix defined for the general three-dimensional case as: where 0 2 and 0 4 are null second-order and fourth-order tensors respectively. In what follows, the dependency of the Jacobian matrix on Q is omitted to simplify the notation. The problem (21), and hence (22), is known to be hyperbolic if and only if [25]: which is known as the Legendre-Hadamard or ellipticity condition. The prerequisite that the stored-energy function is polyconvex ensures the satisfaction of the Legendre-Hadamard condition (25) [17] hence, hyperbolicity is guaranteed. Moreover, the Jacobian matrix (24) has in this case 6 non zero and 6 null eigen-values with associated independent eigen-vectors [25]. The M eigen-values of the Jacobian matrix (24) give rise to M characteristic fields given by the resolution of the following eigen-system: R K Q being the right eigen-vector associated to the K th eigen-value λ K of J Q N . Each of these characteristic fields corresponds to one wave emerging from the origin of the (ξ, t) plane. Equation (26) leads to the two independent systems: and by introducing equation (28) in (27), one gets the acoustic tensor eigen-system: where (ω p , Y p ) are the acoustic tensors eigen-values (ω = ρ 0 λ 2 ) and associated eigen-vectors. Since Legendre-Hadamard condition (25) is satisfied, the acoustic tensor admits three positive eigen-values leading to six characteristic speeds λ 2p, 2p−1 = ± ω p /ρ 0 ,(p = 1, 2, 3). We then see that one eigen-vector Y p of the acoustic tensor is associated to two eigen-values λ K and two velocity vectors v K , thus it comes: Finally, the stress components of the Jacobian eigen-vectors Π K are determined with equation (27). Then, we have for the non zero eigen-values the following characteristic fields: The last step consists in building a complete basis of the kernel of J Q N by finding six independent vectors solving equation (28) for the null eigen-value:

Approximate-state Riemann solver
The solution of problem (22) consists of M discontinuous waves emanating from the origin of the (ξ, t) plane, from which the stationary solution reads: where W m = δ m R m Q are the waves formed of the wave strength coefficients δ m weighting the right eigen-vectors R m Q associated to the M eigen-values λ m of the Jacobian matrix (24). By rearranging equation (32), one can compute the wave strength coefficients by solving: and then calculate the stationary solution Q * and associated normal fluxes F N with assembly (23). Those fluxes are assumed to be constant on the face and the boundary integrals of the semi-discrete equation (18) can be computed.

Transverse Riemann solver
The method derived above for the computation of normal fluxes can be viewed as the donor-cell Upwind (DCU) method [8] meaning that only contributions from the upwind cells sharing an edge (in two dimensions) with the current one are considered. For multidimensional problems waves can travel in several directions such that contributions coming from corner cells must be taken into account in order to improve accuracy and stability of the numerical scheme. The corner transport upwind (CTU) method [26] consists in considering contributions coming from upwind cells sharing only a node (in two dimensions) with the considered grid cell and propagating in bias. This approach allows to improve the Courant condition especially for solid mechanics problems for which strain components are coupled through Poisson's effect. One can define at each cell interface left-going and right-going fluctuations defined as: Let's consider the patch of grid cells shown in figure 3, and focus on the edge denoted (i) which local normal vector  (21) defined on the adjacent edge, hence the name transverse Riemann solver. Spectral analysis of corresponding Jacobian matrix is carried out in [27]. The negative normal fluctuation is, for instance, decomposed on the characteristic basis associated to edge ( j) as: where R m j is based on the normal vector N j but also on different tangent moduli between grid cells L and its neighbor T. The calculation of transverse fluctuations is performed by solving equation (35) for coefficients β m , and by multiplying each wave W m j by the associated characteristic speed. Since only waves with positive celerities with respect to the orientation defined by outward normal vector to considered edge will contribute to the transverse fluctuation, only the positive operator B + is used: An additional numerical flux defined at edges is hence built from these transverse fluctuations: which contributes to the grid cell T (∆X j being the length of edge ( j)). With duplicated nodes introduced by the DG approximation this contribution must be counted negatively for nodes belonging to cell L (outgoing fluctuation) and positively for nodes belonging to cell T (incoming fluctuation).

Boundary conditions enforcement
As in MPM [21,28], boundary conditions (BCs) are treated at nodes that compose the domain boundary (see figure 1). Introducing ghost nodes on boundary interfaces allows to widen the use of the Approximate Riemann solver. Similarly to finite volumes [8], a vector Q G is enforced on those nodes so that the stationary solutions of Riemann problems (22) are consistent with boundary conditions. This vector is determined for one ghost node by (i) setting the free quantities (i.e: velocity (resp. stress) components for Neumann (resp. Dirichlet) BCs) equal to that of the associated internal node and (ii) solving system (32) for enforced quantities on ghost nodes.
Notice that the use of conservation laws (21) requires to inverse constitutive equations in order to enforce Neumann boundary conditions, which is not convenient. Hence, the use of the quasi-linear form (22) also allows to easily handle boundary conditions. Moreover, the transverse Riemann solver needs to introduce corner ghost nodes equivalent to finite volume corner cells [8].

Numerical procedure
Let us assume that the vector of specific conserved quantitiesŪ n is known at every material points that discretize a continuum body Ω in a grid made of N nodes nodes at a given time t n . We can now derive the computational steps to obtain the DGMPM solution, which requires to (i) solve the hyperbolic system (5) in Ω with associated boundary conditions and (ii) update the vectorŪ n on material points.
The scheme has been established with a total Lagrangian formulation therefore, the weak form (16) is written on the reference configuration. Hence, the lumped mass matrix M L i and the pseudo-stiffness matrices K α ij are computed once and for all at the beginning of the calculation: Then, the procedure follows these steps: U known at material points (a) Convective phase: the discrete equation (20) and Riemann problems (22) require a projection of fields onto the grid: to be solved for eachŪ i and Q i respectively. This weighted least squares interpolation is similar the one introduced in FLIP (Fluid Implicit Particle method) to avoid a diffusion-like error [1,29,30].
(b) Compute time step: in the general case the tangent modulus depends on the deformation gradient as well as the waves speeds. It is therefore needed to compute the time step so that the fastest wave can at most cross the smallest cell of the mesh according to Courant condition.
(c) The specific flux vectorsF α,i involved in the discrete equation (20) are computed from Q i knowing ρ 0 , thus avoiding the computation of constitutive equations .
(d) Enforce boundary conditions on mesh nodes (see section 4.5).
(e) Computation of interface fluxes: 1-Build Q ξ ± states based on Q n i where i denotes the nodes belonging to the face on side L or R.
2-Compute the stationary solution by means of the approximate Riemann solver.
3-Calculate the corresponding normal flux.
(f) Advance solution in time by solving the discrete equation (20) at each node.
(g) Map the updated solution to material points with a classical interpolation: At this point the mesh has virtually moved, but since fields have been transferred back to particles, the underlying grid can be discarded and rebuilt for the next time step for computational convenience, thus involving the convective phase (a) at the next time step. Notice that adaptive grid algorithms can be applied in the reference configuration, so that to improve wave front tracking in the current one. allows to increment the deformation ϕ(X, t) and to update stress components which will be used in the auxiliary vector for the next time step, through hyperelastic constitutive equations: This summary highlights significant differences with respect to the MPM. First, while the DGMPM uses a classical interpolation in order to map back the updated solution to material points (step f), FLIP method and MPM require an additional time integration on material points. Indeed, in the original schemes the changes in grid values (i.e: accelerations and strain rates) are mapped back to particles rather than the new value in order to avoid numerical dissipation [30]. Furthermore, the use of conservation laws (5) instead of momentum equation in the weak form implies that both velocity and gradients are solved at nodes making this new approach close to finite volume methods. Finally, the deformation gradient is no longer calculated with shape functions gradients so that the task of choosing between Update Stress First (USF) and Update Stress Last (USL) algorithms [31,32] vanishes. The USL formulation consists in integrating the constitutive equations at material points with the updated nodal velocity field given by the resolution of the MPM discrete equation thanks to shape functions gradients. Within the USF algorithm, the material points stress tensors are computed with the nodal velocity field resulting from the convective phase at the beginning of each time step. Since the stress tensors at particles are used in the discrete equation through the internal forces vector, these two formulations modify the accuracy of the MPM. In that sense the DGMPM scheme is simpler though additional cost is linked to Riemann solutions. Fortunately, the use of discontinuous Galerkin approximation makes this numerical method highly parallelizable [14].

Stability analysis for one-dimensional linear problems
A stability analysis of the one-dimensional DGMPM scheme is performed on a scalar linear hyperbolic equation. Following [33], the discretized equation is written in a finite difference sense in order to perform a von Neumann linear stability analysis.

Model equation -Space discretization
We consider the scalar linear advection equation for an arbitrary quantity q = ρq moving at the constant speed a ∈ R + * in a homogeneous one-dimensional medium of length l: The specific flux function isF = aq and equation (43) is discretized with the discontinuous Galerkin material point method. Thus, it is assumed that the medium has been divided with P material points arbitrarily distributed in E two-nodes elements of constant length ∆X ( figure 4). The grid mesh is such that at least one particle lies in every cell during the computation in order to ensure that there is no hole in the bar. Moreover, periodic boundary conditions are considered to simplify the analysis. Figure 4: One-dimensional mesh made of E elements of constant length ∆X = l E .
Since fields are carried by particles, we seek for the scheme equation that gives the solutionQ n+1 α at material point α and time step n + 1, with respect to the solutions of other particles at previous time step. In this section, Latin and Greek indices denote respectively nodes and material points, S iα is thus the shape function of node i evaluated at the position of material points α. Linear shape functions defined in element k yield: Furthermore, the cell containing a given particle β will be denoted by c(β) and this notation will be also used for nodes, for instance: c(2k − 1) = c(2k) = k.

DGMPM scheme equation
The method followed to write the scheme equation is to trace backward the numerical procedure described in section 5 in order to get an expression for material point α: where f is a function representing the DGMPM scheme. Quantities at time t n+1 are obtained by interpolating nodal solutions of the discrete equation (20) thanks to mapping (17): Interfaces fluxes are F N = (aq * )N, where q * is the stationary solution of Riemann problem at cell interface and N = ±1 the outward unit normal. Nodal solutions at time n + 1 in equation (46) are the result of the discrete system (20): Provided linear shape functions, the lumped mass and the pseudo-stiffness matrices are: where m α is the mass of α th material point. Since the pseudo-stiffness matrix component K ij is zero if i and j are not connected by an element, due to the discontinuities of shape functions, the sum over j in equation (47) reduces to: . We can thus write: Since waves travel from left to right (a > 0), the stationary solution is equal to q L on an interface, that is: Moreover, in a homogeneous medium the mass density can be replaced by ρ = N c(α) m α /∆X where N c(α) is the number of particles in the cell that contains the α th material point. Then, we have: Equations (51) and (54) are introduced into discrete system (47) so that equation (46) is rewritten as: In equation (55) the solutions at nodes result from the convection step (39): Thus, introduction of mapping (56) in equation (55) and permutation of sums over µ and i lead after some simplifications to the scheme equation: This recurrence relation can be further simplified by viewing that the term in parentheses is non zero and equal to one, due to partition of unity, only for element i = c(µ). On the other hand, S 2i−1,α S 2i−2,µ is non zero only if cell i contains material point α and is such that i = c(µ) + 1. Hence, the scheme equation reads: which can be written for simplicity:Q The numerical solutionQ n+1 α depends on the Courant number a∆t/∆X and on the position of material points within the grid mesh through the shape functions in the expression of coefficients D αµ . Furthermore, the first (resp. second) brackets in equation (58) involve shape functions that are non zero if material point µ and α lie in the same cell (resp. previous cell). This property means that the numerical domain of dependence of the DGMPM for scalar linear advection with a > 0 covers two cells regardless of the number of material points.

The von Neumann linear stability analysis
The computational domain is now repeated periodically by mapping it to the domain [−l, 0] so that the solution at material point β and time step n is expanded into a discrete Fourier basis of 2E + 1 harmonics over the domain with A n j , the magnitude of the j th harmonic at time step n, I = √ −1, and k j a wave number. Introduction of this expansion in equation (59) yields: The amplification factor between two time steps at a given point is defined as: A necessary condition to ensure the stability of a numerical scheme is that the amplification factor must be lower than one in modulus: A n+1 /A n < 1. This upper bound allows to ensure that eventual errors do not increase during the computation. For expression (62), this leads to: where the triangle inequality, and the unit modulus of the complex number e I(µ−α)k j ∆X have been used. Hence, the DGMPM scheme is stable for a given discretization if the Courant number is set so that the following condition is satisfied: According to the scheme equation (58), such a stability condition can be very hard to deal with analytically for general discretizations. However, the particular configuration for which one single particle lies in each cell can be studied analytically. In this case, condition (64) reads: which is satisfied for all a ∆t ∆X ≤ 1 regardless of the positions of particles. Hence, this space discretization leads to a stable scheme for the optimal Courant condition while the classical DGFEM scheme developed in [22] is restricted to condition ∆t/∆X = O( √ ∆X). Recall that this DGMPM space discretization can be viewed as a DGFEM with reduced integration of the weak form (see section 3). Hence, combination of reduced integration with projections between mesh nodes and integration points allows to improve the restrictive stability condition of original DGFEM [22,34]. Note that by considering space-time DG formulations [35,36], it is possible to relax constraints of pure space DGFEM and obtain a critical CFL at 1. This optimal CFL condition for DGMPM is however no longer valid for discretizations involving more than one material point per cell. Indeed, figure 5 shows the evolution of the amplification factor with respect to the Courant number for several (regular) discretizations. Since the scheme equation (58) depends on the particle considered, those curves are plotted for each material point contained in a given cell. In figures 5(e) and 5(f) we can see that the amplification factor satisfies the stability condition (64) for CFL numbers lower than one. Two material points in each cell are considered in figure 5(e), with either a symmetrical or an offset configuration, . It also appears that the more restrictive condition on CFL number is given by the upwind (first) material point (CFL = 0.43 for symmetrical configuration and CFL = 0.34 for shifted one). Moreover, the material points positions influence the stability condition as shown by the shifted configuration for which the result is more restrictive. Similar observations can be made with figure 5(f) in which three particles are considered. In that case, Courant number decreases to CFL = 0.29 for configuration 5(c) and CFL = 0.26 for configuration 5(d).
The von Neumann linear stability analysis performed in this section provides an explicit (nonlinear) condition (64) that must be satisfied in order to ensure the stability of the DGMPM scheme. First, this relation allows to fully exploit the ability of the method to rebuild the grid mesh arbitrarily. Indeed, after such a procedure, the number and positions of material points in grid cells can change and one must be able to adapt properly the CFL number in order to ensure stability. An advantage on the original MPM is hence highlighted since no stability condition exists for the latter scheme. Second, it has been shown that some DGMPM discretizations provide an improvement of the restrictive CFL number that applies to the DGFEM scheme. In particular, the optimal condition allowing to capture discontinuities (CFL = 1) can be reached when one particle lies in each cell of the grid mesh. This property turns out to be a strength for the DGMPM since it aims at follow waves accurately as shown with the following numerical results.

Numerical results
In this section, comparisons between DGMPM, MPM and FEM are performed. The first simulation considers one-dimensional linear elasticity and is purely used for analysis since MPM has little interest for small strains. Next, finite deformations are considered for one-dimensional and two-dimensional problems.

Compression impact on a linear elastic bar
We consider a bar made of a linear elastic material with density ρ 0 = 7800kg.m −3 , Young's modulus E = 2×10 11 Pa and arbitrary cross section S , suddenly subjected to a compression stress on its left end T d = 1 × 10 9 Pa while the right one is fixed. We assume that the impact is such that the linearized geometrical framework holds and mass density is constant ρ(x, t) ≡ ρ 0 (X). The analytical solution of this problem is known [37] and consists of an elastic disturbance propagating rightward in the bar at constant speed c = E/ρ before reflection on the right end. The length of the bar is l = 6 m and it is discretized with 150 regular grid cells containing one single centered material point. A comparison is performed between the DGMPM, the MPM and analytical solution. The MPM has been used for this comparison with the USL and USF implementation in order to see which procedure gives the most suitable results for impact problems. The Cauchy stress profile along the bar given by the three numerical schemes are plotted in figure 6 at two different time steps. Figure 6  fits perfectly the analytical profile due to the condition CFL = 1 while the MPM-USF, which suffers from diffusion and oscillations, can not capture a discontinuity though the same Courant number is used. On the other hand, the USL formulation of the material point method is unstable for CFL ≥ 0.7, although this result is not proved analytically. The scheme is then unable to capture discontinuities. However, the oscillations resulting from the MPM-USL procedure decrease with time and the correct stress level is obtained faster than for the USF solution. Furthermore, numerical diffusion is quite similar for the two MPM algorithms. Hence, the USL implementation of the material point method seems to be better suited to impact problems than the USF. It is of interest to see the influence of the number of material points lying in each element. Thus, we consider now the same problem as before, with 150 material points regularly spaced in 75 cells as space discretization. The particles are placed symmetrically with respect to cell centers ( figure 7). In that case, the optimal stability condition no longer holds. A comparison is performed between the same numerical methods as before and the analytical solution. In  CFL number introduces more diffusion than in the MPM solutions but no spurious oscillations appear. On the other hand, MPM solutions are quite similar to that of the previous discretization, though the USF algorithm leads to more oscillations. Furthermore, figure 8 shows that using a second order time integrator within the DGMPM scheme allows to obtain the optimal stability condition and hence, to fit the analytical solution. However, this property no longer holds when material points are moved in each cell so that the symmetry with respect to cells centers disappears. Thus, a condition equivalent to (64) is required for second order time discretization. This will be the object of future works.

Impact problems on a hyperelastic one-dimensional medium
In this section we consider a one-dimensional medium made of a hyperelastic Saint-Venant-Kirchhoff (SVK) material in order to see how the discontinuous Galerkin material point method behaves for nonlinear hyperbolic problems. The following deformation gradient tensor occurs in this medium: The Saint-Venant-Kirchhoff constitutive equation leads to the following PK1 tensor (see Appendix A): where Π= 2µ+λ 2 F(F 2 − 1), Π 2 = λ 2 (F 2 − 1) and (µ, λ) are the Lamé's coefficients. The SVK model is known to not be polyconvex however, the problem treated remains hyperbolic if condition (A.10), which depends on the loads, holds. Hence, we will only consider loading conditions that satisfy hyperbolicity condition. Moreover, the flux function resulting from this constitutive model is concave (− ∂ 2 Π ∂F 2 < 0) and leads to monotonically increasing characteristic speeds with the deformation gradient. Since the evolution of wave celerities with respect to the deformation gradient governs the wave pattern (i.e: either a rarefaction or a shock wave), the structure of the solution will differ from other constitutive laws with convex flux function such as the neo-Hookean model. As a consequence, the Saint-Venant-Kirchhoff model can lead to non-physical solutions. Nevertheless, the use of SVK model in spite of its limitations is here justified since an analytical solution can be found (see Appendix A), which is not possible for a neo-Hookean material.

Compression load
The one-dimensional medium is first submitted to a compression load on its left end T d = 2 × 10 10 Pa so that the solution of this problem consists of a right-going rarefaction wave as shown in Appendix A. The domain is discretized with 150 material points lying in 150 uniform grid cells and a comparison between DGMPM and MPM is performed. Figure 9 shows the stress profiles and the analytical solution along a horizontal line at different time steps. Slight oscillations appear in the DGMPM solution at the beginning of the computation ( figure 9(a)). However, these oscillations are much lower than the MPM ones. We then see that the DGMPM stress profile is close to the analytical solution but does not fit it exactly. Indeed, the approximate Riemann solver does not account for the  information contained in the rarefaction fan and cannot provide the exact solution. Figure 10 shows the error made on Riemann problems at cells interfaces, characterized by the distance between the intersection of integral curves on the one hand, and the intersection of approximated curves on the other hand. The 1-integral curve (resp. 2-integral curve) corresponds to all the states U connected through a shock wave (resp. rarefaction wave) to left (resp. right) initial condition at the left end of the bar. The stationary state corresponds to the intersection of those curves in the phase plane (F,ρ 0 v) and we see that the approximate and analytical stationary solutions are different. One can imagine to use an exact Riemann solver based on the analytical solution in order to reduce this error. However, such an implementation requires to solve a nonlinear problem at each cell interface and is very costly. Moreover, an analytical solution is available for SVK material but it is not the case for other constitutive models, which prevents the use of an exact solver.

Tensile load
We now consider a tensile load of same magnitude that gives rise to a shock for which the celerity of the fastest wave is given by the deformation gradient on the upwind side. The MPM solution of this problem shows the same limitations as before, it is hence not presented. The Figure 11 shows results obtained for the same discretization as before with the DGMPM compared to the analytical solution. We see that slight oscillations appear in the DGMPM solution after the shock. Note that those oscillations (figure 11(a)) decrease with time ( figure 11(b)). However, these oscillations can be removed with a Courant number smaller than one. Moreover, the numerical scheme do not perfectly capture the discontinuity, though a good behavior is shown, due to the wave celerity used for the time step evaluation which is an upper bound of the shock speed ω s (A.13). Nevertheless, the loading conditions considered so far are deliberately extreme for the purpose of visualization. A lower tensile load will lead to a numerical solution closer to analytical results as can be seen in figure 12 which shows the comparison for a load magnitude T d = 2 × 10 10 Pa.

Impact problem on a hyperelastic two-dimensional medium
Let's now move to multi-dimensional simulations by considering an infinite medium in directions e 2 and e 3 , and of width l in direction e 1 (see figure 13), undergoing a tensile load. The domain is modeled by a finite medium, with  In addition zero out-of-plane velocity and strain components are assumed so that plane waves can be simulated. The solid is discretized with 61 × 11 material points such that only one particle lies in each grid cell. The simulation has been performed with the DGMPM-DCU, since BCs reproduce a plane wave, which is compared to the analytical solution. Figures 14(a) and 14(b) show the comparison between numerical and analytical results for longitudinal and transversal stresses Π 11 and Π 22 , with a tensile load set to T d = 4 × 10 9 Pa. First, top figures in subplots show the longitudinal stress isovalues within the medium. The results have been post-treated with the software Paraview [38] which filters allow to represent material points with spheres deliberately big for the purposes of visualization. We see in those stress maps that the one-dimensional dependency of the solution is not disturbed within the numerical scheme. Next, bottom graphs in subplots 14(a) and 14(b) show analytical and numerical stress components along the bottom boundary of the medium. As for finite volume and finite element methods, two-dimensional simulations require a decrease of Courant number with the DGMPM. Since no two-dimensional stability analysis is available so far, the CFL number was set to 0.7 by testing various values. Hence, the shock wave is described less accurately than for one-dimensional computations. However, no oscillations appear in the numerical solutions and the analytical  stress levels are obtained for both components.

Impact problem on a hyperelastic plate
Previous simulations showed that the DGMPM provides good results for one-dimensional problems compared to an analytical solution. In this section, an extension to more general two-dimensional problems is made, for which no analytical solution is available, thus comparison with a finite element solution (bilinear approximation) is performed. The two-dimensional plate is made of the compressible hyperelastic neo-Hookean material, which is known to be a polyconvex model, and considered here in the plane strain case. Geometry, boundary and loading conditions, and material parameters are given in figure 15(a). The solid is discretized such that DGMPM material points are equivalent to finite element nodes, thus the plate is represented with l×h ≡ 33×25 material points and l×h ≡ 32×24 Q1 elements for FEM (see figure 15(b) for a coarse example). Once again, only one material point lies in each cell of the arbitrary computational grid. The finite element computation is performed with the software Abaqus [39], using an explicit time discretization with no artificial viscosity added . These numerical results are compared to those obtained from DGMPM computations using both CTU and DCU methods. The Courant number is set to unity in DGMPM-CTU and to 0.7 in DGMPM-DCU leading to average time steps ∆t CTU = 1.57 × 10 −5 s and ∆t DCU = 1.13 × 10 −5 s, whereas the constant time step used in the FEM simulation is ∆t FEM = 1.43 × 10 −5 s. Figure 16 shows numerical results in terms of the Cauchy stress tensor σ = 1 J ΠF T exported from Abaqus to Paraview with the code developed in [40]. Numerical Cauchy stress components (σ 11 ,σ 22 ) along bottom boundary are gathered in the left plots and longitudinal stress σ 11 isovalues in the domain resulting from CTU and FEM computations are compared in right colored maps for various time steps. Note that the deformed shape of finite element method is used as a background for the CTU stress map in order to compare shapes at different time steps. At the beginning of the computation ( figure 16(a)), stress profiles are quite similar despite slight oscillations in the FEM solutions and more numerical diffusion in the DCU solution than in the CTU one. Then, interactions between shear and tension-compression waves occurring in the left part of the plate give rise to a decrease of the stress level on the wave front which is managed differently by the two methods. Thus, additional oscillations are introduced in the FEM solution leading to a more diffuse wave profile (right colored map in figure 16(b)). On the other hand, the deformed shapes of the plate computed with the CTU and the FEM remain close. After the reflection of waves on the fixed boundary, significant differences can be seen in numerical solutions (figure 16(c)). We also see that final shapes differ, which can be explained by the velocity fields plotted in figure 17.
Those curves show the horizontal component of velocity along bottom boundary at different times. First, figures 17(a) and 17(b) show that spurious oscillations also appear in FEM velocity solution before reflection of the wave on fixed boundary. Those oscillations remain in the solution until the end of the computation ( figure 17(d)). On the other hand, the prescribed velocity boundary condition at the right end of the plate is not enforced exactly in the DGMPM solutions as shown in figure 17(d). Indeed, the enforcement of boundary conditions in MPM is still a challenging question [28] which does not spare the DGMPM. However, this point has not been investigated yet. Nonetheless, boundary conditions do not disturb much waves profiles as seen with all the simulations presented above. Hence, the final DGMPM stresses plotted in figure 16(d), as well as previous ones, do not exhibit overshoots. The absence of such over-stresses in hyperelastic regime is expected to allow a better management of history-dependent constitutive models including threshold such as elastoplasticity. In those models, non physical overshoots can lead to a wrong assessment of residual stresses and strains.

Conclusion
A new discontinuous Galerkin material point method was presented for the numerical simulation of impacts on hyperelastic media undergoing finite deformations. This approach is based on the weak form of a system of conservation laws written in a solid domain represented by particles in an arbitrary computational grid, as in the original material point method. The background mesh is used as a support for the discontinuous Galerkin approximation leading to a weak form element-by-element that involves numerical fluxes at cells faces. These are computed from an approximate Riemann solver leading to the donor-cell-upwind method. In addition, transverse corrections can be added for multidimensional problems so as to improve accuracy of the scheme by adapting the corner-transport-upwind method. The space discretization is based on particles mass density and an explicit time discretization is performed leading to the discrete equations. The stability analysis of this numerical scheme has been carried out in section 6 leading to an explicit condition on CFL number with respect to the discretization used. Moreover, it was shown that when one single particle lies in each cell of the grid mesh for one-dimensional problems, the optimal Courant number is allowed. This result highlights that the DGMPM stability condition is improved compared to that of the DGFEM. The enhanced stability combined with the use of the characteristic structure of the solution lead to a scheme that provides the ability to accurately follow waves in media undergoing finite deformations under total Lagrangian approach. Comparisons with an analytical one-dimensional solution in sections 7.2 and 7.3, and with the finite element method for two-dimensional problems in section 7.4 confirmed the good behavior of the method.
The approach can also be adapted to arbitrary high order approximations. However, well-suited polynomials com-  bined with reconstruction based quadrature such as Moving Least Squares or Spline interpolation will be required in order to avoid numerical instabilities due to the nonphysical negative nodal masses [41,42]. Moreover, the arbitrariness of the grid can be fully exploited by employing adaptive mesh techniques on the reference configuration in order to track accurately waves in the current one for problems involving complex geometries.

Appendix A. Analytical solution of the one-dimensional problem
We focus here on the analytical structure of the solution of the hyperbolic system written in the conservative form (5) for one-dimensional strain states considering a hyperelastic Saint-Venant-Kirchhoff material (SVK). In that case, the deformation gradient and the PK1 tensor are respectively: which corresponds to a plane wave solution. The stress tensor is given by the following SVK constitutive equations: where λ, µ are the Lamé's coefficients, Σ is the second Piola-Kirchhoff stress tensor and E the Green-Lagrange strain tensor defined as: These definitions for the one-dimensional strain problem lead to the stress components: Thus, only the axial stress will be considered since the other component can be deduced from it. For simplicity in what follows, the stress will read: The Jacobian matrix is then: with eigen-values and associated right eigen-vectors: The two characteristic fields of the system are genuinely nonlinear since the following condition is satisfied for each field: Hyperbolicity of the system is ensured by real eigen-values, leading to the following condition on the deformation gradient: strict hyperbolicity being ensured with the strict inequality. Moreover, the dependency on deformation gradient of characteristic speed can lead to two types of wave in the solution. Indeed, considering the Riemann problem: If U(X, 0) is such that F L > F R , then ω 1 L >ω 1 R and the first characteristic field corresponds to a centered rarefaction wave meaning that the characteristics move away from each other. Furthermore, ω 2 L >ω 1 R implies that the second characteristic fields is related to a shock wave in which the characteristics collide and generate a shock (see figure A.18). The same reasoning leads to an opposite configuration for the case F L < F R (see figure A. 19). The two structures of the solution (1-rarefaction,2-shock or 1-shock,2-rarefaction) are similar to the wave pattern that result from the Dam-Break problem of shallow water equations. Hence, the analytical solution of the hyperelastic bar impact problem for the SVK material is close to the dam-break solution developed in [8]. Two states U lying on each side of a shock wave are connected by the so-called Rankine-Hugoniot condition: where • denotes a jump across the i th discontinuous wave traveling at speed s i . The jump condition must be satisfied between a known state U and an unknown state on the other side of the shock U * : Substitution of the shock speed in the second equation with its expression given by the first one leads to a second order polynomial in v * : It is easy to verify that the Lax entropy condition which prevents from an expansion shock (ω(U L ) > s i >ω (U R )), ensures that the discriminant of this polynomial is positive and the solution v * then reads: 16) corresponding to the two families of shock. Finally, by assuming a weak shock (U ≃ U * ) one gets the conditions to be satisfied: Centered rarefaction waves A centered rarefaction wave is a self-similar solution defined in genuinely nonlinear fields such that: U(X, t) =Ũ(ξ = X/t) (A. 19) whereŨ(ξ = X/t) traces out a curve which tangent vector at each point is collinear to a given eigen-vector R i . Such a smooth curve is an integral curve defined byŨ ′ (ξ) = α(ξ)R i (Ũ(ξ))). With the parameter ξ = X/t defined before, one shows that the solution across a rarefaction wave must satisfy: where the denominator is nonzero in a genuinely nonlinear field due to condition (A.9). The integration between ξ = ω i (U L ) and ξ = ω i (U R ) of this system of ordinary differential equations, after some algebra, leads to the expression of integral curves. These expression can be used to compute a state U * connected to a known state U through a centered rarefaction wave: for a 1-rarefaction (A.21)