François Vilar
email: francois-vilar@brown.edu

Pierre-Henri Maire
email: maire@celia.u-bordeaux1.fr

Rémi Abgrall
email: remi.abgrall@math.u-bordeau1.fr

A total Lagrangian discontinuous Galerkin discretization of the two-dimensional gas dynamics equations

Keywords: Discontinuous Galerkin discretization, total Lagrangian formulation, updated Lagrangian formulation, cell-centered scheme, Godunov-type method, unstructured moving grid, curvilinear grid, gas dynamics PACS: 47.11.Df, 47.10.ab, 47.40.Nm 2000 MSC: 76N15, 65M06

Based on the total Lagrangian kinematical description, a discontinuous Galerkin (DG) discretization of the gas dynamics equations is developed for two-dimensional fluid flows on general unstructured grids. Contrary to the updated Lagrangian formulation, which refers to the current moving configuration of the flow, the total Lagrangian formulation refers to the reference fixed configuration, which is usually the initial one. In this framework, the Lagrangian and Eulerian descriptions of the kinematical and the physical variables are related by means of the Piola transformation. Here, we describe a cell-centered high-order DG discretization of the physical conservation laws. The geometrical conservation law, which governs the time evolution of the deformation gradient, is solved by means of a finite element discretization. This approach allows to satisfy exactly the Piola compatibility condition. Regarding the DG approach, it relies on the use of a polynomial space approximation which is spanned by a Taylor basis. The main advantage in using this type of basis relies on its adaptability regardless the shape of the cell. The numerical fluxes at the cell interfaces are computed employing a node-based solver which can be viewed as an approximate Riemann solver. We present numerical results to illustrate the robustness and the accuracy up to third-order of our DG method. First, we show its ability to accurately capture geometrical features of a flow region employing curvilinear grids. Second, we demonstrate the dramatic improvement in symmetry preservation for radial flows.

Introduction

We aim at describing a high-order discontinuous Galerkin (DG) method for solving the twodimensional total Lagrangian form of the gas dynamics equations on general unstructured grids.

It is well known that fluid dynamics relies on two kinematics descriptions: the Eulerian or spatial description and the Lagrangian or material description, refer for instance to [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF]18]. In the former, the conservation laws are written using a fixed reference frame whereas in the latter they are written through the use of a time dependent reference frame that follows the fluid motion. The Lagrangian representation is particularly well adapted to describe the time evolution of fluid flows contained in regions undergoing large shape changes due to strong compressions or expansions. Further, in this approach, there is no mass flux across the boundary surface of a control volume moving with the fluid velocity. Thus, Lagrangian representation provides a natural framework to track accurately material interfaces in multi-material compressible flows. This paper is primarily concerned with the development of a Lagrangian method whose main feature relies on the use of the total Lagrangian formalism. In this approach, the physical conservation laws are written employing the Lagrangian coordinates which refer to the initial configuration of the fluid flow. Moreover, in these equations the divergence and gradient operators are expressed by means of the Piola transformation [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF], which requires the knowledge of the deformation gradient tensor, i.e., the Jacobian matrix associated to the Lagrange-Euler flow map. The deformation gradient tensor characterizes the time evolving deformation and is governed by a partial differential equation named the geometric conservation law (GCL). To ensure the consistency between the initial and the current configurations , the deformation gradient tensor has to satisfy an involutive constraint [START_REF] Plohr | A Conservative Eulerian Formulation of the Equations for Elastic Flows[END_REF],which implies the Piola compatibility condition. The total Lagrangian approach is very well known in the solid mechanics community wherein it is extensively used to model solid dynamics undergoing large deformations [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF]. The first application of the total Lagrangian approach to the gas dynamics equations has been undertaken in [START_REF] Abgrall | A Lagrangian Discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems[END_REF][START_REF] Loubère | Une Méthode Particulaire Lagrangienne de type Galerkin Discontinu[END_REF] by means of a DG type discretization. However, the use of the aforementioned method is restricted to a representation on the initial configuration since it cannot be rigorously re-interpreted on the current configuration. We also note that the theoretical properties of the gas dynamics equations written under the total Lagrangian formulation have been thoroughly studied in [START_REF] Després | Lagrangian Gas Dynamics in Two Dimensions and Lagrangian systems[END_REF][START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique Lagrangienne bidimensionelle[END_REF]. On the contrary to the total Lagrangian formulation, the updated Lagrangian formulation is a moving domain method, which is widely employed. In this approach, the gas dynamics equations are written employing the Eulerian coordinates, which refer to the current configuration of the fluid flow. The time derivative of the physical variables is taken following the fluid particles paths: this the material derivative. The integral formulation of the conservation laws is readily obtained by employing the Reynolds transport formula over an arbitrary moving control volume. The time rate of change of a zone volume is governed by a partial differential equation which is the updated Lagrangian form of the GCL. It is worth mentioning that at the discrete level the zone volume computed from its vertices coordinates must rigorously coincide with the zone volume deduced from the numerical solution of the GCL. This critical requirement is the cornerstone on which any proper multi-dimensional updated Lagrangian scheme should rely. Two approaches are mainly employed to solve the updated Lagrangian formulation of the gas dynamics equations. The first one, which is called the staggered grid hydrodynamics, consists in using a staggered discretization wherein the kinematic variables (vertex position, velocity) are located at nodes whereas the thermodynamic variables (density, pressure, internal energy) are defined at the cell centers. The conversion of kinetic energy into internal energy through shock waves, consistently with the second law of thermodynamics, is ensured by adding an artificial viscosity term. The staggered grid schemes employed in most hydro-codes have been remarkably successful over the past decades in solving complex multi-dimensional compressible fluid flows, refer for instance to [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF][START_REF] Caramana | Elimination of artificial grid distorsion and hourglass-type motions by means of Lagrangian subzonal masses and pressures[END_REF]. However, they clearly have some theoretical and practical deficiencies such as mesh imprinting and symmetry breaking. In addition, the fact that all variables are not conserved over the same space can lead to serious difficulties in the perspective of an arbitrary Lagrangian-Eulerian (ALE) extension. The second approach, known as cell-centered hydrodynamics, employs a cell-centered placement of all hydrodynamic variables including the momentum. This approach consists of a moving mesh finite volume method wherein the numerical fluxes (multi-valued nodal pressures and nodal velocity) are computed through the use of node-centered approximate Riemann solvers. In this framework, momentum and total energy are conserved and an entropy inequality is satisfied at the semi-discrete level to ensure the thermodynamic consistency of the numerical method. Moreover, the numerical fluxes are constructed to satisfy the GCL compatibility. The interested readers may refer to the following papers [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme in arbitrary dimension[END_REF][START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF][START_REF] Cheng | A high order eno conservative lagrangian type scheme for the compressible euler equations[END_REF][START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered Lagrangian hydrodynamics on polygonal grids[END_REF][START_REF] Barlow | A high order cell centred dual grid Lagrangian Godunov scheme[END_REF][START_REF] Burton | A cell-centered Lagrangian Godunov-like method for solid dynamics[END_REF] for a more detailed description of this approach and its variants. Let us point out that work has been done to investigate the relationships between the staggered and the cell-centered discretizations, refer to [START_REF] Maire | Staggered Lagrangian discretization based on cellcentered Riemann solver and associated hydrodynamics scheme[END_REF][START_REF] Loubère | 3D staggered Lagrangian hydrodynamics scheme with cell-centered Riemann solver-based artificial viscosity[END_REF]. It is worth mentioning that the above approaches belong to the class of either finite volume or finite difference based on a piecewise constant representation of the hydrodynamic variables and are thus restricted to a second-order accuracy. Up to our knowledge, the interpretation of the staggered schemes of Goad [START_REF] Goad | WAT: A Numerical Method for Two-Dimensional Unsteady Fluid Flow[END_REF] and Wilkins [START_REF] Wilkins | Methods in Computational Physics, volume 3, chapter Calculation of elasticplastic flows[END_REF] by means of a finite element method has been initially introduced by Lascaux at the beginning of the 70's [START_REF] Lascaux | Application de la méthode des éléments finis en hydrodynamique bi-dimensionnelle utilisant les variables de Lagrange[END_REF][START_REF] Lascaux | Application of the Finite Element Method to 2D Lagrangian hydrodynamics[END_REF]. This finite element approach has been further developed producing various interesting staggered schemes. For instance, a compatible finite element Lagrangian hydrodynamics algorithm used in a multi-material ALE strategy has been described in [START_REF] Barlow | A compatible finite element multi-material ALE hydrodynamics algorithm[END_REF]. We also note the development of a variational multiscale stabilized approach in finite element computation of Lagrangian hydrodynamics where a piecewise linear approximation was adopted for the variables [START_REF] Scovazzi | Stabilized shock hydrodynamics: I. A Lagrangian method[END_REF][START_REF] Scovazzi | Stabilized shock hydrodynamics: II. Design and physical interpretation of the SUPG operator for Lagrangian computations[END_REF]. The case of Q1/P0 finite element is studied in [START_REF] Scovazzi | Multi-scale Lagrangian shock hydrodynamics on Q1/P0 finite elements: Theoretical framework and two-dimensional computations[END_REF] where the kinematic variables are represented using a piecewise linear continuous approximation while the thermodynamic variables utilize a piecewise constant representation. Except the pioneering work of [START_REF] Abgrall | A Lagrangian Discontinuous Galerkin-type method on unstructured meshes to solve hydrodynamics problems[END_REF][START_REF] Loubère | Une Méthode Particulaire Lagrangienne de type Galerkin Discontinu[END_REF], all the aforementioned approaches are characterized by an accuracy which is at most of second order. This accuracy restriction is a natural consequence of the spatial discretization of the Lagrange-Euler flow map employed. Namely, the gas dynamics equations are discretized on a moving grid made of polygonal cells whose edges remain straight lines throughout the motion. This amounts to claim that the Lagrange-Euler flow map admits a linear continuous representation with respect to Eulerian coordinates over the deforming computational grid. Further, the kinematic velocity field also admits a linear continuous representation. Therefore, as noticed in [START_REF] Cheng | A third-order conservative Lagrangian type scheme on curvilinear meshes for the compressible Euler equations[END_REF], this approximation of the grid motion implies a second-order error in the numerical method. To reach a higher order of accuracy, one has to take into account a higher order discretization of the kinematics of the flow. This task has been successfully undertaken in a series of papers [START_REF] Dobrev | Curvilinear finite elements for Lagrangian hydrodynamics[END_REF][START_REF] Dobrev | High Order Curvilinear Finite Elements for Lagrangian Hydrodynamics[END_REF][START_REF] Dobrev | High Order Curvilinear Finite Elements for axisymmetric Lagrangian Hydrodynamics[END_REF]. In [START_REF] Dobrev | High Order Curvilinear Finite Elements for Lagrangian Hydrodynamics[END_REF], the authors describe a high-order finite element framework for solving the gas dynamics equations on curvilinear moving grids. Their method relies on the introduction of a continuous high-order representation of the flow kinematics. The spatial discretization is obtained by means of a variational formulation of the gas dynamics equations. The kinematic variables are expanded over a basis of continuous high-order polynomial functions, whereas the thermodynamic variables are expanded over a basis of discontinuous high-order polynomial functions. The dissipation of kinetic energy into internal energy is ensured by the adding of a high-order tensorial artificial viscosity [START_REF] Tz | A tensor artificial viscosity using a finite element approach[END_REF]. This approach, which can be viewed as the high-order extension of the staggered grid hydrodynamics produces very impressive numerical results. Let us point out that this method is able to more accurately capture geometrical features of a flow region, to maintain robustness with respect to mesh motion and to improve the symmetry preservation in symmetric flows. More recently, an isogeometric analysis of Lagrangian shock hydrodynamics has been proposed in [START_REF] Bazilevs | Isogeometric analysis of lagrangian hydrodynamics[END_REF]. In this approach, the Euler equations of compressible hydrodynamics in the weak form are discretized using NURBS (Non Uniform Rational B-splines) in space. This discretization provides another high-order extension of the staggered grid hydrodynamics which also requires the adding of an artificial viscosity. It has the main advantage of being characterized by an exact symmetry preservation. All these very promising results have motivated us to present another contribution to the new domain of high-order numerical methods for Lagrangian hydrodynamics. Here, we present a cellcentered DG high-order discretization of the total Lagrangian formulation of the gas dynamics equations. The GCL, which governs the time rate of change of the deformation gradient, is discretized by means of a finite element approximation. This allows to satisfy exactly the Piola compatibility condition. The knowledge of the deformation gradient allows to represent consistently the curvilinear grid on the actual configuration. The DG discretization employs a cell-based expansion of the physical variables (including the momentum) over a Taylor basis, i.e., the polynomial terms correspond to those obtained when performing a Taylor expansion at the cell centroid. The interest of this type of basis, which has been already used in [START_REF] Luo | A DG method based on a Taylor basis for the compressible flows on arbitrary grids[END_REF], relies on its adaptability regardless the shape of the cell. The numerical fluxes at the cell interfaces are computed by means of the node-centered solver which has been introduced in [START_REF] Maire | A cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF]. The numerical algorithm not only satisfies the GCL compatibility condition but also conserves momentum and total energy. Further, it satisfies a local entropy inequality at the semi-discrete level. The time discretization employs the classical third-order TVD (Total Variation Diminishing) Runge-Kutta method [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. It is worth pointing out that in the case of a piecewise constant polynomial basis, our DG method boils down to the classical first-order finite volume cell-centered scheme developed in [START_REF] Maire | A cell-centered Lagrangian scheme for twodimensional compressible flow problems[END_REF]. Let us mention, that following the approach of [START_REF] Loubère | Une Méthode Particulaire Lagrangienne de type Galerkin Discontinu[END_REF], a high-order DG spectral finite element method has been introduced in [START_REF] Jia | A new high-order discontinuous Galerkin spectral finite element for Lagrangian gas dynamics in two-dimensions[END_REF]. However, this numerical method does not ensure properly the Piola compatibility condition required to solve the GCL. We would like to mention the recent works [START_REF] Boutin | Extension of ALE methodology to unstructured conical meshes[END_REF][START_REF] Pino | A curvilinear finite-volume method to solve compressible gas dynamics in semi-Lagrangian coordinates[END_REF] which describe respectively a curvilinear finite volume method and curvilinear finite element method for solving gas dynamics equations. The remainder of this paper is organized as follows: In Section 2, we recall the multi-dimensional gas dynamics equations written under the total Lagrangian formulation and their main properties. In Section 3, we present not only the DG discretization of the total Lagrangian formulation of the gas dynamics but also the finite element approximation of the deformation gradient. Finally, Section 4 provides a series of numerical experiments which assess the robustness and the accuracy of our DG method up to third-order. The numerical results obtained demonstrate the ability of our method to dramatically improve the symmetry preservation of symmetric flows.

Multi-dimensional Lagrangian and Eulerian descriptions for fluid flows

In this section, we recall the multi-dimensional gas dynamics equations written under total Lagrangian form. This system is a bit unusual. We show its connections to the Eulerian form and more importantly, what are the ingredients to achieve a full equivalence. Let us point out that a similar work has been done in [START_REF] Plohr | A Conservative Eulerian Formulation of the Equations for Elastic Flows[END_REF] for the equations of elastic flows. Indeed this is not fully straightforward because the two systems are obtained from one an other via a change of variable.

Since we expect discontinuous solutions, the main questions is whether or not we recover the same weak solutions. In this section we conduct the analysis and show a complete equivalence, provided some geometrical results are satisfied. The results obtained also give some hints on how to conduct the numerical approximations of the gas dynamics system written under the total Lagrange form.

Kinematics

To describe the main notions of kinematics that will be useful to describe fluid motions, we follow the presentation given in [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF]. Let X be the position vector of a point of the fluid in its initial configuration. The evolution of the fluid is then characterized by a time-dependent motion Φ, named the mapping, which defines the location of the point X at time t > 0, such as

x = Φ(X, t). (1)
At this point, we can introduce the two usual descriptions of the flows, namely the Lagrangian description and the Eulerian description. The Lagrangian description, otherwise called material description, consists in observing the fluid by following the motion of fluid particles from their initial location. The independent variables used for this description are (X, t). On the other hand the Eulerian description, otherwise called spatial description, consists in observing the fluid at fixed locations in the space. The independent variables used for this description are (x, t). Similarly, in the remainder all the quantities expressed using Lagrangian coordinates (X, t) will be defined by capital letters, while the quantities using the Eulerian coordinates (x, t) will be defined through lower case letters. We notice that through the use of the transformation (1) any fluid quantity f which is expressed in terms of Eulerian variables can also be expressed in terms of Lagrangian variables, and conversely. To emphasize the used variables and for the sake of conciseness, the same notation is used to denote the value of the physical quantity regardless the employed description

f = f (x, t) = f (Φ(X, t), t) = f (X, t). (2)
To be more precise, f (X, t) is the value of the physical quantity experienced at time t by the fluid particle initially located at X, whereas f (x, t) is the value of f experienced by the fluid particle which is located at position x at time t. Obviously, if the particle located at x at time t was initially located at X, both definitions of f express the same quantity, hence equation [START_REF] Barlow | A compatible finite element multi-material ALE hydrodynamics algorithm[END_REF]. Now, let f be a fluid variable with a sufficient smoothness to allow the computation of its first partial derivatives with respect to both Lagrangian and Eulerian variables. First, we introduce the material derivative of f which measures the rate of change of f following a fluid particle along its motion as

d dt f (x, t) ≡ ∂ ∂t f (X, t). (3)
Note that it corresponds to the partial time derivative in the Lagrangian description. The velocity of a fluid particle is denoted U and is defined as

U (X, t) = ∂ ∂t Φ(X, t). (4)
As defined, U is a function of the Lagrangian variables. However it is possible to also express it in terms of the Eulerian variables using definitions (1), ( 2) and (3). It is thus possible to determine the fluid flow by solving the system of ordinary differential equations

d x dt = U (x, t), x(X, 0) = X, (5)
which corresponds to the trajectory equations of the fluid flow. Considering the physical quantity f expressed in terms of the Eulerian variables, that is f = f (x, t), we compute its material derivative employing the chain rule of composite derivative and the trajectory equation ( 5) to get

d dt f (x, t) = ∂ ∂t f (x, t) + U (x, t) • ∇ x f (x, t), (6)
where ∇ x denotes the gradient operator with respect to Eulerian coordinates, while ∇ X denotes the gradient operator with respect to Lagrangian coordinates. This last equation may be interpreted as expressing the time rate of change of an arbitrary physical quantity f = f (x, t) apparent to an observer located on the moving particle instantaneously at the position x. By definition, Φ(X, 0) = X, further for t fixed, equation ( 1) characterizes the deformation of the fluid at time t. Namely, through the time-dependent map, neighboring points in the initial configuration X and X + dX transform through the flow motion into x and x + dx with dx = FdX.

The second-order tensor F is called the deformation gradient tensor. It is nothing but the Jacobian matrix associated to the flow map Φ

F = ∇ X Φ. (8)
Knowing that F(X, 0) = I d , we assume that for all t > 0, the determinant of F satisfies J = det F > 0 so that the flow map is invertible. Let Ω be a region occupied by the fluid in its initial configuration. Its image by the time-dependent map Φ at time t is denoted ω. It is such that

ω = {x | x = Φ(X, t), X ∈ Ω} .
Now, we briefly recall two fundamental relationships for transforming area and volume elements from the initial configuration to the actual one, and conversely. Let dX 1 and dX 2 be two linearly independent line elements in Ω. We define the area element dS as

dS = dX 1 × dX 2 . (9)
We define by dS the algebraic measure of dS, i.e., dS = |dS|. If N denotes the unit outward normal to the surface element, we have N dS = dX 1 × dX 2 . Let dx 1 and dx 2 be the images of line elements dX 1 and dX 2 by the flow map, they are defined by dx i = FdX i , for i = 1, 2. Since, J > 0, the Eulerian line elements are linearly independent and we can define the corresponding surface element as

nds = dx 1 × dx 2 .
Substituting dx i = FdX i in the above equation leads to

nds = FdX 1 × FdX 2 .
At this point, it is interesting to recall that the co-factor of F, i.e., F ⋆ = det(F)F -t , satisfies

F ⋆ (V 1 × V 2 ) = FV 1 × FV 2 ,
for all linearly independent vectors (V 1 , V 2 ), see [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF]. Employing this result, the Eulerian surface element rewrites

nds = F ⋆ (dX 1 × dX 2 ).
Finally, using (9), we express the Eulerian surface element in terms of the Lagrangian surface element through the Nanson's formula

nds = F ⋆ N dS. ( 10
)
Now, let us consider a third line element dX 3 such that the set {dX 1 , dX 2 , dX 3 } is a basis with a positive orientation. We define the Lagrangian volume element as

dV = (dX 1 × dX 2 ) • dX 3 .
Introducing dx 3 = FdX 3 leads to the following definition of the Eulerian volume element dv = (dx 1 × dx 2 ) • dx 3 . Expressing the Eulerian line element in terms of their Lagrangian counterparts yields

dv =(FdX 1 × FdX 2 ) • FdX 3 , =F t F ⋆ (dX 1 × dX 2 ) • dX 3 , = det(F)(dX 1 × dX 2 ) • dX 3 .
Hence, we obtain the classical formula which expresses the transformation of a volume element through the flow map dv = det(F)dV. (11)

Piola identities

In this paragraph, we briefly recall the Piola identities which are fundamental to relate the Eulerian form and the Lagrangian form of the gas dynamics equations. To this end, let us introduce some notations. In the following, the " " symbol is used to define the divergence operator in Lagrangian or Eulerian coordinates. Consequently, if V is a vector, the divergence of V , namely ∇ V , is a scalar such as ∇ V = tr (∇V ). In the case of a tensor T, the divergence of T, namely ∇ T, is the vector defined such that for any fixed vector A the following relation holds

(∇ T) • A = ∇ T t A .
Note that in this context, (.) t represents the transpose operator, i.e., if A is a matrix, A t is its transpose and A -t is a short hand notation for A -1 t = A t -1 . Now, let Ω be a fluid region in the initial configuration and ω = Φ(Ω, t) its image by the flow map at time t > 0. Applying the divergence theorem to ∇ X I d = 0 and

∇ x I d = 0 yields ∂Ω N dS = 0, ( 12a
) ∂ω n ds = 0. (12b)
These formulas correspond to the integral form of the Piola identities which have a simple geometric interpretation. Namely, the summation of the unit normal over a closed surface is equal to zero. Substituting the Nanson's formula [START_REF] Cheng | A high order eno conservative lagrangian type scheme for the compressible euler equations[END_REF] into (12a) and (12b) yields

∂ω (F ⋆ ) -1 n ds = 0, ∂Ω F ⋆ N dS = 0.
Then, applying the divergence theorem, we get

ω ∇ x ( 1 det(F) F t ) dv = 0, Ω ∇ X F ⋆ dV = 0.
Knowing that these formula hold for any arbitrary fluid regions, we finally obtain the Piola identities written under local form

∇ x ( 1 det(F) F t ) = 0, (13a)
∇ X F ⋆ = 0. (13b)
In what follows, we shall derive the gas dynamics equations written under the Lagrangian form and thus we shall utilize extensively the second Piola identity (13b).

Transformation formulas for the divergence and gradient operators

Let T be an arbitrary second-order tensor being a smooth function with respect to the Lagrangian coordinates (and hence the Eulerian ones via the mapping Φ). We shall express the divergence of this tensor with respect to the Eulerian coordinates in terms of the divergence to the Lagrangian coordinates. To this end, we first recall the divergence theorem over the Eulerian region

ω ω ∇ x T dv = ∂ω Tn ds. ( 14
)
Transforming the right hand-side by means of the Nanson's formula [START_REF] Cheng | A high order eno conservative lagrangian type scheme for the compressible euler equations[END_REF]

leads to ω ∇ x T dv = ∂Ω TF ⋆ N dS.
Then, applying the divergence theorem to the right hand-side, we get ω

∇ x T dv = Ω ∇ X (TF ⋆ ) dV.
Rewriting the left hand-side in terms of the Lagrangian coordinates and using [START_REF] Cheng | A third-order conservative Lagrangian type scheme on curvilinear meshes for the compressible Euler equations[END_REF] yields

Ω ∇ x T det(F) dV = Ω ∇ X (TF ⋆ ) dV.
Since the above equation holds for any arbitrary region Ω, we finally obtain the formula expressing the transformation of the divergence operator through the flow map

∇ x T = 1 det F ∇ X (TF ⋆ ). ( 15
)
For all U and V , arbitrary vectors which are smooth functions with respect to both Lagrangian and Eulerian coordinates, the following tensorial identity holds

∇ (U ⊗ V ) = (∇U )V + U (∇ V ), (16)
where the symbol "⊗" denotes the tensor product of vectors defined by (U ⊗ V ) ij = U i V j , for all i, j = 1, . . . , d. Using the above identity and substituting in equation ( 15) the tensor T = A ⊗ V , where A is a fixed vector, yields

A(∇ x V ) = 1 det F ∇ X ((A ⊗ V )F ⋆ ), = 1 det F ∇ X (A ⊗ (F ⋆ t V )), = 1 det F A(∇ X (F ⋆ t V )).
And since this relation holds for any constant vector A, it follows that

∇ x V = 1 det F ∇ X (F ⋆ t V ), (17)
which expresses the relation between the divergence of a vector function written using both Lagrangian and Eulerian coordinates. Finally, substituting T = ϕ I d in [START_REF] Dobrev | Curvilinear finite elements for Lagrangian hydrodynamics[END_REF], where ϕ is a smooth scalar function with respect to both Lagrangian and Eulerian coordinates, and using the following tensorial identity

∇ (ϕ T) = ϕ(∇ T) + T∇ϕ, (18)
one easily obtains the following formula relating to the transformation of the gradient operator

∇ x ϕ = 1 det F ∇ X (ϕ F ⋆ ). ( 19
)
We have now expressed the relations between the gradient and divergence operators written in terms of Eulerian coordinates and its counterpart written in terms of Lagrangian coordinates. Using the Piola identity, this last relation also reads

∇ x ϕ = 1 det F F ⋆ ∇ X ϕ = F -t ∇ X ϕ. (20)

Geometrical conservation law

We have seen in the previous paragraphs that the expressions of the differential operators in terms of the Lagrangian coordinates ( 15), ( 17) and ( 19) require the knowledge of the deformation gradient F.

Here, we derive the conservation law which governs the time evolution equation of the deformation gradient. Using the trajectory equation ( 4), the deformation gradient tensor definition F = ∇ X Φ yields an equation on the time rate of change of F(X, t)

∂ F ∂t -∇ X U = 0, (21)
supplemented with the initial condition F(X, 0) = I d . It is important to note that the solution of (21) defines a deformation gradient which derives from a motion Φ provided that F satisfies the compatibility condition

∇ X × F = 0, (22)
where the tensor ∇ X × F is the curl of the deformation gradient defined for all constant vector A by

(∇ X × F)A = ∇ X × (F t A).
The compatibility condition ( 22) is an involutive constraint, i.e., if the compatibility condition is satisfied at time t = 0 then it is satisfied for all time t > 0. To show this result, we compute the time rate of change of the curl of F knowing that the deformation gradient tensor satisfies the conservation law [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF], as

∂ ∂t (∇ X × F) = ∇ X × ∂ F ∂t = ∇ X × ∇ X U .
By definition of the curl operator, for all constant vector

A (∇ X × ∇ X U ) A = ∇ X × (∇ X U ) t A . Knowing that (∇ X U ) t A = ∇ X (A • U ), we readily obtain that ∇ X × (∇ X U ) t A = ∇ X × (∇ X (A • U )) = 0.
Hence, ∂ ∂t (∇ X × F) = 0, which writes equivalently

(∇ X × F)(X, t) = (∇ X × F)(X, 0), for all t > 0. ( 23
)
Next, we show that if the involutive constraint [START_REF] Jia | A new high-order discontinuous Galerkin spectral finite element for Lagrangian gas dynamics in two-dimensions[END_REF] is satisfied then the Piola identity (13b) holds.

To this end, let A and B be two arbitrary constant vectors. Recalling that the cofactor of F is characterized by

F ⋆ (A × B) = FA × FB,
and knowing that

(A × B) • (∇ X F ⋆ ) = ∇ X F ⋆ t (A × B) ,
we readily obtain

(A × B) • (∇ X F ⋆ ) =∇ X (F t A × F t B), = F t B • ∇ X × (F t A) -F t A • ∇ X × (F t B) , = F t B • ((∇ X × F)A) -F t A • ((∇ X × F)B) .
Hence, if F is such that ∇ X × F = 0, then for all A and B, (A × B) • ∇ X F ⋆ = 0, thus ∇ X F ⋆ = 0 and the Piola identity (13b) is satisfied. Now, denoting by J the determinant of the deformation gradient, i.e., J = det F, and knowing that F satisfies the conservation law [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF], we compute the time rate of change of J. Let us point out that J represents the ratio of the Eulerian volume element to the Lagrangian volume element, i.e. J = dv dV . Further, if ϕ is a scalar function expressed in terms of F then the chain rule of composed derivative reads

∂ ϕ ∂t = ∂ϕ ∂F : ∂ F ∂t , (24)
where the symbol ":" denotes the inner product of second-order tensors defined by S : T = tr(S t T).

In [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cellcentered Lagrangian scheme[END_REF], ∂ϕ ∂F is the second-order tensor whose components read ∂ϕ ∂F ij , where F ij denotes the generic component of F. Applying the above chain rule for ϕ(F) = det F leads to

∂ J ∂t = F ⋆ : ∂ F ∂t . (25)
Here, we have used the fact that the derivative of the determinant of F with respect to its argument is equal to the cofactor of F, see [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF], that is

∂ ∂F (det F) = F ⋆ .
By virtue of [START_REF] Tz | A tensor artificial viscosity using a finite element approach[END_REF], taking the inner product of the conservation law ( 21) by the cofactor F ⋆ yields the time rate of change of the Jacobian

∂ J ∂t -F ⋆ : ∇ X U = 0, (26)
which is written under a non-conservative form. However, thanks to the tensorial identity

∇ (T t V ) = V • (∇ T) + T : ∇V , (27)
and also to the Piola identity (13b), equation ( 26) is rewritten under the conservative form as

∂ J ∂t -∇ X (F ⋆ t U ) = 0. ( 28
)
Integrating the above equation over the Lagrangian region Ω and applying the divergence theorem leads to

d dt Ω J dV - ∂Ω U • F ⋆ N dS = 0.
Introducing ω = Φ(Ω, t) and using the Nanson's formula [START_REF] Cheng | A high order eno conservative lagrangian type scheme for the compressible euler equations[END_REF] and [START_REF] Cheng | A third-order conservative Lagrangian type scheme on curvilinear meshes for the compressible Euler equations[END_REF], allows to rewrite the above equation in the current configuration

d dt ω dv - ∂ω U n ds = 0. ( 29
)
This last equation is nothing but the time rate of change of the volume of the Eulerian region ω, that is why it is called the geometric conservation law (GCL). Equation [START_REF] Lascaux | Application of the Finite Element Method to 2D Lagrangian hydrodynamics[END_REF] expresses the time rate of change of the Jacobian expressed in terms of the Lagrangian coordinates. To define the counterpart equation in the actual configuration, we make use of the material derivative definition (3) and divergence relation [START_REF] Dobrev | High Order Curvilinear Finite Elements for axisymmetric Lagrangian Hydrodynamics[END_REF], and finally get

d J dt -J ∇ x U = 0, (30)
where all the quantities involved are defined on the actual configuration, using (x, t) as coordinates.

Physical conservation laws

Here, we recall the formulation of the conservation laws of mass, momentum and total energy expressed with respect to the initial and the actual configurations. For any initial configuration Ω, ω(t) = Φ(Ω, t) is its image by the flow map at time t > 0.

Conservation of mass

Mass conservation is expressed by

d dt ω(t) ρ(x, t) dv = 0.
Since ω(0) = Ω, rewriting this last relation in terms of the Lagrangian coordinates yields

d dt Ω ρ(X, t) J(X, t) dv = 0. (31)
Since Ω and thus ω(t) is arbitrary, and denoting by ρ 0 the initial density field, equation [START_REF] Luo | A DG method based on a Taylor basis for the compressible flows on arbitrary grids[END_REF] leads to

∂ ∂t (ρ J) = 0, (32a)
ρ(X, t) J(X, t) = ρ 0 (X), for all t > 0,

recalling that accordingly to relation (2), ρ(X, t) (resp. J(X, t)) expresses the density (resp. the Jacobian) of a particle at time t, initially located at X, i.e., ρ(X, t) = ρ(x(X, t), t). In equations [START_REF] Maire | A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes[END_REF], we have derived two form of the mass conservation equation written under total Lagrangian form.

The use of identities (32b) and equation ( 28) yields an equation on the time rate of change of the specific volume, as

ρ 0 ∂ ∂t ( 1 ρ ) -∇ X (F ⋆ t U ) = 0. ( 33
)
By means of relation (32b) and of the transformation formulas presented, the counterpart of this last equation expressed in terms of the Eulerian coordinates naturally writes

ρ d dt ( 1 ρ ) -∇ x U = 0. ( 34
)

Conservation of momentum

If t represents a force defined per unit area acting on the boundary surface ∂ω(t), one knows that t = Tn, where n is the local unit outward normal and T the second-order tensor named the Cauchy stress tensor, see [18] for more details. Hence, the Newton's law applied to ω(t), in the absence of volumic forces, writes

d dt ω(t) ρU dv = ∂ω(t)
Tn ds.

This can be transformed into volume integral as

d dt ω(t) ρU dv = ω(t) ∇ x T dv.
Using Reynolds transport theorem [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF] and a change of variables, we get

d dt ω(t) ρU dv = ω(t) ρ d U dt dv, = Ω ρ ∂ U ∂t J dV, = Ω ρ 0 ∂ U ∂t dV. thanks to (32b)
Finally, using [START_REF] Dobrev | Curvilinear finite elements for Lagrangian hydrodynamics[END_REF], we obtain

Ω ρ 0 ∂ U ∂t dV - Ω ∇ X (TF ⋆ ) dV = 0.
Since this relation is valid for any volume Ω, we obtain the momentum equation

ρ 0 ∂ U ∂t -∇ X (TF ⋆ ) = 0. ( 35
)
Its integral version in the reference element is

d dt Ω ρ 0 U dV - ∂Ω TF ⋆ N dV = 0.
Using the same technique employed previously, the momentum equation ( 35) written using Eulerian coordinates (x, t) reads

ρ d U dt -∇ x T = 0. (36)

Conservation of total energy

We consider a fluid which does not conduct heat and without volumic source of heat. The conservation of total energy E, i.e., the sum of the specific internal energy ε and the specific kinetic

energy 1 2 U 2 , writes d dt ω(t) ρE dv = ∂ω(t) t • U ds, = ∂ω(t) (TU ) • n ds, since T is symmetric = ω(t) ∇ x (TU ) dv.
The symmetry of matrix T is a classical consequence of the conservation of the angular momentum, refer to [START_REF] Gurtin | The Mechanics and Thermodynamics of Continua[END_REF]. Using [START_REF] Dobrev | High Order Curvilinear Finite Elements for axisymmetric Lagrangian Hydrodynamics[END_REF] and the same arguments as for the momentum, we arrive at

Ω ρ 0 ∂ E ∂t dV - Ω ∇ X F ⋆ t TU dV = 0.
Finally, the local form of the total energy equation in the reference configuration writes

ρ 0 ∂ E ∂t -∇ X F ⋆ t TU = 0, (37)
while its integral form is

d dt Ω ρ 0 E dV - ∂Ω TU • F ⋆ N dV = 0.
Total energy equation [START_REF] Noh | Errors for calculations of strong shocks using artificial viscosity and an artificial heat flux[END_REF] expressed in terms of Eulerian coordinates writes

ρ d E dt -∇ x (TU ) = 0. ( 38
)

Gas dynamic equations written under total Lagrangian form

The fluid under consideration is inviscid and then the Cauchy stress tensor reduces to T = -P I d , where P > 0 is the thermodynamic pressure. In that case, the local form of the physical equations is (recall

F ⋆ = JF -t ) ∂ F ∂t -∇ X U = 0, (39a)
ρ 0 ∂ U ∂t + ∇ X (P F ⋆ ) = 0, (39b)
ρ 0 ∂ E ∂t + ∇ X (F ⋆,t P U ) = 0. ( 39c
)
The thermodynamic closure is given by the equation of state P = P (ρ, ε) where ε = E -1 2 U 2 denotes the specific internal energy. It is worth mentioning that the deformation gradient tensor equation (39a) implies the specific volume conservation equation [START_REF] Maire | A high-order one-step sub-cell force-based discretization for cell-centered Lagrangian hydrodynamics on polygonal grids[END_REF], using the mass conservation (32b) and the equation on the Jacobian (28). In the above system, the mass density is obtained by means of the mass conservation written under total Lagrangian form as ρ det F = ρ 0 .

Entropy consideration

Assuming enough smoothness, we derive the time rate of change of the specific entropy η. We first recall the Gibbs identity

θdη = dε + P d( 1 ρ ),
where θ > 0 is the absolute temperature. Since ε = E -1 2 U 2 , and 1 ρ = J ρ 0 , we can express the time rate of change of the specific entropy as

θ ∂ η ∂t = ∂ E ∂t -U • ∂ U ∂t + P ρ 0 ∂ J ∂t . Since ∂ J ∂t = F ⋆ : ∂ F ∂t , we finally get ρ 0 θ ∂ η ∂t = ρ 0 ∂ E ∂t -ρ 0 U • ∂ U ∂t + P F ⋆ : ∂ F ∂t . (40)
Replacing the time rate of change of F, U , and E, and using the tensorial identity [START_REF] Lascaux | Application de la méthode des éléments finis en hydrodynamique bi-dimensionnelle utilisant les variables de Lagrange[END_REF], one obtains

ρ 0 θ ∂ η ∂t = -∇ X (F ⋆ t P U ) + U • (∇ X (F ⋆ P )) + P F ⋆ : ∇ X U = 0.
Let us point out that this results has been obtained without recourse to Piola identity. Recalling that ∂ ∂t η(X, t) = d dt η(x, t), we conclude that for smooth flows the specific entropy is conserved along the trajectory equation.

Jump relations at a surface of discontinuity

Though the derivation of the jump relations for systems of conservation law is very classical, the system (39) is not very common, so we prefer to derive the jump relations from scratch. Moreover, the results obtained are at the core of the discretization method. Hence we believe it is useful to know exactly where they come from. In particular, we show that the jump relation obtained on F as well as those obtained from the Piola compatibility condition (13b) play a central role to connect the jump relations in Eulerian and Lagrangian coordinates. We can therefore anticipate that the discretization of theses equations will be also at the center of our method. The system [START_REF] Plohr | A Conservative Eulerian Formulation of the Equations for Elastic Flows[END_REF] is in conservation form and may admit discontinuous solutions. We use the notations of figure 1, in which Σ(t) is a discontinuity moving in the initial domain Ω, W Σ being its local velocity (W Σ = W Σ • N Σ being the normal velocity). As usual, • represents the jump operator defined, for any locally smooth term f (scalar, vector or tensor), and any point X ∈ Σ, by

Ω + Ω - Σ(t) N Σ
f = f + -f -, where f ± = lim h→0 ± f (X + hN Σ ) .
The Rankine Hugoniot relations are obtained from a classical analysis, note however that relation (39a) is a little unusual because being a relation between matrices.

Geometrical conservation law: the Hadamard relations

Each subdomain Ω ± being fixed except for the moving boundary Σ(t), the Reynolds transport formula applied on the deformation gradient tensor F yields

d dt Ω ± F dV = Ω ± ∂ F ∂t dV ∓ ∂Ω ± ∩Σ(t) W Σ F dS,
where W Σ is the (local) normal velocity of the surface Σ(t). Now, adding together the two terms corresponding respectively to Ω + and Ω -leads to

d dt Ω F dV = Ω ∂ F ∂t dV - Σ(t) W Σ F dS. ( 41
)
Similarly, using the divergence theorem, one gets

Ω ± ∇ X U dV = ∂Ω ± \Σ(t) U ⊗ N dS ∓ Σ(t) U ⊗ N Σ dS.
Again, adding together the contributions of Ω + and Ω -, it easily follows that

∂Ω U ⊗ N dS = Ω ∇ X U dV + Σ(t) U ⊗ N Σ dS. ( 42
)
In the end, subtracting relation ( 42) from ( 41), one obtains

d dt Ω F dV - ∂Ω U ⊗ N dS = Ω ∂ F ∂t -∇ X U dV - Σ (W Σ F + U ⊗ N Σ ) dS.
Using equation (39a) as well as its integral form, we finally obtain the Hadamard relation

W Σ F + U ⊗ N Σ = 0. (43)
One can apply this relation to the vectors N Σ and T Σ , which is tangent to Σ (and thus orthogonal to N Σ ). Doing so, one gets

W Σ F N Σ + U = 0, (44a)
and

W Σ F T Σ = 0. (44b)
Let us further discuss these relations. In the case where

W Σ = 0, from (44a) it is clear that U is continuous across Σ. Now, if W Σ = 0, one naturally finds that F N Σ = -U W Σ .
Furthermore, for any vector T orthogonal to N Σ , the following relation holds

F T = 0, (45)
which means that FT has no jump provided that T is normal to N Σ . Consequently, one can finally states that

F W Σ + U = 0, (46)
which is another form of the Hadamard compatibility condition. This last relation can also be derived from the definition of the discontinuity velocity w σ in the Eulerian frame, where σ is nothing but the image of the discontinuity Σ through the fluid flow. Thus, let X Σ (t) be the vector position of a point attached to the discontinuity surface Σ during its motion in the Lagrangian frame. Let x σ (t) be the image of X Σ (t) through the flow map Φ, i.e., x σ = Φ(X Σ (t), t). Knowing that the Eulerian discontinuity σ is itself the image of the Lagrangian discontinuity Σ in the flow map yields x σ ∈ σ. Bearing this in mind, the respective speeds of the Lagrangian and Eulerian discontinuities read as

W Σ = d X Σ dt , w σ = d x σ dt . (47)
Now, recalling that x σ = Φ(X Σ (t), t) and employing the chain rule leads to the two following relations, respectively in Ω + and Ω -

d x σ dt = ∇ X Φ ± d X Σ dt + ∂ Φ ± ∂t ,
which, by means of the flow velocity and deformation gradient tensor definitions, immediately yields the following relations on the discontinuity Eulerian and Lagrangian velocities

w σ = F ± W Σ + U ± . (48)
One can see these last relations are perfectly consistent with the previous result stated in [START_REF] Vilar | Cell-centered discontinuous Galerkin discretizations for two-dimensional scalar conservation laws on unstructured grids and for one-dimensional Lagrangian hydrodynamics[END_REF]. Now, let us show that relation (44b) has a geometrical interpretation related to the second Piola identity.

Considering T 1 and

T 2 such that T 1 × T 2 = N Σ and recalling that F ⋆ (T 1 × T 2 ) = FT 1 × FT 2 , it follows that F ⋆ N Σ = FT 1 × FT 2 .
Applying the jump operator to this inequality, and denoting by T = 1 2 T + + T -, leads to

F ⋆ N Σ = F T 1 × F T 2 + F T 1 × F T 2 = 0,
since F T = 0 for any vector orthogonal to N Σ . Hence, in the case of a non-steady discontinuity, the Hadamard compatibility relation implies the continuity of the normal component of the cofactor of F. This is also a consequence of the second Piola condition (13b), ∇ X F ⋆ = 0. Indeed, using the same arguments as those used to get [START_REF] Sod | A survey of several finite difference methods for systems of non-linear hyperbolic conservation laws[END_REF], one gets

Ω ∇ X F ⋆ dV + Σ F ⋆ N Σ dS = 0,
which leads to the following fundamental result

F ⋆ N Σ = 0. ( 49
)
This relation has been obtained without assuming anything more on Σ, and expresses the continuity of the normal in the Eulerian frame.

We end this section by a final comment. The geometric conservation law (21) reduces to

W Σ J + F ⋆ t U • N Σ = 0.
Thanks to (49), this relation is equivalent to

W Σ J + U • F ⋆ N Σ = 0. ( 50
)
This can be recovered from the Hadamard relation. Noticing that F ⋆ t F = J I d , we dot multiply (44a) by F ⋆ N Σ , and get

W Σ ( F ⋆ N Σ ) • ( F N Σ ) + ( F ⋆ N Σ ) • U = 0.
Since F ⋆ N Σ = 0, the first term is indeed W Σ F ⋆ t F and thus we recover (50).

Physical conservation laws

Standard results applied to the mass conservation (32a) leads to

W Σ ρJ = 0. ( 51
)
The mass swept by the discontinuity is

M = ρ ± J ± W Σ .
The relation ( 51) is of course strictly equivalent to the jump relation on mass in Eulerian coordinates. Indeed, defining m = ρ ± (U ±w σ ) • n σ as the Eulerian mass flux, where n σ is the unit normal to the discontinuity surface σ in the actual configuration, the following relation holds

M dS = -m ds. ( 52
)
This relation states that the mass crossing the discontinuity surface is identical regardless the configuration employed. The minus sign in the above equation is due to the fact that m stands for the mass flux crossing the discontinuity in the Eulerian frame whereas M stands for the mass flux swept by the moving discontinuity in the Lagrangian frame. This relation (52) can also be seen as a consequence of Hadamard's relation. Indeed, by means of relation (48), it yields

ρ ± (U ± -w σ ) • n σ ds = -ρ ± F ± W Σ • n σ ds. ( 53
)
According to Nanson formula, we have

n σ ds = F ⋆ N Σ dS. ( 54
)
And recalling that the normal in the Eulerian frame is continuous, as F ⋆ N Σ = 0, two new relations arise

F ± t n σ ds = J ± N Σ dS.
Substituting this last result in equation (53) leads to

ρ ± (U ± -w σ ) • n σ ds = -ρ ± J ± W Σ dS. ( 55
)
Thus, by means of Hadamard's relation and Nanson formula, we recover the previous relation (52).

For the momentum and total energy conservation equations, the jump relationships write

ρJU -P F ⋆ N Σ = 0, (56)
ρJE -P U • F ⋆ N Σ = 0. ( 57
)
We also recall the jump relation for the Jacobian

J W Σ + U • F ⋆ N Σ = 0.
Using (51), the above equations become

M 1 ρ + U • F ⋆ N Σ = 0, (58a)
M U -P F ⋆ N Σ = 0, (58b)
M E -P U • F ⋆ N Σ = 0. (58c) Dot-multiplying (58b) by F ⋆ N Σ , one gets M U • F ⋆ N Σ -P ( F ⋆ N Σ ) 2 = 0, (59)
and thus, employing (58a) multiplied by M leads to

M 2 = - P 1 ρ ( F ⋆ N Σ ) 2 . ( 60
) Since E = ε + 1 2 U 2 , one gets E = ε + U • U ,
and by means of this last expression as well as relations (58a) and (58c), it follows that

-M ε -M P 1 ρ = 0.
Assuming that the discontinuity is a shock wave, M = 0, we finally conclude that

ε + P 1 ρ = 0. ( 61
)
This is the Hugoniot relationship which defines the shock wave in the thermodynamic plane.

We end this section by some well-known statements concerning the jump of velocity across a discontinuity. First, we recall relation (59)

M U • F ⋆ N Σ -P ( F ⋆ N Σ ) 2 = 0. ( 62
)
Knowing that F ⋆ N Σ is nothing but the normal in the actual configuration, this last relation states that the jump of normal velocity through a discontinuity is proportional to the jump in pressure.

Now, let T be a vector orthogonal to N Σ , and a vector t defined as t = F T . The use of the Nanson formula (54) and of the continuity relations expressed through F T = 0 and F ⋆ N Σ = 0, allows us to write

t • n σ dS = F T • F ⋆ N Σ dS, = J T • N Σ dS, = 0,
which states that t = F T is a vector in the actual configuration orthogonal with the normal n σ . This result also pictures the previous continuity relation F T = 0. Finally, dot-multiplying (58b) by t, it follows that

M U • t = J P T • N Σ = 0. ( 63
)
This final statement permits us to conclude that the tangential velocity is continuous on a discontinuity. These two last results will be very helpful in the design of the numerical flux responsible of the motion of the computational grid.

In this section, we have touched on the main considerations inherent of the two type of descriptions of flows, namely the Lagrangian description and the Eulerian description. Doing so, some crucial identities have arisen, as the Piola compatibility conditions stated in [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems[END_REF]. We have then derived the geometrical and physical conservation laws governing the time evolution of the different geometrical and physical quantities involved, in both Lagrangian and Eulerian frameworks. We have seen that under a fundamental assumption on the Jacobian, i.e., det F > 0, ensuring that the flow map is invertible, both formulations are perfectly consistent. Finally, this section has been ended by the establishment of jump relations through a moving discontinuity. All these fundamental results provide us with a complete framework and with leading constraints which will help us in the implementation of a numerical scheme approximating the solutions of the gas dynamics equations written under a total Lagrangian formalism.

High-order discontinuous Galerkin cell-centered Lagrangian schemes

Here, we present a general high-order discontinuous Galerkin discretization of the gas dynamics equations written under the total Lagrangian formulation, while the Lagrangian-Eulerian map is spanned by a finite element representation. This provides an approximation of the deformation gradient tensor satisfying the Piola identity. The velocity field having a high-order polynomial representation with respect to the space variables, the mesh edges in the actual configuration are parametrized by means of one-dimensional finite element basis functions, consistent with the mapping representation. The DG discretization of the physical conservation laws for the specific volume, the momentum and the total energy are performed ensuring the respect of the GCL.

Discontinuous Galerkin general framework

We intend to solve the gas dynamics equations written under the total Lagrangian form

∂ F ∂t -∇ X U = 0, (64a)
ρ 0 ∂ ∂t ( 1 ρ ) -∇ X (F ⋆ t U ) = 0, ( 64b
)
ρ 0 ∂ U ∂t + ∇ X (P F ⋆ ) = 0, (64c)
ρ 0 ∂ E ∂t + ∇ X (F ⋆ t P U ) = 0. ( 64d
)
Here, J = det F is the Jacobian, and F ⋆ = JF -t is the cofactor matrix of the deformation gradient tensor F. We have seen previously that the specific volume equation (64b) of a Lagrangian cell is a consequence of the mass conservation and the geometrical equation (64a), see section 2.5.1. Nonetheless, in the development of our numerical scheme these two equations will be treated in different ways, equation (64a) being considered as purely geometric will be discretized in some continuous manner, whereas equation (64b) will be treated in a consistent way with the momentum and total energy equations, by mean of a discontinuous Galerkin approach. Anyhow, we shall demonstrate later that the discretizations of these two equations are consistent. We also will be able to design a scheme ensuring a local entropy inequality, and characterized by a high-order accuracy.

The thermodynamical closure of system (64) is obtained through the use of an equation of state, which writes P = P (ρ, ε), where ε is the specific internal energy, ε = E -1 2 U 2 . These equations are valid provided that the Lagrangian-Eulerian flow map exists, that is J > 0. In this framework, the computational grid is fixed, however one has to follow the time evolution of the Jacobian matrix associated to the Lagrange-Euler flow map. We made the choice of working on the initial configuration of the flow to avoid some difficulties inherent to the moving mesh scheme, as dealing with curvilinear geometries, in the case of very high-order scheme. Furthermore, in this frame the basis functions are time independent and defined on the initial mesh. This mesh being perfectly known, geometrical quantities required by the DG discretization can be computed and stored initially and used during the whole calculation.

Let {Ω c } c be a partition of the domain Ω into non-overlapping polygonal cells. We also partitionate the time domain in intermediate times (t n ) n with ∆t n = t n+1 -t n the n th time step. We use a DG discretization in order to develop on each cells our unknowns onto P s (Ω c ), the set of polynomials of degree up to s. This space approximation leads to a (s + 1) th space order accurate scheme. Let φ c h be the restriction of φ h , the polynomial approximation of a function φ, over the cell Ω c

φ c h (X, t) = K k=0 φ c k (t) σ c k (X), ( 65
)
where the φ c k are the K + 1 successive components of φ h over the polynomial basis, and σ c k the polynomial basis functions. Recalling that in the two-dimensional case the dimension of the polynomial space P s (Ω c ) is (s+1)(s+2)

2

, one has to select a set of (s+1)(s+2) 2 = K + 1 basis functions. The establishment of a DG discretization is based on a local variational formulation of the equations.

For sake of simplicity, we introduce the following generic equation

ρ 0 ∂ φ ∂t + ∇ X (F ⋆ t f ) = 0, ( 66
)
where f is the flux function associated to the unknown φ, where both functions are defined using the Lagrangian coordinates. Then, we first focus our study on the discretization of such an equation, using DG method. A local variational formulation of this equation and substituting the function φ by its piecewise polynomial approximation φ c h results in

Ωc ρ 0 ∂ φ c h ∂t σ c j dV = K k=0 d φ c k dt Ωc ρ 0 σ c j σ c k dV, = Ωc f • F ⋆ ∇ X σ c j dV - ∂Ωc f • σ c j F ⋆ N dS, (67)
where σ c j is a function picked into the chosen basis {σ c k } k=0...K of P K , and f is the continuous numerical flux function defined at cell interfaces. Here, we can identify Ωc ρ 0 σ c j σ c k dV = (σ c j • σ c k ) c as the generic coefficient of the symmetric positive definite mass matrix, where the scalar product is defined as follow

(φ • ψ) c = Ωc ρ 0 (X) φ(X) ψ(X) dV. ( 68
)
Let us note that this scalar product is weighted by the initial mass density. Consistent with this scalar product definition, we introduce the mass averaged value as

φ c = 1 m c Ωc ρ 0 (X) φ(X) dV. ( 69
)
In our work, the choice of using Taylor basis has been made to construct the polynomial approximation of our numerical solutions, (65). These particular basis functions, introduced in [START_REF] Luo | A DG method based on a Taylor basis for the compressible flows on arbitrary grids[END_REF], are based on a Taylor expansion at the centroid of the cells. The main interest of this basis lies in the fact that it does not depend on the shape of the cell. Namely, we shall employ the same basis for triangular and quadrangular cells. Moreover, with this basis, we will be able to construct a DG discretization over general polygonal grids. Practically, the Taylor basis used here are slightly different from the one presented in [START_REF] Luo | A DG method based on a Taylor basis for the compressible flows on arbitrary grids[END_REF]. Knowing that the mass matrix is weighted by the initial density, we derive our basis function by a Taylor expansion with respect to the center of mass,

X c = (X c , Y c ) t , of the cell c φ(X) = φ(X c ) + s k=1 k j=0 (X -X c ) k-j (Y -Y c ) j j!(k -j)! ∂ k φ ∂X k-j ∂Y j (X c ) + o( X -X c s ), ( 70
)
where the center of mass is defined by

X c = 1 m c Ωc ρ 0 (X) X dV, (71)
m c being the constant mass of cell Ω c . Let us point out that in the case of a uniform initial density, the center of mass and the centroid coincide. In equation (70), we identify the term (X-Xc) k-j (Y -Yc) j j!(k-j)!

to the successive basis functions. To ensure a conservative discretization, we set the first moment equal to the mass averaged value, φ c 0 = φ c , the first basis function is equal to one, σ c 0 = 1, and we impose that σ c k c = 0, ∀k = 0. Consequently, the k th -order polynomial components of the approximated function and the associated basis function read

φ c k(k+1) 2 +j = (∆X c ) k-j (∆Y c ) j ∂ k φ ∂X k-j ∂Y j (X c ), (72)
σ c k(k+1) 2 +j = 1 j!(k -j)! X -X c ∆X c k-j Y -Y c ∆Y c j - X -X c ∆X c k-j Y -Y c ∆Y c j c , (73)
where 0 < k ≤ s, j = 0 . . . k, and ∆X c = Xmax-X min

• σ c k ) c = m c δ 0k
, where δ ij , the Kronecker symbol, is equal to 1 if i = j and null otherwise. Recalling that (σ c j •σ c k ) c identifies with the generic coefficient of the mass matrix, the equations corresponding to mass averaged values does not depend on the other polynomial basis components equations, and correspond to a first-order finite volume scheme.

Recalling the variational equation (67), one can see the contribution of two different terms, the interior term Ωc f • F ⋆ ∇ X σ c j dV , and the boundary term ∂Ωc f • σ c j F ⋆ N dS. In [START_REF] Cockburn | The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case[END_REF], it has been demonstrated that to design a (s + 1) th order numerical scheme, a quadrature rule over the faces being exact for polynomials of degree 2s + 1 is needed, as a quadrature rule over the elements being exact for polynomials of degree 2s. Consequently, to evaluate the interior terms, a two-dimensional high-order quadrature rule is employed. Nevertheless, a specific treatment of the boundary terms is required to ensure the compatibility with the geometrical conservation law (GCL), equation [START_REF] Loubère | Une Méthode Particulaire Lagrangienne de type Galerkin Discontinu[END_REF]. Indeed, substituting φ by 1 ρ and f by -U in equation (67), for j = 0, yields

m c d dt ( 1 ρ ) c 0 = Ωc ρ 0 ∂ ∂t ( 1 ρ ) dV = ∂Ωc U • F ⋆ N dS, (74)
where U , the continuous numerical flux, is nothing but the velocity responsible for the motion of the grid. We have to ensure that the cells volume computed by means of the moving grid geometric informations are exactly coincident with the cells volume comupted solving the specific volume equation. This shall require a careful treatment of the numerical fluxes attached to cell interfaces.

Let us recall that in the particular case of an uniform flow, the gas dynamics equations collapse to the simple advection of the density field. Considering (67) in the uniform case, i.e., f is constant, the left-hand side of the equation has to be equal to zero, assuming a consistent numerical flux f . This means that the following relationship must hold for all basis functions

Ωc F ⋆ ∇ X σ c j dV = ∂Ωc σ c j F ⋆ N dS, (75)
which is nothing but a generalization of the weak form of the Piola compatibility condition equation (13b) as

Ωc σ c j (∇ X F ⋆ ) dV = 0. (76)
Assuming that the Nanson formula [START_REF] Cheng | A high order eno conservative lagrangian type scheme for the compressible euler equations[END_REF] and gradient relation [START_REF] Gresho | On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via finite-element method that also introduces a nearly consistent mass matrix[END_REF] hold after discretization of the deformation gradient tensor F, and then rewriting identity (75) in the actual configuration yields a well-known result, consequence of the Green-Gauss formula, which reads ωc

∇ x σ c j dv = ∂ωc σ c j nds. ( 77
)
Observing identity (75), one clearly sees that some consistency is required on the discretization of the deformation gradient tensor F inside the cell and on its boundaries. Furthermore, one knows that at the continuous level the consistency between the two configurations lies, among other things, on two fundamental relations, the Nanson formula [START_REF] Cheng | A high order eno conservative lagrangian type scheme for the compressible euler equations[END_REF] and the Piola compatibility condition (13b). Consequently, particular considerations on the discretization of F are required to ensure that these relations hold at the discrete level. Obviously, making use anew of a DG discretization to approximate the deformation gradient tensor will failed to ensure these essential requirements on the geometry. These are the reasons that have motivated our choice of discretizing the tensor F by means of a mapping using finite element basis on triangular cells.

Semi-discrete equation for the deformation gradient tensor

Let Ω c be a generic polygonal cell in the initial configuration which is paved into non-overlapping triangles T c i as Ω c = ntri i=1 T c i , refer to Figure 2. Recalling the mapping formulation expressed in (1), we aim at developing a (s + 1) th order continuous approximation of this flow map function, in a generic triangle T i . To this end, we shall use the following finite element approximation of the mapping s as

Ω c T c i

Φ i h (X, t) = q∈Q(i) Λ i q (X) Φ q (t). (78)
Here, Λ i q is a finite element basis function of degree s and Q(i) is the control points set of triangle T i , including its vertices {p -, p, p + }. We identify Φ q (t) = Φ(X q , t) as the position at time t of the control point initially located at X q . The functions Λ i q being finite element basis functions, they identify with the P 1 barycentric coordinate basis functions in the case of our second-order scheme, and with the P 2 finite elements basis functions in the third-order case.

Using the continuous polynomial mapping approximation (78) and the trajectory equation ∂ ∂t Φ q = U q we obtain our semi-discretization of the deformation gradient tensor

∂ ∂t F i h (X, t) = q∈Q(i) U q (t) ⊗ ∇ X Λ i q (X), (79)
where U q is the velocity of the control point q. In [START_REF] Kluth | Discretization of hyperelasticity on unstructured mesh with a cellcentered Lagrangian scheme[END_REF], to develop their cell-centered Lagrangian scheme for the hyperelasticity, the authors make use of a similar discretization of the deformation gradient tensor. This particular definition of F ensures naturally the Piola compatibility condition.

Knowing that initially F ⋆,i h (X, 0) = F i h (X, 0) = I d , the Piola condition being satisfied, it sufficient to prove that the condition is involutive. This result follows immediately in the two-dimensional case, the cofactor matrix F ⋆ writing

∂ ∂t F ⋆,i h = q∈Q(i) V q ∂ Y Λ i q -V q ∂ X Λ i q -U q ∂ Y Λ i q U q ∂ X Λ i q , ( 80
)
where U and V are respectively the X and Y components of the velocity U . And taking the divergence of equation ( 80), one gets

∂ ∂t (∇ X F ⋆,i h ) = q∈Q(i) (∂ XY Λ i q -∂ Y X Λ i q ) -V q U q = 0. (81)
The same conclusion can be easily obtained in the three-dimensional case, the cofactor matrix being defined as ∂ ∂t F ⋆,i h = p,q Ψ pq ⊗ ∇ X Λ i p × ∇ X Λ i q , where the components of vector Ψ pq reads Ψ X pq = V p W q , Ψ Y pq = W p U q and Ψ Z pq = U p V q , where W being the Z components of the velocity U .

Using such a discretization (79), it is clear that even if the basis function Λ i q are continuous over the whole domain, the deformation gradient tensor being expressed through the gradient of these functions, F will be piecewise continuous over the triangular cells. Anyhow, thanks to the mapping formulation and hence to equation (79), the requirements enlightened previously on the consistency and continuity of F ⋆ N on the triangles boundaries will be ensured by construction of the scheme. Thanks to these properties, we are now able to show the consistency between the two different discretizations of the geometrical conservation law, (64a) and (64b). Let us define the continuous flow motion velocity U , using the approximated flow map (78) and the trajectory equation ∂ ∂t Φ q = U q , such as

U | T i (X, t) = q∈Q(i) Λ i q (X) U q (t), (82)
where U | T i is the restriction of the continuous velocity U on T i . This continuous velocity identifies on the polygonal cell boundaries with the numerical flux velocity introduced in equation (74). This definition allows us to rewrite (79), as

∂ ∂t F i h (X, t) = ∇ X U | T i (X, t). ( 83
)
Making use of ( 25) and (83), we can now express the time rate of change of the discrete Jacobian J i h (X, t) mean value over the triangle T i as

d dt T i J dV = T i F ⋆ : ∇ X U dV.
In this last equation, we dropped the subscript h and superscript i for the sake of simplicity.

Recalling that the Piola compatibility condition is ensured at the discrete level, everywhere in the domain, this last expression rewrites

d dt T i J dV = T i ∇ X (F ⋆ t U ) dV, = ∂T i U • F ⋆ t N dS. ( 84
)
Now, making use of equation ( 74), and thanks to the continuity of F ⋆ t N at the triangle boundaries, it easily follows that

d v c dt = ∂Ωc U • F ⋆ t N dS, = ntri i=1 ∂T i U • F ⋆ t N dS,
and we finally get the following relationship ensuring the consistency of the two geometric conservation law discretizations such as

d v c dt = ntri i=1 d dt T i J dV. ( 85
)
Regarding (67), the last points which remains to be addressed is the definition of the numerical fluxes and the integration of the boundary terms.

Entropic analysis

We design the numerical fluxes of our scheme in such a way that a local entropy inequality is satisfied at the semi-discrete level. This approach, which ensures kinetic energy conversion into internal energy through shock waves, is similar to the one used in our previous paper [START_REF] Vilar | Cell-centered discontinuous Galerkin discretizations for two-dimensional scalar conservation laws on unstructured grids and for one-dimensional Lagrangian hydrodynamics[END_REF][START_REF] Vilar | Cell-centered discontinuous Galerkin discretization for two-dimensional Lagrangian hydrodynamics[END_REF]. Let η be the specific entropy and θ > 0 the absolute temperature defined by means of the Gibbs identity as follows

θdη = dE -U • dU + P d( 1 ρ ). ( 86
)
We aim at expressing the time rate of change of the specific entropy. Here, for the sake of conciseness we identify the functions 1 ρ , U , E and P to their polynomial approximation over the considered cell Ω c , respectively ( 1 ρ ) c h , U c h , E c h and P c h . Firstly, multiplying the gas dynamics equations (64b),(64c) and (64d) by respectively P , U and 1 as test functions and integrating by parts leads to the local variational formulations which express the pressure work, the kinetic energy and the total energy. Secondly, combining these variational formulations according to the Gibbs formula (86) yields

Ωc ρ 0 θ ∂ η ∂t dV = ∂Ωc P U + P U -P U • F ⋆ N dS - Ωc [P (F ⋆ : ∇ X U ) + U • F ⋆ ∇ X P ]dV.
Here, the overbars above the variables P , U and P U denote the numerical fluxes which are continuous across cell interfaces. It remains to simplify the second term in the right-hand side of the above equation using the Piola compatibility condition, ∇ X • (F ⋆ ) = 0, and recalling the identity

∇ X (F ⋆,t P U ) = U • ∇ X • (F ⋆ P ) + P F ⋆ : ∇ X U .
Finally, we get obtain the following expression of the time rate of change of the specific entropy

Ωc ρ 0 θ ∂ η ∂t dV = ∂Ωc P U + P U -P U -P U • F ⋆ N dS. ( 87
)
At this point, it remains to express the numerical fluxes in such way that an entropy inequality is satisfied. To this end, we first make the following fundamental assumption

P U = P U . ( 88
)
This assumption allows us to factorize the right-hand side of equation ( 87)

Ωc ρ 0 θ ∂ η ∂t dV = ∂Ωc (P -P )(U -U ) • F ⋆ N dS. (89)
Finally, to enforce a local entropy inequality at the semi-discrete level, we prescribe the following sufficient condition on the numerical fluxes

P -P = -Z (U -U ) • F ⋆ N F ⋆ N = -Z (U -U ) • n, ( 90
)
where Z is a positive constant which has the physical dimension of a density times a velocity. For the numerical applications, we use Z = ρ a, a being the thermodynamic sound speed, which corresponds to the acoustic impedance. In the end, we have derived a particular form of the numerical fluxes ensuring a local entropy inequality at the semi-discrete level, which is also perfectly consistent with the jump relation stated in (62). Now, to enforce the respect of the geometric conservation law we detail the particular treatment of the boundary term in (67).

Nodal solvers located at the control points

We start this section related to the boundary terms integration by some geometric considerations. We have seen in section 3.2 how we develop the continuous mapping function onto finite element basis functions. The continuous polynomial function (78) describes the flow motion of the fluid. Thus, a straight line edges triangle T i in the initial configuration will be deformed through the fluid flow into a triangle τ i in the actual configuration, as displayed in Figure 3, in the particular case of a fourth order scheme.

We state here that the definition of the curvilinear edges of triangle τ i in the actual configuration can be defined by means of the trace over ∂T i of the finite element basis functions used in (78). Thus, to parametrize the face f pp + in the actual configuration, we project the continuous mapping function Φ i h over f pp + . It leads to the definition of the curved edge by means of s + 1 control points, which reads

p - p + p + p - p p Φ T i τ i
x | pp + (ζ) = q∈Q(pp + ) x q λ q (ζ) = x p λ p (ζ) + q∈Q(pp + )\{p,p + } x q λ q (ζ) + x p + λ p + (ζ), (91)
where Q(pp + ) is the set of the face control points, ζ ∈ [0, 1] being the curvilinear abscissa and λ q the one-dimensional Lagrangian finite element basis functions of degree s. In the second-order case, s + 1 = 2 and thus the triangle edges are defined only using the vertices of the triangle {p -, p, p + }, which is perfectly natural, the edges remaining straight line edges in this case. In the third-order case, one more control point is required on the face to define the curved edge. The triangle edges then identify with Bezier curves.

Here, we recall the local variational formulation of the gas dynamics equation type

Ωc ρ 0 ∂ φ c h ∂t σ c j dV = - Ωc F ⋆ ∇ X σ c j • f dV + p∈P(c) p + p f • σ c j F ⋆ N dL, (92)
where P(c) is the vertices set of the cell Ω c . The corner stone of the scheme consists in constructing the numerical fluxes through some point contributions. To do so we assume the numerical fluxes to be polynomial functions of the same degree than the piecewise polynomial approximations of the unknowns, and than the approximated continuous mapping function as well. Such an assumption allows us to express the numerical flux f on face f pp + through the use of s + 1 point contributions as

f | pp + (ζ) = f + pc λ p (ζ) + q\{p,p + } f qc λ q (ζ) + f - p + c λ p + (ζ). ( 93
)
It can be noted in this definition of the numerical flux that the control points contribution is local to the cell, and also to the face as there is left and right contributions at the vertices, see Figure 4.

Using the fact that the Taylor basis functions employed are of the same degree than f , one can express their trace on face f pp + by means of their interpolated values at the different face control points, such as And as we know that the discretization of the deformation gradient tensor enforces the relation F ⋆ N dL = ndl, one can write

σ c j| pp + (ζ) = σ c j (X p )λ p (ζ) + q\{p,p + } σ c j (X q )λ q (ζ) + σ c j (X p + )λ p + (ζ). ( 94
) p + p - p N pp + N p -p f - pc f + pc f - p + c f qc Ω c
F ⋆ N dL | pp + (ζ) = n dl | pp + = ∂x ∂ζ dζ | pp + × e z = q∈Q(pp + ) ∂λ q ∂ζ (ζ) (x q × e z ) . (95)
We have finally expressed the different terms present in the boundary integral as functions of the curvilinear abscissa ζ. Consequently, we are now able to integrate analytically the boundary term. This is the general procedure used to derive the semi-discrete scheme. Nevertheless, some further assumptions are needed before applying this procedure to the specific volume, momentum and total energy equations. The first one comes from the fact that the numerical flux U is responsible of the motion of the grid. The different contributions of the numerical flux around a control point are required to be coincident to maintain the integrity of the grid. Such an assumption writes

U ± pc = U p , ∀c ∈ C(p) and U qL = U qR = U q , (96)
where C(p) represents the set of cells surrounding the point p, and Ω L and Ω R the two neighboring cells sharing the face control point q. The other assumption will allow us to treat the total energy equation in a consistent manner with the other ones. According to assumption (88), we prescribe

(P U ) ± pc = P ± pc U p and (P U ) qc = P qc U q . ( 97
)
Gathering the different curvilinear definitions and specific assumptions of the numerical fluxes, we are now ready to integrate analytically the boundary terms present in the local variation formulations of the specific volume, momentum and total energy equations. At the end, performing an index permutation allows us to get face control points solver, as well as nodes solver containing left and right contributions. The semi-discrete scheme finally writes

Ωc ρ 0 ∂ ∂t ( 1 ρ ) c h σ c j dV = - Ωc U • F ⋆ ∇ X σ c j dV + p∈P(c) (U p • l j pc n j pc + q\{p,p + } U q • l j qc n j qc ), ( 98
) Ωc ρ 0 ∂ U c h ∂t σ c j dV = Ωc P F ⋆ ∇ X σ c j dV - p∈P(c) (F j pc + q\{p,p + } F j qc ), ( 99
) Ωc ρ 0 ∂ E c h ∂t σ c j dV = Ωc P U • F ⋆ ∇ X σ c j dV - p∈P(c) (U p • F j pc + q\{p,p + } U q • F j qc ), (100)
where the j th moment of the subcell forces, F j pc and F j qc , reads F j pc = P - pc l -,j pc n -,j pc + P + pc l +,j pc n +,j pc and F j qc = P qc l j qc n j qc .

The nodal weighted normals are defined as l j pc n j pc = l -,j pc n -,j pc + l +,j pc n +,j pc , l +,j pc n +,j pc =

1 0 λ p| pp + (ζ)σ j| pp + (ζ) ∂x ∂ζ dζ | pp + × e z , l -,j pc n -,j pc = 1 0 λ p| p -p (ζ)σ j| p -p (ζ) ∂x ∂ζ dζ | p -p × e z ,
whereas the face control point weighted normals writes

l j qc n j qc = 1 0 λ q| pp + (ζ)σ j| pp + (ζ) ∂x ∂ζ dζ | pp + × e z .
Thanks to this specific treatment of the boundary terms, the semidiscretization of the specific volume equation (98) ensures the respect of the GCL by construction of the scheme. In the end, we make use of the conclusion of the previous entropy analysis presented in section 3.3 to express to control point solvers F j pc , F j qc , and U p , U q . Recalling relation (90)

P -P = -Z (U -U ) • n,
leads to define it follows the subcell forces according to

F j pc = P c h (X p , t) l j pc n j pc -M j pc (U p -U c h (X p , t)), (102)
F j qc = P c h (X q , t) l j qc n j qc -M j qc (U q -U c h (X q , t)), (103)
where the matrices are defined as

M j pc = Z c l -,j pc n -,j pc ⊗ n -,0 pc + l +,j pc n +,j pc ⊗ n +,0 pc , (104)
M j qc = Z c l j qc n j qc ⊗ n 0 qc . ( 105
)
One can clearly see in (104) the two edges contribution at a cell node. We know from section 3.1 that the equations related to the averaged values of the unknowns are independent of the successive derivative ones. This property of the chosen basis permits us to enforce conservation relations just focusing on the simple equations related to the first moments. Consequently, an obvious way to ensure the conservation of momentum and total energy over the domain is to set the following constraints

c∈C(p) F 0 pc = 0 and F 0 qL + F 0 qR = 0. ( 106
)
Thanks to (106), we finally have an explicit expression of the nodal velocity U p

M p U p = c∈C(p) P c h (X p , t) l 0 pc n 0 pc + M 0 pc U c h (X p , t) , (107)
where M p = c∈C(p) M 0 pc is a positive definite matrix. An identical procedure leads the following definition of the face control point velocity

M q U q = M q Z L U L h (X q ) + Z R U R h (X q ) Z L + Z R - P R h (X q ) -P L h (X q ) Z L + Z R l 0 qL n 0 qL , (108)
where

M q = 1 Z R M 0 qR = 1 Z L M 0 qL = l 0 qL n 0 qL ⊗ n 0 qL is positive semi-definite.
Actually, this matrix is non-invertible. So the only information deriving from this equation is the normal velocity definition

(U q • n 0 qL ) = Z L U L h (X q ) + Z R U R h (X q ) Z L + Z R • n 0 qL - P R h (X q ) -P L h (X q ) Z L + Z R . ( 109
)
This is nothing but the solution of the one-dimensional Riemann problem at the cell interface located at the face control point q. To define the tangential velocity of the face control points, we make use of the results of the jump relations presented in section 2.8. It has been stated in equation (62) that at a discontinuity the gap in normal velocity is proportional to the gap in pressure. The definition of the normal velocity (109) is consistent with this statement, the first term identifying with an average value and the second one expressing the gap in pressure. The Rankine-Hugoniot relations also tell us that the tangential velocity is continuous at a discontinuity, see equation ( 63). This is the reason why we define the tangential velocity of the face control point using the same continuous part that in (109), as

(U q • t 0 qL ) = Z L U L h (X q ) + Z R U R h (X q ) Z L + Z R • t 0 qL . (110)
Gathering the above results, the face control point velocity writes

U q = Z L U L h (X q ) + Z R U R h (X q ) Z L + Z R - P R h (X q ) -P L h (X q ) Z L + Z R n 0 qL . ( 111
)
Regarding the semi-discrete equation of the deformation gradient tensor in the triangles constituting the polygonal cells, equation (79), we obviously need to define the velocity of the interior points, refer to Figure 5. Any low order assumption on the definition of these interior points velocity will lead to a decrease of accuracy. Consequently, we set the velocity of an interior point q of cell Ω c to the interpolated value at this point of the polynomial approximation of the velocity inside the cell, i.e., U q = U c h (X q ), where X q is the initial position of point q. In this section, we have presented a general way to design high-order cell-centered Lagrangian semi-discrete schemes, allowing to cope with curvilinear geometries, and ensuring the GCL as well as the Piola compatibility condition by construction of the scheme. This particular discretization also enforce a local entropy inequality and global conservation of moment and total energy. Now, we show that under specific choices and assumptions one can obtain a consistent discretization on the actual moving configuration.

Compatibility between discretizations on initial and actual configurations

We recall the Lagrangian gas dynamics equation expressed in terms of Eulerian coordinates derived in section 2.5

ρ d dt ( 1 ρ ) -∇ x U = 0, ( 112a
)
ρ d U dt + ∇ x P = 0, (112b)
ρ d E dt + ∇ x (P U ) = 0. ( 112c
)
Let us introduce the generic equation

ρ d φ dt + ∇ x f = 0, ( 113
)
where f is the flux function associated to the unknown φ, where both functions are defined using the Eulerian coordinates. The use of DG discretization aims at approximating the function φ by

φ c h (x, t) = K k=0 φ c k (t) ς c k (x, t), (114)
where φ c h is the restriction on the cell ω c of the polynomial approximation φ h of function φ, and ς c k the chosen basis functions defined on the moving grid. A local variational formulation of this equation and substituting the function φ by its piecewise polynomial approximation φ c h results in

ωc ρ d φ c h dt ς c j dv = ωc f • ∇ x ς c j dv - ∂ωc f • ς c j nds, ( 115
)
where f is the continuous numerical flux function defined at cell interfaces. We have seen in section 2 that the compatibility of both formulations at the continuous level relies mainly on the Nanson formula [START_REF] Cheng | A high order eno conservative lagrangian type scheme for the compressible euler equations[END_REF] and the Piola compatibility condition (13b), namely on the deformation gradient tensor definition. Regarding the semi-discretization presented in this paper, it has been proven that these two identities are ensured by construction of the scheme. Furthermore, thanks to equation ( 83) one can see that, under a correct time integration, the discrete deformation gradient tensor derives directly from the motion of the grid. Thus, assuming a polynomial fluid flow motion of degree s, the relations obtained at the continuous level transforming points, surfaces, volumes, as well as gradient and divergence operators, from the actual configuration to the initial one hold at the discrete level. In the framework of moving grid schemes, this assumption on the fluid flow is implicit and is governed by the grid edges parametrization. Finally, the remaining ingredient to ensure an equivalence between the discretizations is a wise choice of basis functions based on the moving configuration to ensure that the successive moments φ c k identify with the ones obtained in our discretization. Actually, the single situation ensuring such a consistency is the case where the basis functions are conserved along the trajectories, i.e., d ς c j dt = 0. This property is naturally ensured in the case of point-based finite elements functions, the control points being advected by means of the fluid velocity. But in the case of the chosen Taylor basis, if we want the functions to follow the fluid flow we have to define ς c j such as

ς c j (x, t) = σ c j (X) = σ c j (Φ -1 (x, t)), (116)
where Φ -1 is the invert function of the mapping. It seems that such a choice would not be relevant to develop a scheme on the actual configuration, this method requiring the storing of the initial position of each points during the whole calculation. Actually, only a finite number of points is required. It is sufficient to store the initial position of each control points inside and on the boundaries of the polygonal cells, plus a certain number of quadrature points needed in the evaluation of the interior terms. Anyhow, concerning this latter point, any high-order DG scheme on moving mesh would require a particular treatment of the interior term to maintain the accuracy.

In the end, we conclude that it is possible to derive a high-order DG scheme on moving mesh approximating the solutions of the Lagrangian gas dynamics equations (112), which leads to the same results than the one obtained using the initial configuration scheme presented in this paper. This consistency lies on a specific definition of the basis functions in the Eulerian frame, see definition (116), and is possible only because in the presented scheme the deformation gradient tensor derives precisely from the motion of the grid. It is worth mentioning that in the whole scheme implementation presented in this paper, a straight line edges initial grid has been considered. Nevertheless, one may want to start with curvilinear cells. Next, we present the procedure which allows to take into account such grids.

Specific procedure devoted to curvilinear grids

Let Ω c be a cell of the initial configuration Ω at time t = 0, and ω c its corresponding cell in the actual configuration domain ω at time t. We assume that the initial domain is paved using a curvilinear grid. In that case, the scheme discretization remain correct, one just has to initialize correctly the deformation gradient tensor and the mass matrix. Let the straight polygonal cell Ω r c be the referential cell of Ω c , in the referential domain Ω r , refer to Figure 6.

Φ H (X, t) Φ T (X r , t) Ω c ω c Ω r c X r X x Φ 0 (X r )
Figure 6: Reference, initial and actual configuration.

X denotes the initial position of the moving point located at time t at the position x. The mathematical transformation Φ H (X, t) represents the displacement due to the flow motion of the fluid. X r is the referential position of the point initially located at X. Φ 0 (X r ) represents the initial transformation of the domain and Φ T (X r , t) corresponds to the total deformation mapping, refer to Figure 6. Gathering the above notation leads to define Φ T as the composition of Φ H and Φ 0 as follows

Φ T (X r , t) = Φ H (X, t) • Φ 0 (X r ).
The use of the chain rule of composed derivatives and of the deformation gradient tensor definition yields

F T = ∇ Xr Φ T (X r , t), = ∇ X Φ H (X, t) • ∇ Xr Φ 0 (X r ), = F H F 0 ,
where

F 0 = ∇ Xr Φ 0 (X r ) and F H = ∇ X Φ H (X, t).
Taking the determinant of the tensors on both sides of the above equation yields

J T (X r , t) = J H (X, t) J 0 (X r ), (117)
where J T = det F T , J H = det F H and J 0 = det F 0 . To recover the Lagrangian equation of continuity, we recall the mass conservation principle Ωc ρ 0 (X) dΩ = ωc ρ(x, t) dω.

Expressing the integrands in terms of the referential coordinates yields

Ω r c ρ 0 (Φ 0 (X r )) J 0 (X r ) dΩ r = Ω r c ρ(Φ T (X r , t), t) J T (X r , t) dΩ r .
From this equation, it follows easily that

ρ 0 J 0 = ρ J T . (118)
This new continuity equation only affects the construction of the mass matrix. Indeed, trying to calculate ωc ρ dψ dt σ c j dω, the time rate of change of the successive moments of the function ψ, one gets

ωc ρ d ψ dt σ c j dω = Ω r c ρ J T ∂ ψ ∂t σ c j dΩ r , = Ω r c ρ 0 J 0 ∂ ψ ∂t σ c j dΩ r .
And substituting the function ψ by its polynomial approximation,

ψ c h (X r , t) = K k=0 ψ c k (t)σ c k (X r ), it yields ωc ρ d ψ c h dt σ c j dω = K k=0 d ψ k dt Ω r c ρ 0 J 0 σ c j σ c k dΩ r . (119)
Finally, in equation ( 119) we identify Ω r c ρ 0 J 0 σ q σ k dΩ r as the coefficients of the mass matrix. In the case of a non-deformed initial geometry, J 0 (X r ) = 1, the former discretization is recovered.

So far, a new high-order cell-centered Lagrangian semi-discrete scheme, allowing to cope with curvilinear grids, and ensuring the GCL as well as the Piola compatibility condition by construction, has been presented.

Limiting procedure based on the characteristic variables

In the case of discontinuous problem, without a specific treatment, high-order numerical schemes produce solutions containing spurious oscillations. On the other hand, if we apply a limiting procedure directly to the polynomial approximation of the physical variables ( 1 ρ ) h , U h and E h , we cannot enforce totally the monotonicity of the solutions. To correct this flaw, we shall use a limitation procedure which relies on the characteristic variables following the approach originally introduced by Cockburn, Lin and Shu in [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems[END_REF]. Let us assume that the fluid variables are sufficiently smooth to compute their partial derivatives. In section 2.7, it has been demonstrated that in the case of smooth flows the specific entropy is conserved along the trajectory equation, i.e., d η dt = 0. Thanks to this result, one gets the following relation between the pressure and the density material derivatives

d P dt = -ρ 2 a 2 d dt ( 1 ρ ). (120)
Gathering this last relation with the Lagrangian gas dynamics equations written in the actual configuration (112), we obtain the new non-conservative form of the gas dynamics equations

d P dt + ρa 2 ∇ x U = 0, (121a) d
U dt + 1 ρ ∇ x P = 0, (121b)
d η dt = 0. ( 121c
)
Using the material derivative definition [START_REF] Burton | A cell-centered Lagrangian Godunov-like method for solid dynamics[END_REF], and setting W = (P, U, V, η) t where U and V are respectively the x and y components of U , one can rewrite the above system in the two dimensional Cartesian frame as

∂ W ∂t + A x ∂ W ∂x + A y ∂ W ∂y = 0, (122)
where 4 × 4 matrices A x and A y are given by

A x =     U ρ a 2 0 0 1 ρ U 0 0 0 0 U 0 0 0 0 U     and A y =     V 0 ρ a 2 0 0 V 0 0 1 ρ 0 V 0 0 0 0 V     .
If n denotes a unit vector, we define

A(n) = A x n x + A y n y A(n) =     U n ρ a 2 n x ρ a 2 n y 0 1 ρ n x U n 0 0 1 ρ n y 0 U n 0 0 0 0 U n     . (123)
This matrix admits four real eigenvalues:

λ 1 = U • n -a, λ 2 = λ 3 = U • n and λ 4 = U • n + a.
One knows that the system of equations ( 122) is hyperbolic if for all n, A(n) admits four real eigenvalues. Thus, the above system is hyperbolic. Now, we are able to introduce the four differential Riemann invariants associated to unit direction n as

dα E = dE -U • dU + P d( 1 ρ ), (124)
dα -= dP -ρa dU • n, (125)
dα + = dP + ρa dU • n, (126)
dα 0 = dU • t, (127)
where t = e z × n. Recalling that for an isentropic flow the differential of the pressure expresses as dP = -ρ 2 a 2 d( 1 ρ ), the Riemann invariants differentials dα ± rewrite

dα ± = d( 1 ρ ) ± 1 ρa dU • n, (128)
Following the same procedure than in the one-dimensional case presented in [START_REF] Vilar | Cell-centered discontinuous Galerkin discretizations for two-dimensional scalar conservation laws on unstructured grids and for one-dimensional Lagrangian hydrodynamics[END_REF], we define the polynomial approximation of the Riemann invariants by linearizing the differential definitions (124), ( 127) and (128) on each cells around the mean values in the cells

α c E,h = K k=0 α c E,k σ c k = E c h -U c 0 • U c h + P c 0 ( 1 ρ ) c h , (129)
α c ±,h = K k=0 α c ±,k σ c k = ( 1 ρ ) c h ± 1 Z c U c h • n, (130)
α c 0,h = K k=0 α c 0,k σ c k = U c h • t, (131)
where φ i h is the polynomial approximation of φ on the cell Ω c , and φ c 0 its mass averaged value. The scalar Z c = a c 0 ρ c 0 is the acoustic impedance. This procedure is equivalent to linearize the gas dynamics equations, on each cells, around an averaged state. We can now apply a high-order limitation procedure on the Riemann invariant polynomial approximations. In our work, the choice has been made to use the vertex-based slope limiter presented in [START_REF] Kuzmin | A vertex-based hierarchical slope limiter for p-adaptative discontinuous Galerkin methods[END_REF], which allows us to preserve smooth extrema. In the end, one can obtain the limiting coefficients for the Riemann invariants polynomials. Hence, using definitions (129), ( 130) and (131), we recover the limiting coefficients corresponding to the system of variables polynomial approximations

( 1 ρ ) c k = 1 2 (α c +,k + α c -,k ), U c k = 1 2 Z c (α c +,k -α c -,k )n + α c 0,k t, E c k = α c E,k + 1 2 Z c (α c +,k -α c -,k )U c 0 • n + α c 0,k U c 0 • t - 1 2 P c 0 (α c +,k + α c -,k ).
Concerning the unit vector n and t, we could define these projection vectors as the two orthogonal directions e x and e y . This option in the choice of n and t would not ensure symmetry preservation in the case of radial flows on polar grids. Finally, we decide to use the velocity mean value direction over the cell and its orthogonal vector, i.e., n = U c 0 \||U c 0 || and t = e z × U c 0 \||U c 0 ||. In the next section devoted to the numerical results, we shall show that this limiting procedure ensures the preservation of the cylindrical symmetry. The scheme robustness and accuracy will also be assessed using several relevant test cases.

Numerical results

To demonstrate the accuracy and the robustness of our scheme on the gas dynamics system, we have run test cases taken from the literature. During the whole calculation we are working on the fixed initial grid. However, plotting final solutions on the initial mesh, the results are difficult to analyze. Luckily, knowing the deformation gradient tensor everywhere and at anytime, we are able to observe the solution on the actual, deformed, mesh. For a better understanding of the results, all the problem solutions are displayed on the final mesh. Concerning the time discretization, the choice has been made to use the well-known TVD Runge-Kutta, see [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF]. In the remainder, the two particular cases of second and third order of accuracy are assessed. The Noh problem [START_REF] Noh | Errors for calculations of strong shocks using artificial viscosity and an artificial heat flux[END_REF] is a well known test case used to validate Lagrangian schemes in the regime of infinite strength shock wave. In this test case, a cold gas with unit density is given an initial inward radial velocity of magnitude 1. The initial pressure is given by P 0 = 10 -6 and the polytropic index is equal to 5 3 . A diverging cylindrical shock wave is generated which propagates at speed D = 1 3 . The density plateau behind the shock wave reaches the value 16. The initial computational domain is defined by (X, Y ) = [0, 1] × [0, 1]. The boundary conditions on the X and Y axis are wall boundary conditions whereas a pressure given by P ⋆ = P 0 is prescribed at X = Y = 1. We run the Noh problem on a 50 × 50 Cartesian grid. This configuration leads to a severe test case since the mesh is not aligned with the flow. We note that we have a very smooth and cylindrical solution, and that the shock is located at a circle whose radius is approximately 0.2, refer Figure 7(a). On Figure 7(b), we observe that the second-order plot is very sharp at the shock wave front and very similar to the one-dimensional cylindrical solution. Moreover the density at the shock plateau is not far from the analytical value. This shows the ability of our scheme to preserve the radial symmetry of the flow.

Sedov point blast problem.

We consider the Sedov problem for a point-blast in a uniform medium. An exact solution based on self-similarity arguments is available, see for instance [START_REF] Kamm | On efficient generation of numerically robust Sedov solutions[END_REF]. The initial conditions are characterized by (ρ 0 , P 0 , U 0 ) = (1, 10 -6 , 0), and the polytropic index is equal to 7 5 . We set an initial deltafunction energy source at the origin prescribing the pressure in the cell containing the origin as follows, P or = (γ -1)ρ or ε 0 vor , where v or denotes the volume of the cell containing the origin and ε 0 is the total amount of release energy. By choosing ε 0 = 0.244816, as suggested in [START_REF] Kamm | On efficient generation of numerically robust Sedov solutions[END_REF], the solution consists of a diverging infinite strength shock wave whose front is located at radius r = 1 at t = 1, with a peak density reaching 6. First, we run Sedov problem with the second-order DG scheme with a 30 × 30 Cartesian grid on the domain (X, Y ) = [0, 1.2] × [0, 1.2], refer to Figure 8(a). Then, keeping the same conditions, we make use of a first unstructured grid made of 1110 triangular cells, refer to Figure 9(a), and a second unstructured grid composed of 775 polygonal cells produced by a Voronoi tessellation, refer to Figure 9(b). We point out that the triangular grid is made of completely anisotropic elements. In these three cases, the numerical solution is very close to the one-dimensional analytical solution, refer to Figure 8(b), 10(b) and 11(b). At the end of the computation, the shock wave front is correctly located and is almost cylindrical, refer to Figure 8(a), 10(a) and 11(a). Further, the density peak almost reaches 6. These results demonstrate the robustness and the accuracy of our DG scheme.

Taylor-Green vortex problem

To assess the accuracy of our DG scheme we perform a convergence analysis utilizing the smooth solution of the Taylor-Green vortex test case, initially proposed by [START_REF] Dobrev | High Order Curvilinear Finite Elements for Lagrangian Hydrodynamics[END_REF]. This test case is derived by considering an analytical solution of the incompressible Navier-Stokes equations. Here, we use an adapted version to handle the compressible Euler equations. The Taylor-Green vortex is characterized by the following conditions. The computational domain is defined by (X, Y ) = [0, 1]×[0, 1]. The initial density is uniform and denoted by ρ 0 . The initial velocity field is divergencefree and reads U 0 = U 0 sin(πx) cos(πy)

-cos(πx) sin(πy) . Let us point that with these two fields, the volume equation is automatically satisfied. Now, to satisfy momentum equation, the pressure field, P 0 , is obtained by balancing the inertia term Substituting the velocity expression in the above equation leads to

∇ x P 0 = -ρ 0 d U 0 dt .
P 0 = 1 4 ρ 0 (U 0 ) 2 [cos(2πx) + cos(2πy)] + C 0 , (132)
where C 0 is a constant that allows to define a non-negative pressure. Using these definitions of the density, velocity and pressure, the volume equation and the momentum equation are automatically satisfied. However, since we are computing this solution by solving the compressible Euler equations, it remains to check whether or not the total energy equation is satisfied. To this end, let us investigate the time rate of change of the specific internal energy ε

ρ 0 d ε 0 dt + P 0 ∇ x U 0 = 1 γ -1 U 0 • ∇ x P 0 ,
since ε 0 = P 0 (γ-1)ρ 0 and ∇ x U 0 = 0. Replacing U 0 and P 0 by their analytical expressions, we finally obtain

ρ 0 d ε 0 dt + P 0 ∇ x U 0 = π 4 ρ 0 (U 0 ) 3 γ -1 [cos(3πx) cos(πy) -cos(3πy) cos(πx)]. (133)
Finally, adding the kinematic energy equation to the previous equation leads to the total energy equation

ρ 0 d E 0 dt + ∇ x (P 0 U 0 ) = π 4 ρ 0 (U 0 ) 3 γ -1 [cos(3πx) cos(πy) -cos(3πy) cos(πx)]. (134)
The numerical simulation of this test case solving the Lagrangian hydrodynamics equations requires the addition of the above source term in the energy (total or internal) equation. For the numerical applications, we set ρ 0 = 1, C 0 = 1, U 0 = 1 and γ = 7 5 .

To illustrate the accuracy of the second-order scheme, we first compare the numerical results obtained to the exact solution, at the final time t = 0.75, refer to Figure 12. Secondly, using the analytical solution we compute the global truncation errors, refer to (a) Second-order scheme without limitation. in Figure 12 exhibit a satisfying accuracy. One can also observe how the straight line edges assumption restricts the grid motion. Concerning the rates of convergence gathered in Table 1, the results

L 1 L 2 L ∞ h E h L 1 q h L 1 E h L 2 q h L 2 E h L∞ q h L∞ 1 10
5.06E-3 1.94 6.16E-3 1.93 2.20E-2 1.84

1 20
1.32E-3 1.98 1.62E-3 1.97 5.91E-3 1.95 confirm the expected second-order rate of convergence. Now, to overcome the grid motion restriction observed in Figure 12, we shall present the numerical results obtained with the third-order extension of our scheme.

4.2. Third-order scheme 4.2.1. Polar Sod shock tube problem.

We consider the extension of the classical Sod shock tube [START_REF] Sod | A survey of several finite difference methods for systems of non-linear hyperbolic conservation laws[END_REF] to the case of the cylindrical geometry. This problem consists of a cylindrical shock tube of unity radius. The interface is located at r = 0.5. At the initial time, the states on the left and on the right sides of the interface are constant. The left state is a high pressure fluid characterized by (ρ 0 L , P 0 L , U 0 L ) = (1, 1, 0), the right state is a low pressure fluid defined by (ρ 0 R , P 0 R , U 0 R ) = (0.125, 0.1, 0). The gamma gas law is defined by γ = 7 5 . The computational domain is defined in polar coordinates by (r, θ) ∈ [0, 1] × [0, π 2 ]. We prescribe symmetry boundary conditions at the boundaries θ = 0 and θ = π 2 , and a wall boundary condition, i.e., the normal velocity is set to zero, at r = 1. The aim of this test case is to assess the scheme accuracy and its ability to preserve the radial symmetry. Thanks to Section 3.6, we start the computation with an initial curvilinear grid such as the ones displayed in Figure 13. In the first case, the domain defined in polar coordinates by (r, θ)

∈ [0, 1] × [0, π 2
] is made of 100 cells in the radial direction and 3 non-uniform cells in the angular direction, characterized by an angle of Π 20 , 7Π 60 and Π 3 , refer to Figure 13(a). This non-uniformity of the grid could cause serious problems in the preservation of the symmetry. Indeed, the results obtained with the first and second-order schemes, on the initial mesh displayed in Figure 13(a), clearly exhibit this loss of symmetry, refer to Figure 14. Nonetheless, running the same problem with the third-order scheme, on the same initial mesh, one gets excellent results, refer to Figure 15. In Figure 15(b), we have plotted the density in all cells versus the cell center radius. We observe that the numerical solution is almost superimposed on the analytical one. This shows the great accuracy of our third-order scheme and its ability to preserve cylindrical symmetry. In this case, the symmetry preservation is due to the high accuracy of the scheme. We can thus conclude that dealing with curvilinear grids, a high-order numerical scheme is required to preserve symmetry, without a specific treatment. Now, we present by Bezier curves, the grid is not perfectly circular. Nevertheless, at the end of the computation, the numerical solution performed by the third-order DG scheme exhibits a very good symmetry preservation, refer to Figure 16(a). We also note on Figure 16(b) that the numerical solution is very close to the one-dimensional cylindrical solution. This result proves the strong accuracy and robustness of the third-order scheme without limitation in this particular case of a grid made of only one cell in the radial direction.

Variant of the Gresho vortex problem

In this section we propose a variant of the initial vortex problem defined by Gresho in [START_REF] Gresho | On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via finite-element method that also introduces a nearly consistent mass matrix[END_REF]. The analytical solution is obtained by considering a steady vortex which is solution of the incompressible Euler equations. This solution is characterized by a balance between inertia and pressure gradient into the momentum equation. We recall hereafter the procedure to derive such a solution. Before proceeding any further we start by recalling the writing of the differential operators using cylindrical polar coordinates. Let (x, y, z) denote the usual Cartesian coordinates. The polar cylindrical coordinates are (r, θ, z) where r ≥ 0 and θ ∈ [0, 2π[, and we have the following relationships x = r cos θ, y = r sin θ.

Let (e x , e y , e z ) be the Cartesian orthonormal basis of the three dimensional space R 3 . The corresponding basis in cylindrical polar geometry is (e r , e θ , e z ), where e r = cos θ e x + sin θ e y , e θ = -sin θ e x + cos θ e y .

For P = P (r, θ, z) and U = u r e r + u θ e θ + u z e z , the gradient and the divergence operators are expressed as follows in cylindrical polar coordinates

∇P = ∂ P ∂r e r + 1 r ∂ P ∂θ e θ + ∂ P ∂z e z , ∇ U = 1 r ∂ (ru r ) ∂r + 1 r ∂ u θ ∂θ + ∂ u z ∂z .
The curl of the velocity field U is given by

∇ × U = ( 1 r ∂ u z ∂θ - ∂ u θ ∂z )e r + ( ∂ u r ∂z - ∂ u z ∂r )e θ + 1 r ( ∂ (ru θ ) ∂r - ∂ u r ∂θ )e z . (135)
Let ψ = ψ(r, θ, z, t) be a physical variable attached to the fluid, its material derivative writes as

d ψ dt = ∂ ψ ∂t + U • ∇ψ = ∂ ψ ∂t + u r ∂ ψ ∂r + u θ r ∂ ψ ∂θ + u z ∂ ψ ∂z .
Being given the velocity vector U , the acceleration reads as

d U dt = ∂ U ∂t + u r ∂ ∂r + u θ r ∂ ∂θ + u z ∂ ∂z U . (136)
Let us assume that the fluid flow obeys the following properties:

• The velocity field is orthoradial and its orthoradial component depends only on the radius r, that is,

U = u θ (r)e θ . (137)
• The thermodynamics variables, ρ and P are scalar valued functions with respect to r.

• The thermodynamic closure is provided by a gamma gas law, P = (γ -1)ρε.

Using the previous assumptions, we readily obtain

∇ρ = d ρ dr e r , ∇P = d P dr e r , ∇ε = d ε dr e r . (138)
We also point out that the flow is incompressible since from (137) we deduce that ∇ U = 0. Moreover, the material derivative of a physical variable such that ψ = ψ(r, t) coincides with its time derivative since the velocity is orthogonal to the gradient of this variable

d ψ dt = ∂ ψ ∂t .
These points involve that the internal energy and the mass continuity equations are automatically satisfied. Replacing the velocity field into (136) leads to the following expression for the acceleration

d U dt = - u 2 θ r e r .
In writing this equation, we have used the fact that d e θ dθ = -e r . Finally the momentum equation reduces to the ordinary differential equation

d P dr = -ρ u 2 θ r . (139)
The time evolution of the specific kinetic energy, k = 1 2 U 2 , reads as

ρ d k dt + ∇P • U = 0. ( 140
)
Knowing that k = k(r) and that ∇P and U are orthogonal leads to the conservation of the kinetic energy. Finally, using (135) and the previous assumptions, the curl of the velocity field reduces to

∇ × U = ( u θ r + d u θ dr )e z . (141)
Being given ρ = ρ(r) and u θ = u θ (r), the pressure field is obtained by solving the differential equation ( 139). In what follows, we propose to derive a smooth vortex solution by setting ρ = ρ 0 and

u θ (r) = 2 2n U 0 ( r rv ) n (1 -r rv ) n if r ∈ [0, r v ], 0 if r ∈]r v , 1]. (142)
Here, r v denotes the radius of the vortex and n is a integer such that n > 1. The factor 2 2n is a normalization factor chosen such that the maximum of the orthoradial velocity is equal to one. Choosing the above form for the orthoradial velocity ensures that it is a smooth function with a compact support over [0, 1]. The integration of the momentum equation (139) gives

P (r) = P (0) + 2 4n ρ 0 U 2 0 h( r rv ) if r ∈ [0, r v ], P (0) + 2 4n ρ 0 U 2 0 h(1) if r ∈]r v , 1], (143)
where P (0) is an arbitrary integration constant and h is the real valued function defined by

h(ξ) = ξ 0 s 2n-1 (1 -s) 2n ds, for ξ ∈ [0, 1]. (144)
The vorticity, ω = |∇ × U | is determined by

ω = U 0 r v g(ξ) ξ + dg dξ , (145)
where ξ = r rv and g(ξ For numerical applications, we define the analytical solution setting the parameters to the following values: U 0 = 1, ρ 0 = 1, P (0) = 5, r v = 0.4 and n = 6. We have displayed the corresponding curves in Figure 17. This variant of the Gresho problem is an interesting validation test case to assess the robustness and the accuracy of a Lagrangian scheme. On the one hand, the vorticity leads to a strong mesh rotation which can cause some problems such as negative Jacobian determinants or negative densities. On the other hand, if the numerical diffusion is too important the flow motion stops very early. We run this Gresho problem on a polar grid made of 40 × 18 cells with our first, second, and third-order DG schemes. The computational domain is defined in polar coordinates by (r, θ) ∈ [0, 1] × [0, 2π]. The solutions are displayed with a zoom in the zone (r, θ) ∈ [0, 0.5] × [0, 2π]. We start with the first-order scheme. In Figure 18(a), we observe that the grid is barely deformed. This is due to the too important numerical diffusion. Obviously, the grid being slightly deformed, the mesh does not present any tangled cells, neither than crossed points nor negative Jacobians in the triangular cells. Now, with the second-order DG scheme, the solution presents another problem. This time, the computation does not stop before the final time. At the end of the computation, the grid is strongly deformed, and one can see the vortex structure inherent to this Gresho problem, refer to Figure 18(b). However, the linear approximation as well as the straightline representation do not allow to follow properly the deformation. Indeed, on Figure 18(b) we note that the final grid is characterized by non-valid cells, wherein some triangles exhibit negative Jacobian determinant. Finally, we perform this Gresho vortex problem with our third-order DG scheme on a curved polar grid made of 40 × 18 as well. We also display in Figure 19(b) the exact motion of this mesh to compare it with our numerical solutions. First, we note that as in the first-order and second-order cases, the third-order scheme preserves symmetry. This scheme is characterized by a very low level of numerical dissipation and the grid is extremely deformed at the end of the computation, more than in the second-order case, refer to Figure 19(a). Thanks to the scheme properties and to the Bezier representation, the solution is very close to the expected one, refer to Figure 19(b). Furthermore, the grid does not contain any non-valid cells, and the Jacobian of the triangular cells remain positive. We have also displayed the plots corresponding to the pressure, the velocity and the density expressed at the centroid of the cells, in the three different cases presented, refer to Figure 20. We note on Figure 20(a) and Figure 20(b) that the numerical solutions corresponding to the first-order schemes are totally smeared. We also observe the huge gain in accuracy between the second-order and the third-order DG scheme. Despite the extreme deformation of the grid, the numerical solutions obtained by the third-order scheme are very close to the exact solutions. This confirm the high accuracy and robustness of the designed scheme. Now, regarding the incompressibility assumption, this test case deriving from a solution of the incompressible Euler equations, the density must remain equal to one during the calculation. We note on Figure 20(c) that the result obtained with the third-order scheme is more accurate than the ones obtained with the first-order and second-order numerical schemes. At the end of the computation, even if the mesh is highly deformed, the incompressibility assumption is very well

) = ξ n (1 -ξ) n if ξ ∈ [0, 1] else g(ξ) = 0. For n = 6,

Sedov point blast problem

Once more we consider the Sedov problem which has been already described. We run this problem, with our third-order DG scheme employing the limitation procedure based on the Riemann invariants polynomial approximation, described in the Section 3.7. The initial computational domain is defined by (X, Y ) = [0, 1.2] × [0, 1.2] and paved by a 30 × 30 Cartesian grid. We note on Figure 21(b) that the numerical solution is very close to the one-dimensional analytical solution. Further, we observe that the shock wave front is cylindrical and well located at the end of the computation, refer to Figure 21(a). These results also demonstrate the robustness and the accuracy of this scheme. Nonetheless, overlapping cells are visible on the grid at the stopping time, refer to Figure 21(a). Further improvements are certainly needed to cure this weakness of the scheme. This cells overlapping phenomenon may result from several reasons. The first one is the local criteria of our discretization. The scheme has been built such as the normals in the actual configuration are continuous on cells sharing a common face. But there is no continuity in the normals nor in the tangents between neighboring edges. The choice of Bezier curves to parametrize the mesh edges has been done to be consistent with the definition of the mapping using P2 finite elements basis functions, these curves being the trace of the P2 polynomials functions on the edges. It is clear that in some cases, this representation is not accurate enough, and thus bring geometric discontinuities. Considering the polar mesh depicted in Figure 13, the Bezier representation does not allow us to obtain a circular mesh. The normals are discontinuous between edges in the angular direction. In [START_REF] Boutin | Extension of ALE methodology to unstructured conical meshes[END_REF], the authors use conics to parametrize their cell edges. These functions allows a perfect representation of polar grids, and the normals would be continuous in this case. However, these functions do not correspond to the trace of the P2 finite elements basis functions on the edges, another discretization of the mapping and thus of the deformation gradient tensor would be needed. Nevertheless, we think that this phenomenon may likely derive from the limitation procedure. The loss of accuracy due to the limitation may downgrade the approximation of the fluid flow velocity, and therefore the deformation of the edges. Another possible explanation of this phenomenon may come from the face control point velocity definition. The use of the Hadamard compatibility condition [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing schemes[END_REF] and helped by the fact that the deformation gradient tensor is discontinuous between triangles inside the polygonal cells could bring more diffusion and thus more stability in the definition of the interior points velocity.

Taylor-Green vortex problem

We make use of the smooth Taylor-Green vortex test case described in the previous second-order section to assess the accuracy of the third-order scheme. First, we compare the solution obtained with the third-order discontinuous Galerkin scheme with the exact solution, on a Cartesian grid made of 10×10 cells, see Figure 22. The results displayed in Figure 22 accuracy compared to the second-order scheme, see Figure 12, the numerical solution corresponding to the third-order scheme being very close to the exact solution. Now, knowing the analytical solution of this Taylor-Green problem, we compute the global truncation error corresponding to our third-order DG schemes and display it in Table 2. Comparing the results gathered in Tables 1 and 2, we conclude that as expected, the third-order scheme is a lot more accurate than the secondorder scheme. We also notice in Table 2 that the asymptotic regime is reached at 3, which proves as expected the third-order accuracy.

L 1 L 2 L ∞ h E h L 1 q h L 1 E h L 2 q h L 2 E h L∞ q h L∞ 1 10

Computational efficiency

Finally, we study the efficiency of our numerical method by comparing the first, second and thirdorder DG schemes, keeping approximately constant the number of degrees of freedom (DOF) in the DG discretization, i.e. the number of moments required in the polynomial approximation of the unknowns. To assess the efficiency, we have displayed the numerical error and the computational time. The corresponding test case is the Taylor-Green problem. We see with the results gathered that for approximately 600 degrees of freedom, the global truncation errors are divided by 10 from the first-order to the second-order, and again from the second-order to the third-order scheme. We also see that even if the computational time required by the third-order scheme is greater than for the lower order ones, the evaluated times remain small. In the case of 2400 DOF, the error is now divided by 20 increasing the order of accuracy. We also observe that the computational time of the first and second-order schemes are very close, while the measured time in the third-order case is three times greater. It is worth mentioning that the scheme is not parallelized and not optimized at all. We can hope with simple improvements in the code to regain an equivalent computational time to the lower orders. Anyhow, for a three times greater computational time, the third-order scheme is 20 times more accurate than the second-order scheme, and 400 times more accurate than the first-order numerical method.

Conclusion

We have presented a cell-centered high-order DG discretization devoted to the numerical simulation over general unstructured grids of the compressible Euler equations written under the total Lagrangian formulation. In this work, the GCL is discretized by means of a finite element approximation, which fulfills the Piola compatibility condition, whereas the physical conservation laws are approximated employing a local variational formulation and an expansion of the physical variables over a Taylor basis. The representation of the flow in the current configuration is performed by means of a curvilinear grid whose motion is rigorously consistent with the GCL. The numerical method ensures momentum and total energy conservation and satisfies an entropy inequality at the

Figure 1 :

 1 Figure 1: Notations for the Rankine Hugoniot relations.

2 and

 2 ∆Y c = Ymax-Y min 2 are scaling factors where X max / min = max / min p∈P(c) {X p } and Y max / min = max / min p∈P(c) {Y p }, defining P(c) as the vertices set of the cell Ω c . The purpose of these scaling factors is to improve the condition number of the mass matrix. Thanks to this particular basis function definition, it easily follows that (σ c 0

Figure 2 :

 2 Figure 2: Triangular subdivision of a generic polygonal cell Ωc.

Figure 3 :

 3 Figure 3: Nodes arrangement for a cubic Lagrange finite element mapping.

Figure 4 :

 4 Figure 4: Different control points contribution of the numerical flux in the cell Ωc.

Figure 5 :

 5 Figure 5: Transformation of Ωc to ωc through the flow. • boundary cell nodes, boundary cell control points, interior cell points.

4. 1 .

 1 Second-order scheme 4.1.1. Noh problem.

  Density in all cells versus cell center radius.

Figure 7 :

 7 Figure 7: Second-order DG scheme with limitation for the Noh problem on a 50 × 50 Cartesian grid.

  Density in all cells versus cell center radius.

Figure 8 :

 8 Figure 8: Second-order DG scheme with limitation; Numerical results for the Sedov problem on a 30 × 30 Cartesian grid.

Figure 9 :

 9 Figure 9: Unstructured grids for computing the Sedov problem.

Figure 10 :

 10 Figure 10: Second-order DG scheme with limitation; Numerical results for the Sedov problem on a triangular grid made of 1110 cells.

  Density in all cells versus cell center radius.

Figure 11 :

 11 Figure 11: Second-order DG scheme with limitation; Numerical results for the Sedov problem on a polygonal grid made of 775 polygonal cells.

Figure 12 :

 12 Figure 12: Taylor-Green vortex deformation of a Cartesian grid made of 10 × 10 cells, at time t = 0.75.

  100 × 1 cells.

Figure 13 :

 13 Figure 13: Initial curvilinear polar grids defined in polar coordinates by (r, θ) ∈ [0, 1] × [0, π 2 ].

  Second-order scheme.

Figure 14 :

 14 Figure 14: Polar Sod shock tube problem. Density maps obtained with the first and second-order DG schemes on the domain defined in polar coordinates by (r, θ) ∈ [0, 1] × [0, π 2 ] made of 100 × 3 non-uniform cells.

Figure 15 :

 15 Figure 15: Polar Sod shock tube problem. Solution obtained with the third-order DG scheme without limitation on the domain defined in polar coordinates by (r, θ) ∈ [0, 1] × [0, π 2 ] made of 100 × 3 non-uniform cells.

Figure 16 :

 16 Figure 16: Polar Sod shock tube problem. Solution obtained with the third-order DG scheme without limitation on the domain defined in polar coordinates by (r, θ) ∈ [0, 1] × [0, π 2 ] made of 100 × 1 cells.

Figure 17 :

 17 Figure 17: Smooth vortex analytical solution corresponding to n = 6.

  Second-order scheme without limitation.

Figure 18 :

 18 Figure 18: Gresho problem variant on a polar grid made of 40 × 18 cells at the final time t = 1. Grid deformations obtained with the first and second order scheme.

Figure 19 :

 19 Figure 19: Gresho problem variant on a polar grid made of 40 × 18 cells at the final time t = 1. Grid deformations obtained with the third-order scheme without limitation and the analytical solution.

Figure 20 :

 20 Figure 20: Gresho problem variant on a polar grid made of 40 × 18 cells at the final time t = 1: comparison between the analytical solution and those obtained with the first, second and third-order numerical schemes, plotted at the centroid of the cells.

Figure 21 :

 21 Figure 21: Point blast Sedov problem on a Cartesian grid made of 30 × 30 cells with the third-order scheme with limitation.

  clearly show the huge gain in 51

Figure 22 :

 22 Figure 22: Taylor-Green vortex deformation of a Cartesian grid made of 10 × 10 cells, at time t = 0.75.

Table 1 .

 1 The results displayed

	1
	0.9
	0.8
	0.7
	0.6
	0.5
	0.4
	0.3
	0.2
	0.1
	0
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Table 1 :

 1 Rate of convergence computed on the pressure in the case of the Taylor-Green vortex at time t = 0.1, with the second-order DG scheme without limitation.

	1 40 1 80 1 160	3.33E-4 1.99 4.12E-4 1.99 1.53E-3 1.98 8.35E-5 2.00 1.04E-4 2.00 3.86E-4 1.99 2.09E-5 -2.60E-5 -9.69E-5 -

  the function h is given by

	h(ξ) =	1 24	ξ 24 -	12 23	ξ 23 + 3ξ 22 -	220 21	ξ 21 +	99 4	ξ 20 -	792 19	ξ 19 +	154 3	ξ 18
		-	792 17	ξ 17 +	495 16	ξ 16 -	44 3	ξ 15 +	33 7	ξ 14 -	12 13	ξ 13 +	1 12	ξ 12 .

Table 2 :

 2 Rate of convergence computed on the pressure in the case of the Taylor-Green vortex at time t = 0.1, with the third-order DG scheme without limitation.

		2.67E-4 2.96 3.36E-7 2.94 1.21E-3 2.86
	1 20 1 40 1 80 1 160	3.43E-5 2.97 4.36E-5 2.96 1.66E-4 2.93 4.37E-6 2.99 5.59E-6 2.98 2.18E-5 2.96 5.50E-7 2.99 7.06E-7 2.99 2.80E-6 2.99 6.91E-8 -8.87E-8 -3.53E-7 -

Table 3 :

 3 First-order DG scheme at time t = 0.1.

	D.O.F	N	E h L 1	E h L 2	E h L∞	time (sec)
	600	24 × 25 2.67E-2 3.31E-2 8.55E-2	2.01
	2400 48 × 50 1.36E-2 1.69E-2 4.37E-2	11.0
	D.O.F	N	E h L 1	E h L 2	E h L∞	time (sec)
	630	14 × 15 2.76E-3 3.33E-3 1.07E-2	2.77
	2436 28 × 29 7.52E-4 9.02E-4 2.73E-3	11.3

Table 4 :

 4 Second-order DG scheme without limitation at time t = 0.1.

	D.O.F	N	E h L 1	E h L 2	E h L∞	time (sec)
	600	10 × 10 2.67E-4 3.36E-4 1.21E-3	4.00
	2400 20 × 20 3.43E-5 4.36E-5 1.66E-4	30.6

Table 5 :

 5 Third-order DG scheme without limitation at time t = 0.1.

Acknowledgments

R. Abgrall has been funded by the European Research Council (ERC) Advance Grant "ADDECCO" (contract # 226316). This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the Investments for the future Programme IdEx Bordeaux (ANR-10-IDEX-03-02)

semi-discrete level which guarantees its consistency with the Second Law of thermodynamics. The numerical results display an accuracy up to third-order for smooth solutions. The gain in accuracy provides a dramatic improvement of the symmetry preservation for symmetric flows. In the future, we intend to improve the formulation of the characteristic variables based limiting procedure. We also plan to develop a DG discretization of the gas dynamics equations written under the Lagrangian updated formulation and to extend its capability to the non-linear elasticity equations.