An Arbitrary Lagrangian-Eulerian strategy to solve compressible fluid flows

Keywords: Lagrangian Fluid Dynamics, ALE, Rezoning, Nodal Mesh Quality, Non linear local mesh relaxation, Conservative Local Remapping, Discrete Maximum Principle (local bound preserving) enforcement

In this paper, we consider an Arbitrary Lagrangian Eulerian (ALE) alternative to the computation of bi-dimensional compressible fluid dynamics written in Lagrangian frame. In this formulation, the conservation laws are solved on a mesh moving at the speed of the flow. In the context of multi-material flows, this Lagrangian description allows the exact preservation of the interface between different species.

.

Recently, new schemes have been designed as a full centered version [7] [22]. These two schemes can be seen as an extension of Godunov scheme to fully multidimensional flux (at nodes). Each one verifies a local discrete vertex centered entropy inequality.

Despite all efforts to obtain robust simulations for Inertial Confinement Fusion (ICF) computations, nowadays, the pure Lagrangian frame mode does not ensure a good cell quality during the calculation: for example non convex cell and/or tangled cell may appear.

. A pure Lagrangian phase is followed by a two step process: the rezoning of the Lagrangian mesh followed by a remapping step.

.

Introduction

The purpose of the paper is to describe numerical methods for the computation of two dimensional compressible Euler flows written in Lagrangian framework. Lagrangian schemes was first based on a staggered version (degrees of freedom for velocity are located at nodes and thermodynamics variables are considered as cell centered [START_REF] Wilkins | Calculation of elastic-plastic flow[END_REF], [START_REF] Campbell | A compatible lagrangian hydrodynamics algorithm for unstructured grids[END_REF], [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF]). Recently full cell centered schemes were proposed see [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF] [START_REF] Maire | A cell-centered lagrangian scheme for two-dimensional compressible flow problems[END_REF]. We consider an Arbitrary-Lagrangian-Eulerian (ALE) method formulated as a splitting between the pure Lagrangian step and the association of rezoning [START_REF] Winslow | Numerical solution of the quasilinear poisson equation in a non-uniform triangular mesh[END_REF] [29] [START_REF] Jun | A modified equipotential method for grid relaxation[END_REF] [17] [START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF] plus a remapping step [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF]. We will propose new methods for each part of this second step. First, for the rezoning process we describe a strategy in order to obtain a smoothed mesh close to Lagrangian grid. Roughly speaking, the method is derived by the definition of a nodal mesh quality. We will also explain three main consequences of this definition, such as non linear mesh relaxation, global mesh quality measurement, and the possible mixing of different smoothing grid algorithms. Second of all, for the remapping step, we introduce an approximate version of polygon-polygon intersection that computes exact intersection between two arbitrary meshes. Our approach is intermediate between a pure exact intersection algorithm see [START_REF] Weiler | Hidden surface removal using polygon area sorting[END_REF] [START_REF] Greiner | Efficient clipping of arbitrary polygons[END_REF] and the local volume swept region method in [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF], [START_REF] Loubère | A subcell remapping method on staggered polygonal grids for arbitrary-lagrangian eulerian methods[END_REF] where no intersection is needed. In our case, only a very simple segment-segment intersection algorithm is needed, offering still a linear complexity algorithm (w.r.t the number of cells). Finally, we discuss for first and second order accuracy, the local maximum principle for density and any specific quantities. The paper is organized as follows. First, we briefly recall the pure Lagrangian schemes for gas dynamics. Second, we deal with a commonly used alternative to pure Lagrangian scheme : the non direct ALE formalism constituted by two steps. Concerning the rezoning step, we introduce the definition of nodal mesh quality. As a consequence, we can naturally define non linear mesh relaxation, and adaptive local and mixed smoothing mesh algorithms. These tools help to control on one hand the quality of the final smoothed mesh and on the other hand its "proximity" to the Lagrangian grid. Some numerical examples are given.

To conclude the rezoning part, we give also a variant of Escobar et al [START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF] for the smoothing and untangling of meshes with arbitrary connectivity. In our case, we take into account the local mesh connectivity around the node, the degree of the node and each individual cell around this node. In the fourth part, we recall the remapping step needed to recompute the new values in the cells (and/or nodes) of the new smoothed mesh. The exact problem is based on the intersection of the Lagrangian and rezoned grids. A solution to circumvent the computation of polygon-polygon intersection, is to compute the "swept area" see in [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF], which is an intersection free process.

Here, we propose a modification based on "self" line clipping. For the first order, we obtain a positive scheme under a condition on the node displacement, and guarantee a discrete local maximum principle (local bound preserving). For the "second order" version of our local approximate remapper, we use an iterative process to enforce maximum principle to extensive quantities (such as density). For intensive quantities such as internal energy, see [START_REF] Vanderheyden | Compatible fluxes for van leer advection[END_REF], we use a non linear representation. Some numerical examples using these tools are then presented.

Lagrangian Hydrodynamics

We are mainly interested in numerical approximation of Euler equations. In Eulerian frame, these are the conservation of mass, momentum and total energy:

             ∂ t ρ + ∇.(ρU) = 0, ∂ t (ρU) + ∇.(ρU ⊗ U) + ∇P = 0, ∂ t (ρE) + ∇.((ρE + P )U) = 0. (1)
where ρ, U, E represents the density, velocity, and specific total energy. P stands for the pressure given by a thermodynamical relation P=P (ρ, ǫ), where ǫ is the specific internal energy of the fluid, for a perfect gaz P = (γ -1)ρǫ. To obtain equations on a moving grid, we note that for any quantity q, dq dt = ∂ t q + U.∇q

In the pseudo Lagrangian frame (spatial derivatives are still those of (1)), the equations are:

             ρ d dt τ -∇.U = 0, ρ d dt U + ∇P = 0, ρ d dt E + ∇.(P U) = 0.
(

) 2
where τ = 1 ρ . A commonly used spatial discretization of this system is done after integrating (2) on a moving cell C(t):

             d dt C(t) ρdx = 0,
At this stage, there are basically three different kinds of spatial discretization: staggered, collocated in cell or at nodes.

Staggered schemes

Two staggered schemes are nowadays in use:

(1) Wilkins [START_REF] Wilkins | Calculation of elastic-plastic flow[END_REF]. This scheme is commonly used by many codes and can be considered as a reference. The main features are: (a) Mixed Q1/P0 (Velocity/Thermodynamics). (b) Internal energy formulation ρ d dt ǫ + P ∇.U = 0. The scheme is not conservative for the total energy. (c) Hourglass correction (non physical oscillation for non-triangular mesh can appear otherwise). (d) Artificial viscosity (to add entropy for capturing shocks in velocity/pressure). (2) Shashkov et al. [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF] [START_REF] Scovazzi | Multi-scale Lagrangian shock hydrodynamics on Q1/P0 finite elements: theoretical framework and twodimensional computations[END_REF].

The scheme makes improvement over the previous scheme. The authors proposed to add more robust tools to correct some defects of Wilkins schemes.

(a) Sub-zonal pressure forces control more precisely the hourglass modes. (b) The scheme is conservative for the total energy. (c) Artificial viscosity have been improved.

Centered schemes

We briefly recall the centered scheme of Després/Mazeran [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF] and Maire et al [START_REF] Maire | A cell-centered lagrangian scheme for two-dimensional compressible flow problems[END_REF]. Both of them are finite volume schemes with the following properties:

(1) Momentum and total energy conservation.

(2) Approximate linear Riemann solver: averaged normal to node for Després/Mazeran and normal to face for Maire. (3) Nodal fluxes (compatibility of mesh movement and continuity equation).

The two schemes differ from the intermediate pressure discretization and exhibit some side effects :

• Després/Mazeran nodal solver (first order and more) depends on the aspect ratio of the mesh while Maire et al does not. • The method of Després/Mazeran can develop hourglass modes (first order and more) and needs to be coupled with some hourglass filter or ALE process.

It is an open question to know if Maire et al scheme possesses such modes. • Maire's nodal solver is very dissipative and needs to be second order for some flows.

Nodal Based Finite element methods

For finite element discretization, see the following approach of Barlow [START_REF] Barlow | A compatible finite element multi-material ALE hydrodynamics algorithm[END_REF], and also recently, G. Scovazzi et al. [START_REF] Scovazzi | Stabilized shock hydrodynamics: I. a lagrangian method[END_REF] and [START_REF] Scovazzi | Stabilized shock hydrodynamics: Ii. design and physical interpretation of the supg operator for lagrangian computations[END_REF] have developed a new approach to Lagrangian shock hydrodynamics based on the variational multi-scale method applied to nodal based finite element discretizations. These methods are conservative for the mass, momentum and total energy, and use an artificial viscosity to improve robustness under shocks conditions.In general way, some drawbacks of pure Lagrangian scheme are:

• interacting shocks can destroy mesh quality and increase the error • the mesh can become very distorded for rarefaction waves from a solid into a light material (vacuum), because high curvature of interfaces and boundaries (with fast variations) may develop.

Shear and vorticity tend to tangle and fold the mesh, and drive the computation to breakdown. A commonly accepted solution in these cases is an ALE approach: it improves the geometrical quality of mesh by a smoothing and remapping process.

Remark 1 For this paper, we will show results with Wilkins Lagrangian scheme but it extends of course to any of these. Because most of mesh complexity comes from cells having more than four nodes, in the following sections, we focus on quadrilaterals, but the approach is valid and described for arbitrary conformal connectivity and may be easily extended to three dimensional case.

Mesh rezoning : Escobar connectivity adaptation algorithm and nodal quality control

The basics of our mesh smoothing approach adapts the Escobar et al algorithm [START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF] to an arbitrary polygonal mesh. The local transformation that leads to a local objective function takes into account the local connectivity of the node and of the cells surrounding it. In a second part, we propose a mesh control monitoring with nodal quality. Mainly, any smoothing algorithm tends to produce square quadrilaterals although the Lagrangian initial mesh may have very large aspect ratio and/or area ratio. In ALE computation, if nodes are systematically displaced, a lot of numerical diffusion implied by remapping step tends to smear shocks, and some physical phenomenon may not be captured. Moreover, we need the smoothed grid to be as close as possible to the Lagrangian grid (at new time level) to minimize mass fluxing and numerical diffusion induced. This means that if this grid is almost optimal, we only need to cure very localized regions. In this way, we propose to give a weight to Lagrangian nodes. In the context of mesh smoothing, up to our knowledge, all the algorithms such as elliptic ones Barycentric, Tipton [START_REF] Tipton | Grid optimization by equipotential relaxation[END_REF], Jun [START_REF] Jun | A modified equipotential method for grid relaxation[END_REF], or optimization based: Reference Jacobian Matrix (RJM) [START_REF] Knupp | Reference jacobian optimization-based rezone strategies for arbitrary lagrangian-eulerian methods[END_REF], feasible set [START_REF] Garimella | Untangling of 2d meshes in ale simulations[END_REF], condition number [START_REF] Knupp | Algebraic mesh quality metrics[END_REF], Escobar [START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF] are based on cell quality definition. This notion describes somehow the regularity of the cell aspect see Knupp [START_REF] Knupp | Algebraic mesh quality metrics[END_REF], Pebay [START_REF] Pebay | Some results about the quality of planar mesh elements[END_REF]. In the context of local mesh modification (h-adaptivity [START_REF] Frey | Mesh Generation Application to Finite Elements[END_REF]) this notion is very natural, nevertheless for constant connectivity (r-adaptivity see for example [START_REF] Frey | Mesh Generation Application to Finite Elements[END_REF]), we propose to define the quality at nodes. It will be justified after setting some basic notation.

Notation for ALE unknowns

In this part, we give the basic notation. The superscript/underscript n indicates the node index.

                  
M prev , Lagrangian mesh at time t, M lag , Lagrangian mesh at time t + ∆t (new time level), M ale , Lagrangian rezoned mesh, M n , defines the position of node n in the mesh M .

(

) 5
Recall the following identities:

     M lag n = M prev n + ∆tU prev n , U prev n
is the velocity of the vertex n,

M ale n = M lag n + ∆ ale n . (6)
The unknown here is ∆ ale n for all nodes in M lag . We need to define some notions for connectivity:

(1) The Cell/Node and Node/Cell view point:

             N (C)
(2) The Nodal/Nodal view point:

     N (n) is the set of nodes belonging to one of C(n), #N (n) : number of nodes inside N (n), (9)
We need to define a more detailed description of different kinds of neighborhood of n in (9): short, medium and long range, see Figure 2.

Short C(n) n = { n ′ ∈ C(n) ; |n C l -n ′ C l | = 1 } = {n n , n m } always two nodes, M edium C(n) n = { n ′ ∈ C(n) ; |n C l -n ′ C l | = 2 } , M edium C(n) n =              {∅} , for simplex, {n d } , (n d = (n nn = n pp ) for quadrilateral cell),
{n nn , n pp } , for pentagonal cell,

Long C(n) n = N (C(n))\Short C(n) n \M edium C(n) n . (10)
For meshes of simplices,

N (n) ≡ Short C(n) n
, and for mesh containing mixed simplices, quadrilaterals and pentagons, only Short Remark 2 (1) The total number of the set of nodes inside the Short range is generally called:

degree(n) (11)
(2) The long range of node n (if it exists) is here neglected because it can not define three consecutive nodes inside the cells.

Most of the time, M lag shows pathologies only in very localized regions of the computational domain, so that, we can decide to smooth only some nodes inside a sub-domain. Moreover, we often need to define three different categories of subset nodes, labelled with status, Previous, Lagrangian or Relax:

(1) Previous set:

P REV = n; M n ∈ M lag ; ∆ ale n := -(M lag n -M prev n ) , (12)
(2) Lagrangian set:

LAG = n; M n ∈ M lag ; ∆ ale n := 0 , (13)
(3) Relax set:

REL = n; M n ∈ M lag ; ∆ ale n is a degree of freedom = {n; n / ∈ {P REV ∪ LAG}} . ( 14
)
For the sets PREV and LAG, the displacement is fixed so that belonging points are not degree of freedom for the rezoning algorithm. Practically, at initial time we set all nodes to LAG and decide to raise a node to REL if at some time, a geometric criteria over angles or area has reached a threshold value. The node will keep this status until the end of the computation (or will be re-analyzed every time step). A global smoothing process is done by an iterative procedure where the maximal number of iteration is small (< 20).

The method used to iterate can be formulated at least in two ways:

• Explicit (Jacobi update): Initialize M 0 n = M LAG n for all n, after for each n ∈ REL, compute M s+1 n with a local smoother with all neighborhood j ∈ N (n) fixed. Pros : Permits to conserve symmetry, does not depend on the numbering of vertices. Cons : For initial poor non simplicial meshes (tangled), the convergence towards untangling is very slow, or even worse that untangling is not obtained.

• Semi-implicit (Gauss-Seidel update): Initialize M 0 n = M LAG n for all n, after for each n ∈ REL, compute M s+1 n with a local smoother with neighborhood j ∈ N updated by previous pass. Pros : For initial poor non simplicial meshes (tangled), the convergence is improved. Cons : Symmetry is lost (if not fully converged), depends on the numbering of vertices.

Escobar Optimization process extended to arbitrary connectivity

Extension of Escobar method

In reference [START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF] Escobar et al. propose a way to define all over the space IR d the functionals found in literature designed to improve the global quality of a mesh. Their approach permits simultaneously to untangle and smooth a tetrahedral mesh. They noticed that existing optimization techniques such as [START_REF] Knupp | Algebraic mesh quality metrics[END_REF] [17] Reference Jacobian Method (RJM) consist in minimizing an objective function. This last is generally a rational function q(x) = n(x)/d(x) for which the denominator d(x) may be zero making it not defined all over the whole space. The main idea is to replace d(x) by [START_REF] Barlow | A compatible finite element multi-material ALE hydrodynamics algorithm[END_REF]. The affine map between T and T R is given by x = x 0 + Au, where A is the Jacobian matrix of the map referenced to node x 0 .

x k = (x 1 k , x 2 k ) ∈ IR 2 , k=0,1,2 and T R the reference simplex u 0 = (0, 0), u 1 = (1, 0),u 2 = (0,
A =    x 1 1 -x 1 0 x 1 2 -x 1 0 x 2 1 -x 2 0 x 2 2 -x 2 0    (15)
The authors considers the equilateral triangle to be an ideal triangle T I whose vertices are v 0 = (0, 0), v 1 = (1, 0) and v 2 = (1/2, √ 3/2). Let v = W u be the linear mapping that transforms T R to T I , in this case, the Jacobian matrix is:

W =    1 1/2 0 √ 3/2    . (16)
Finally, the affine mapping taking T I to T is noted by S = AW -1 . The matrix norm is used: the determinant or trace of S is computed to construct algebraic quality measures of T see [START_REF] Knupp | Algebraic mesh quality metrics[END_REF]. The Frobenius norm of S, given by |S| = tr( t SS) is used to define an algebraic cell quality : q(S) = maximum value is reached for equilateral simplex, and flat triangle has zero quality measure. For each node of the mesh, an objective function can be associated to some simplex in Figure 5, f m = 1/q(S), f m = |Sm||Σm| 2σm , where Σ = σS -1 , with σ = det(S). Unfortunately, it may happen that σ m be 0 (in case of two consecutive flat edges) or becomes negative so that the previous function is not defined. Escobar fixes this problem, by the following modified objective function: fm = |Sm||Σm| σm+ √ σ 2 m +4δ 2 , defined on the whole space. The full Escobar functional F n (x) = N bSimplices m=1 fm is constructed by adding all the simplices contributions. Unfortunately, for arbitrary connectivity meshes, we do not know what is an ideal cell. Moreover, the different simplices (Short or Medium Range) do not possess the same amount of information. For this reason, the application W in ( 16) is only an approximation for general polygons. In fact, for each simplex, we will introduce a fitted parameter for each class, Short and Medium. In terms of linear application, it means that we need to introduce a different matrix (W α ) for each of simplices class, and we will adapt this to the non homogeneous local connectivity around each node see Figure 5. From a generic point of view, we have to deal with the following transformation matrix:

W α =    1 cos(α) 0 sin(α)    , S α = AW α -1 . (17)
Here, we propose an extension of Escobar functional that takes into account non homogeneous connectivity of the mesh.

Proposition 3

The transformation matrix [START_REF] Knupp | Reference jacobian optimization-based rezone strategies for arbitrary lagrangian-eulerian methods[END_REF] is computed in the following way:

(1) For all internal simplices (Short Range), to obtain equilibrium (constant angle) around n, the optimal parameter is given by (see [START_REF] Frey | Mesh Generation Application to Finite Elements[END_REF] for notation):

α i = 2π degree(n) , ∀i = 1, .., degree(n) (18)
(2) For external (diagonal) simplices (Medium Range) belonging to the cell C around n (see notation (7) ( 8)), we look an equilibrium inside C, here parameters are:

α n C l ,1 = α n C l ,2 = 2π degree(n)(#N(C) -2) . ( 19
)
Short Proof.

(1) Short Range simplices. We want to balance all information of simplices Here, C i are the i t h local cell linked to node n, T S i the SHORT range sub-simplex inside C i , and T i,1 , T i,2 the MEDIUM range sub-simplices. The angle α i (resp.

α 2,1 2,2 3,1 3,2 4,1 4,2 1,1 1,2 2,1 2,2 3,1 3,2 4,1 4,2 n C
α i,1 ,α i,2 ) is the angle inside T S i (resp. T i,1 , T i,2 ) from n.
see Figure [START_REF] Caramana | The construction of compatible hydrodynamics algorithms utilizing conservation of total energy[END_REF] so that α i = α * (same angle), with the constraint that degree(n) i=1

α i = 2π. (2) Medium Range simplices. Using the previous result, let us consider a cell for which we want to equilibrate the angle inside the cell (see Figure 6) coming from the node n, here the constraint is

#N (C)-2 j=1 α j = 2π degree(n) .
2π degree(n) Using (17) (18) [START_REF] Lopez | Simultaneous untangling and smoothing of moving and fixed grids[END_REF], we can consider now the local objective function:

f α m = |S α m ||Σ α m | σ α m + (σ α m ) 2 + 4δ 2 (20)
F α n (x) = N bSimplices m=1 f α m (21)
For the numerical approximation of this local functional, we have used the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm, see [START_REF] Mokhtar | Nonlinear programming[END_REF]. The parameter used in BFGS are the stopping criteria |∇F α n (x * )| < 10 -4 , and re-initialization of cumulative approximated hessian direction by opposite gradient if the descent direction is too far from it or any 20 iterations. We re-initialize the line-search any 5 iterations. There is still some improvements to make to determine δ for arbitrary meshes, see [START_REF] Lopez | Simultaneous untangling and smoothing of moving and fixed grids[END_REF] for an answer. Practicaly, we choose to take (after setting (σ min = min σ m )) :

δ =      max( ǫ(ǫ -σ min ), ǫ R |σ min |), if σ min < ǫ, (ǫ = 10 -14 ), (ǫ R = 10 -7 ), 0 otherwise.
In next work, we plan to adapt the algorithm to the anisotropic case. At this point, we have build a method whose goal is to obtain robust smoothing process, and we will show some test cases.

Numerical test cases related to Escobar with connectivity adaptation

In the following, we test the smoothing effects of the functionnal minimization (20) on a variety of static mesh (no coupling with hydrodynamics). In the first test, the initial mesh Figure 7 is composed by four quadrilaterals, and possess poorly positioned nodes. At convergence, the smoothing techniques must produce a mesh made of squares for Test 1 Figure [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF] The three previous tests of the proposed Escobar with connectivity adaptation [START_REF] Loubère | A subcell remapping method on staggered polygonal grids for arbitrary-lagrangian eulerian methods[END_REF] show that the convergence is better than the original Escobar but that it does not preserve aspect ratio (as the original Escobar and also Reference Jacobian Method (RJM) [START_REF] Knupp | Reference jacobian optimization-based rezone strategies for arbitrary lagrangian-eulerian methods[END_REF] and Winslow [START_REF] Winslow | Numerical solution of the quasilinear poisson equation in a non-uniform triangular mesh[END_REF]).

Definition of nodal mesh quality and consequences for rezoning step

We begin this section by general considerations.

• The non-linear optimization algorithms such as the Reference Jacobian Method [START_REF] Knupp | Reference jacobian optimization-based rezone strategies for arbitrary lagrangian-eulerian methods[END_REF], Escobar [START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF], etc, are much more time consuming than elliptic linear ones Tipton [START_REF] Tipton | Grid optimization by equipotential relaxation[END_REF] or Jun [START_REF] Jun | A modified equipotential method for grid relaxation[END_REF], or barycentric (Laplace). Moreover, these schemes does not preserve aspect ratio of the Lagrangian mesh which can be bad in some situation (planar shock on an initial cartesian mesh). • The Lagrangian grid needs to be improved as locally as possible.

• The smoothing process does not ensure that all the cells have been really

upgraded with respect to cell quality itself, only some averaged quality is actually improve. • In case of real application (i.e. coupled with hydrodynamics), we never reach convergence at each time step, but we need to control very closely the quality of the rezoned mesh.

We emphasize that we focus on a process fullfiling the following property:

The smoothed mesh must be as close as possible to the Lagrangian grid [START_REF] Maire | A cell-centered lagrangian scheme for two-dimensional compressible flow problems[END_REF] In pratice, we impose the following constraint on the node displacement:

ClosetoLag(n) := ∆t|U prev, * n | |M lag n M ale n | (23)
                   if 0 ≤ ClosetoLag(n) < 1 then M ale n := M ale n × ClosetoLag(n) else M ale n is accepted. ( 24
)
where U prev, * n is some averaged Lagrangian velocity of node n see [START_REF] Després | The vofire method for multi-component flows on unstructured meshes[END_REF]. After node displacement, we ideally expect that new position tends to produce better surrounding cells, but in practice this is not always true, because of highly non linear process induced by [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF]. Our smoothing process is a discretisation of a non-linear (degenerate) Elliptic operator with compact support. We need to define a vertex quality measure for nodal mesh movement. In this way, the node movement influences simultaneously all the cells around. At this stage, we have to define the properties of this nodal function, and we will give some examples of such functions. Definition 4 A nodal mesh quality function Q n is a measure of the regularity for cells belonging to C(n). It is defined and continuous from IR 2 → [-1,1]. Let's define (see Figure 2):

min sin := min C(n) min(sin(α S n ), sin(α SM np ), sin(α SM nn )), max sin := max C(n) max(sin(α S n ), sin(α SM np ), sin(α SM nn )), min area := min C(n) min(A(T S n ), A(T SM np ), A(T SM nn )), max area := max C(n) max(A(T S n ), A(T SM np ), A(T SM nn )).
Then, we introduce some basic examples of nodal functions:

sin: Q sin n = min sin sin,rel : Q sin,rel n = min sin max(max sin,| min sin |) area,rel : Q area,rel n = min area max(max area,| min area|) aspect ratio weight: Q ar, * n = min C(n) min(Ln p ,Ln n ) max C(n) max(Ln p ,Ln n ) Q r
n , (r is one of the 3 previous functions, sin, sin rel, area rel).

All these nodal quality functions must fullfill the following properties:

(1) Q n ≃ 1 iff the mesh around n is locally optimal, (2) Q n = 0 iff at least a cell around n is locally flat (n is aligned with 2 consecutives j ∈ Short(n)), (3) Q n < 0 iff at least a cell around n is locally degenerate (non convex or tangled).

In our case, we need to define functions with negative values. Indeed for arbitrary polygonal cells, we need to make a difference between convex flat (Figure 11 Case 1 and Case 2) and strictly non convex cell (and a fortiori tangled cell Figure 11 Case 3 and Case 4 for which the mesh is considered as unvalid).

We need to mark degenerated nodes, because ideally, we must maintain cell convexity, at least because the existence of the Jacobian transformation between Lagrangian and Eulerian continuous frames does need it. This previous definition 4 does not depend of Long C(n) n range because it would envolve more that three consecutive nodes, hence the word "locally" in the three properties is only based on Short and Medium neighborhood. The evaluation of these functions is inexpensive because it uses angular and area criteria which have already been computed to determine if a node needs to pass from LAG to REL see [START_REF] Garimella | Untangling of 2d meshes in ale simulations[END_REF] [START_REF] Greiner | Efficient clipping of arbitrary polygons[END_REF].

Remark 5 For nodes that are not degree of freedom in the rezoning process (in PREV or LAG), we decide to give the maximum value Q n := 1 for any function because their displacement is known, so that their new location can be considered as optimal. Remark 6 From the cell quality notion, we can define quality of vertex n by taking the minimum cell quality of all the neighboring cells (which is exactly the same for simplices meshes). Here, we propose an alternative but weaker approach below (with the meaning of Short, Medium and Long range see equation ( 10)). In definition 4, the reader can easily check that the functions verify:

• Q n < 0 ⇒ it exists a non convex cell in C(n). • Q n ≥ 0 ⇒ T S n in Short C(n) n and T SM np , T SM nn in M edium C(n) n are not tangled.
We give some basic properties of previous nodal function definition.

Properties 7 (1) The function Q sin n (see definition 4) is optimal for the following patches: local orthogonality with arbitrary aspect ratio (very inter-esting for symmetry preservation), this function is well suited for quadrilateral (logical connectivity) meshes.

(2) The function Q area,rel n is optimal for the following patches: arbitrary iso area surrounding cells (can preserve symmetry).

(3) The introduction of aspect ratio for the last function is there to control the ratio of the extremum length of adjacent edges. As a first consequence of the previous definition 4, we can introduce a global mesh quality using each individual nodal quality.

Definition 8

Let Q M n be an arbitrary nodal quality function for node n in mesh M, we define two global qualities for mesh M:

(1) Q M ∞ = min n Q M n (2) Q M 1 = 1 Nv Nv n Q M n (N v : Total number of vertices) (25)
Beside this definition, we need to define a parameter to compare two given meshes.

Q ref : global reference quality, Q ref ∈ [0, 1
] (given by the user or given by some estimate see below).

(

) 26
Using Q ref , we can define:

N b M (Q ref ) : the number of nodes such that Q M n < Q ref . (27)
We note that Q M ∞ is linked at least to one vertex (non uniqueness) with the worst nodal quality. Practicaly, we use a real number Q min ∈ (0, 1) such that no rezoning occurs when Q M ∞ ≥ Q min . And we set

Q ref = C 1 Q min , C 1 ≤ 1, (C 1 ≡ 1/2 as instance). But for simplicity we use Q ref ≡ Q min in the following.
We are now able to quantify the whole mesh quality thanks to (25)(26) [START_REF] Scovazzi | Multi-scale Lagrangian shock hydrodynamics on Q1/P0 finite elements: theoretical framework and twodimensional computations[END_REF]. Moreover, we can also compare two meshes M 1 and M 2 .

Definition 9 M 1 is strictly better than M 2 if and only if one of the three situations is verified:

(1) If Q M 1 ∞ ≥ 0 and Q M 2 ∞ ≥ 0, then one of the three following cases: . Q M 1 ∞ > Q M 2 ∞ and N b M 1 (Q ref ) ≤ N b M 2 (Q ref ), OR . Q M 1 ∞ ≥ Q M 2 ∞ and N b M 1 (Q ref ) < N b M 2 (Q ref ), OR . Q M 1 ∞ = Q M 2 ∞ and N b M 1 (Q ref ) = N b M 2 (Q ref ) and Q M 1 1 > Q M 2 1 .
(

) If Q M 1 ∞ > 0 and Q M 2 ∞ < 0. 2
(

) If Q M 1 ∞ < 0 and Q M 2 ∞ < 0, then one of the two following cases: . Q M 1 1 > Q M 2 1 and N b M 1 (Q ref ) ≤ N b M 2 (Q ref ), OR . Q M 1 1 ≥ Q M 2 1 and N b M 1 (Q ref ) < N b M 2 (Q ref ). 3
Remark 10 Another possible definition would be to replace the three inequalities in (1) by the same two inequalities of (3) offering a more symmetric treatment. However, this first definition permits to control the worst quality node. This definition is then a rule for the case where we want to penalize non convex meshes.

A second consequence is a way to generalize the global mesh relaxation process by a local relaxation much more natural and well suited to the proximity requirement.

Definition 11 Let M 1 and M 2 be two meshes with the same connectivity, and Q M n a nodal quality function, we define a third mesh M 3 by : ∀n, n th vertex of M 3 :

M 3 n = M 1 n + w(Q M 1 n , Q M 2 n )(M 2 n -M 1 n ) (28) w(Q M 1 n , Q M 2 n ) is a weight function with values in [0,1]. Examples (1) w(Q M 1 n , Q M 2 n ) = Constant, corresponding to a global mesh relaxation between M 1 and M 2 . (2) w(Q M 1 n , Q M 2 n ) = 1+Q M 2 n 2+Q M 1 n +Q M 2 n ( Q M 1 n +Q M 2 n = -2). This is a kind of average between both, if Q M 1 n and Q M 2 n are equal, then w(Q M 1 n , Q M 2 n ) = 0.5

, and the function is decreasing in

Q M 1 n for Q M 2 n fixed, increasing in Q M 2 n for Q M 1 n fixed.
(3) Heaviside kind of process:

w(Q M 1 n , Q M 2 n ) =      0, if Q M 1 n ≥ Q M 2 n , 1,
else. in this case ( 28) is equivalent to:

M 3 n =      M 1 n , if Q M 1 n ≥ Q M 2 n , M 2 n , else. (29)
this particular mesh is refered as the two "intersected " mesh of M 1 and M 2 .

Properties of non-linear nodal relaxation (28)

(1) Order independent, the construction of M 3 is very fast (only positions differs). ( 2) If M 1 and M 2 are close then M 3 is close to M 1 and M 2 ("between").

(3) We have to be careful because

Q M 1 n > Q M 2 n does not implies Q M 3 n ≥ Q M 1 n > Q M 2
n . We can not deduce a comparison principle for the third mesh. However, in many situations it does, as we will see in numerical examples Figure 14 and Figure 15, and moreover, we can take the best mesh over the three using the previous definition [START_REF] Frey | Mesh Generation Application to Finite Elements[END_REF].

As an illustration of this kind of process, we give some examples of the non linear relaxation below. For the three examples, we choose the Q sin n as nodal function and we choose Q ref = 0. We test equation ( 29) on different static rezoning process (no hydrodynamic coupling). For each one, with a given mesh M1 (the data), we construct a second mesh M2 by using a more or less successful "smoothing" process, and we apply the formula [START_REF] Tipton | Grid optimization by equipotential relaxation[END_REF] to obtain the third mesh M3 call mesh M1 Inter M2.

(1) In the first test Figure 14, mesh M1 gets one tangled cell (see the two nodes just above (0,0)), mesh M2 is obtained using the Jun smoothing algorithm except for some nodes around the origine (0,0). We put away two nodes to tangle the mesh in some other regions than M1. Note that the Jun process has smoothed out the tangled cell from M1. We notice that M 1 ∩ M 2 mesh does not possess any tangled cell anymore, only a slightly non convex cell remains. (2) In a second test see Figure 15, mesh M1 is the same as the previous example, but in this case, the smoothing process has produced a valid and better shaped M2 mesh. Then M 1 ∩ M 2 mesh is quite close to M2.

In practice, we relax the Lagrangian grid using an iterative smoothing process. Therefore, we need to define some parameter such that a node is supposed to

0 1 0 1 Mesh M1 0 1 -1 0 1 0 1 Mesh M2 0 1 -1 0 1 0 1 Mesh M1 inter M2 0 1 Figure 14. First Mesh M1 is characterised by : Q M 1 ∞ =-0.547, N b M 1 (0)=4
, Tangled Cell: 1, the second mesh M2 results from an iteration of Jun's smoothing except for some nodes around the origin (0,0), its characteristics are Q M 2 ∞ =-0.999, N b M 2 (0)=10, Tangled Cells: 2, the third mesh M3 is the intersected mesh resulting from [START_REF] Tipton | Grid optimization by equipotential relaxation[END_REF] with the following characteristic

Q M 3 ∞ =-0.008, N b M 3 (0)=3, Tangled Cell: 0. -1 0 1 0 1 Mesh M1 0 1 -1 0 1 0 1 Mesh M2 0 1 -1 0 1 0 1 Mesh M1 inter M2 0 1 Figure 15.
The first mesh M1 correspond to Mesh M1 in Figure 14 :

Q M 1 ∞ =-0.547, N b M 1 (0)=4, Tangled Cells: 1, the second smoothed mesh M2 is obtained with Jun algorithm, it has Q M 2 ∞ =0.109, N b M 2 (0)=0
, Tangled Cells: 0, the third mesh M3 is the intersected mesh resulting from [START_REF] Tipton | Grid optimization by equipotential relaxation[END_REF] with the following characteristic Q M 3 ∞ =0.117, N b M 3 (0)=0, Tangled Cells: 0.

be nearly optimal (it then becomes fixed for the next iterate). This number depends, of course, on the nodal quality function, we call it Q opt and we associate it to some nodal function (see definition 4) as 10% less than the optimal values:

Q opt =     
0.99, for sin f unction (81 degree) 0.9, for area, rel f unction.

A third consequence is that the nodal quality function gives a priority order that permits to mix algorithms. Basically, for nodes having acceptable quality we use some classical smoothing algorithms (or even do nothing) while for bad quality vertices, we use a more robust (but time consuming) one. For example, if Q ref is fixed in [START_REF] Scovazzi | Stabilized shock hydrodynamics: I. a lagrangian method[END_REF], and consider a mesh M to smooth for one iteration, and let

Q * ∈ [0, 1] = max(0, min(C 2 Q ref , C 3 + Q M ∞ )) (C 2 = 0.25, C 3 = 0.1) (31)
for each node n such that if (Q M n ≥ Q * ), we use an elliptic algorithm with Jacobi update (Tipton or Jun), but if (Q M n ≤ Q * ), we use an optimization algorithm with Gauss-Seidel update (RJM or Escobar type). As instance, we name the scheme mixed Jun-Escobar if Jun algorithm is used for the elliptic part and Escobar for the optimization part. The speed-up for this mixed strategy is between 5 and 10 (w.r.t. optimization applied to every nodes). Moreover, we can conserve symmetry, orthogonality and aspect ratio with the first two nodal quality functions (sin based) in definition 4. We name Smooth, the function that transforms mesh M in a smoothed mesh. It depends on M , proximity constraint [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF], and the nodal quality:

Smooth(M, ClosetoLag, {Q M n }). (32)
In fact, we propose the following strategy : for this part, we use equations ( 23)-( 32) and definitions 4 and 11.

Let M lag and Q ref be given,

M 0 = M lag , compute ClosetoLag 0
Choice of a nodal quality function Evaluation of Nodal quality and Global mesh quality of from one rezoning step to the very next one. To do this, we use some detector for symmetries and orthogonality. Moreover, at the end of the previous algorithm, we can use some time adaptative Q ref (at each rezoning step), depending on the quality numbers

M 0 if Q lag ∞ ≥ Q ref then M ale := M lag ,
Q M ∞ , Q M 1 and N b M (Q ref )
of the previous overall rezoning step. We decide to increase by some factor Q ref between two consecutive rezoning phases after two hydrodynamics step if: Increase Q ref : At least two of three assertions are verified:

Q M lag ∞ ≥ Q M ale ∞ , Q M lag 1 ≥ Q M ale 1 , N b M lag (Q ref ) ≤ N b M ale (Q ref ).
In this situation, this mean that Q ref is perhaps too small, so that the rezoning mesh is too close to Lagrangian grid and then we need to limit the Lagrangian nodal weight for the next step.

On contrary, we decrease by some factor

Q ref if: Decrease Q ref : At least two of three assertions are verified: Q M lag ∞ ≤ Q M ale ∞ , Q M lag 1 ≤ Q M ale 1 , N b M lag (Q ref ) ≥ N b M ale (Q ref ).
An interesting behavior of this algorithm is that it can be seen as a generalisation of ALE-n-cycle (smoothing process is applied after every n hydrodynamic Lagrangian steps), it gives an automatic and adaptive trigger of smoothing procedure. Examples of this strategy coupled with hydrodynamic is shown in section 5.

Remapping : Self intersection, discrete maximum principle and non-linear representation

In this section, we now have to compute the Lagrangian variables (ρ, U, E, ǫ, ...) on the new rezoned mesh in a conservative way. The aim of this part is to propose an approximate geometric scheme verifying a discrete maximum principle (local bound preserving) for the first order under a positivity condition and more precise than [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF]. For second order, we propose to enforce this discrete maximum principle using a local iterative high order reduction. Let us begin to recall the exact geometric problem and some approximation of [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF]. The notation for this chapter is:

M := M ale GOAL: REMAPPING Compute DATA : M C Q , Q N N Q Q Q Q A A N N C C C A A
(1) For a cell C, we note C' a neighbor cell (by edge if no set is precised), and |C| its volume, we will suppose here that all the volumes are strictly positive. ( 2) Let e be an edge (Figure 17), we suppose it is oriented from n1(e) to n2(e) (for example from the lower global node to the greater) and we call left(e) (resp. right(e)) the left cell (resp. right) to e. We note l C,e the local index of e in C. We then define the following sign function:

sign(C,e) =      +1 if e is trigonometric oriented in cell C, - 1 else.
(

) 33
We have of course sign(C,e)=-sign(C',e) for two adjacent cells C and C'. Let recall the following geometric relation (see [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF]) between the old and new mesh cells,

C N = C A ∪ ∪ C ′ ∈N (C A ) C N ∩ C ′A / ∪ C ′ ∈N (C A ) C A ∩ C ′N , ( 34
)
This relation is valid when C N is inside N (C A ), and is equivalent to the extended flux form:

C N QdV = C A QdV + C ′ ∈N (C A ) F Q CC ′ (35)
The quantity

F Q CC ′ def := [ C N ∩C ′A QdV -C A ∩C ′N QdV ]
is the exact flux obtained by the computation of exact intersection between the Lagrangian and rezone polygonal meshes), mainly by using a polygon-polygon intersection algorithm. Note that the local conservation comes from the fact that

F Q CC ′ = -F Q C ′ C .
and we will impose that approximate flux verifies this relation. The first order scheme can be written as:

|C N | QN C = |C A | QA C + C ′ ∈N (C A ) |C N ∩ C ′A | QA C ′ -|C A ∩ C ′N | QA C (36)
The positivity condition :

|C A | ≥ C ′ ∈N (C A ) |C A ∩ C ′N |,
gives a local discrete maximum principle (with respect to edge and corner neighboring cell). In order to avoid any intersection calculation, they approximate the multidimensional problem (35) by edge flux based on swept regions δF e which are algebraic volumes V (δF e ) swept by edges displacement. This procedure is by construction intersection free (Figure 18). In term of flux representation the first order Margolin-Shashkov [START_REF] Margolin | Second order sign preserving remapping on general grids[END_REF] scheme (here we put as superscript MS for their name, and 1 for the order) can be reinterpreted as (35) via fluxes:

F M S,1 CC ′ =      |V (δF CC ′ )| QA C ′ , if V (δF CC ′ ) > 0 -|V (δF CC ′ )| QA C , else. (37)
= max(0, V (δF

CC ′ )) QA C ′ + min(0, V (δF CC ′ )) QA C (38)
We rewrite the scheme in terms of individual edge flux:

F M S,1 e = V (δF e )Q e ,              Q e :=      QA right(e) if V (δF e ) ≥ 0, QA lef t(e)
else.

sign(lef t(e), e) = 1, sign(right(e), e) = -1.

(39)

The reader can then check that (37) rewrites into (39)

F M S,1 CC ′ = sign(C, e)V (δF e )Q e
We propose a modification of this flux which is based on segment-segment intersection (with at most one intersection point to compute per edge), this is an intermediate approach from time consuming polygon-polygon intersection and cheap no intersection computation algorithms. It basically consists in computing the intersection point of right situation in Figure 18. The reason is that the previous flux do not give any exchanges for null algebraic swept volume, see Figure 19: The first order self intersection flux writes:

δ -δ x x F e δ )=0 V(
F self,1 e = nblmt(e) k=1 V (Z e k )Q e k (40)
Where we use the notation:

(1) nblmt(e) is the number of cells in M A that e N goes through. It may takes only three values (0,1 or 2) corresponding to:

             0 ⇔ edge is either static or sliding, 1 ⇔ no self-intersection, (it recovers Margolin-Shashkov (37)), 2 ⇔ self-intersection. (41)
(2) Z e k is the k th sub-volume (Figure 20), and V(Z e k ) his algebraic volumes, note that if nblmt(e)=2, then we have always V(Z e 1 ).V(Z e 2 ) < 0. (3) Q e k is defined by:

Q e k =      QA C ′ if V (Z e k ) ≥ 0, QA C else.
In terms of discrete balance, our scheme writes:

F self,1 CC ′ =              F M S,1 CC ′ , if S * = {∅}, else |V (Z e 2 )| QA C ′ -|V (Z e 1 )| QA C , if V (Z e 2 ) > 0 |V (Z e 1 )| QA C ′ -|V (Z e 2 )| QA C , else. (44)
where |V (Z e k )| denotes the absolute value of the algebraic sub volume fluxed. In fact (44) can also be rewriten:

F self,1 CC ′ = nblmt(e) k=1 sign(C, e)V (Z e k )Q e k . (45
|C A | ≥ e∈C | nblmt(e) k=1,V (Z e k )<0 V (Z e k )| (46)
implies a discret local maximum principle (local bound preserving) with respect to edge neighborhood (and then on node neighborhood).

Short Proof.

( V (Z CC ′ k ) are equal. A sufficient condition is to compute all the volume (exact) with this formula (always possible in two dimension). We recall that the volume of any polygonal cell C given by trigonometric oriented nodes (M 1 , ..., M q ) is calculated by :

|C| = 1 2 q n=1 ( OM n ∧ OM n+1 ). e 3
Let ∆ n the displacement field for node n:

OM N n = OM A n + ∆ n , we compute ( OM N n ∧ OM N n+1 ) = ( OM A n + ∆ n ) ∧ ( OM A n+1 + ∆ n+1 ) = OM A n ∧ OM A n+1 + OM A n ∧ ∆ n+1 + ∆ n ∧ OM A n+1 + ∆ n ∧ ∆ n+1 .
After suming for all the nodes and using anti symmetry property of ∧ operator, that is to say q

n=1 OM A n ∧ ∆ n+1 + ∆ n ∧ OM A n+1
= 0, we obtain:

q n=1 ( OM N n ∧ OM N n+1 ) = q n=1 ( OM A n ∧ OM A n+1 ) + q n=1 ( ∆ n+1 ∧ ∆ n+1 )
which gives the result. This geometric constraint

|C N | = |C A | + e∈C nblmt(e) k=1 V (Z CC ′ k ) (47)
implies that a uniform constant state is preserved by the scheme, and it is a necessary condition to obtain the proof. (2) On second hand, we notice that the discrete balance of fluxes (40) writes:

             |C N | QN C = |C A | QA C + e∈C nblmt(e) k=1 V (Z CC ′ k ) Q * CC ′ with Q * CC ′ =      QA C if V (Z CC ′ k ) < 0, QA C ′ if V (Z CC ′ k ) > 0.
This is a convex combination iff (47) is verified with is given by the positivity condition (46).

Second order and Discrete Maximum principle remapping

Our second order remapping extension for self intersection flux (40) writes (we put 2 as superscript for second order schemes):

F self,2 e = nblmt(e) k=1 Z e k Q e k + ∇Q e k (x -x e k )dx (48)
where ∇Q e k (resp. x e k ) is a gradient evaluation (resp. the centroid) corresponding either to ∇Q C or ∇Q C ′ (resp. the cell C or C') depending on the upwinding rule for Q e k (42). The linear part in (48) is integrated by an exact quadrature formula.

Density and extensive quantities remapping

In the case of an arbitrary limiter of the gradient, it is well known that the second order approximation of density :

ρ A,2 = ρ A C + ∇ρ A C (x -x C ), x C centroid of C x C := 1 |C| C xdx. ( 49
)
does not warranty discrete maximum principle of density on new mesh M N . In [START_REF] Shashkov | The repair paradigm and application to conservation laws[END_REF], [START_REF] Loubère | The repair paradigm: new algorithms and applications to compressible flow[END_REF] the authors propose a conservative post-processing method to impose a discrete maximum principle (local bound preserving). We propose here an other alternative obtained by an always symmetric local iterative process, and based at most on local neighbors N(C).

Here, we use the result for first order (46) that is the corner stone of this approach. For second order scheme, if a cell does not verify the discrete maximum principle, we reduce the initial gradient contribution both inside and on the neighbor donor cells. We apply a more restrictive reduction factor on the limitation of gradients. To achieve it, we iteratively divide the limitation coefficient by a factor (two as instance) until the discrete maximum principle is obtained for cells violating it. For such a cell C and neighbor donor cells, we apply the following iterative process:

ρ A,2,(i) = ρ A C + 1 2 i (∇ρ A C ) 0 (x -x A C ). ( 50
)
Here, (∇ρ A s ) 0 is the initial gradient. We emphasize that this process does always converge (∃i; i < ∞) to a value such that discrete maximum principle is obtained for C. This function is continuous and tends to the first order representation that verifies this property.

Remark 14 [START_REF] Barlow | A compatible finite element multi-material ALE hydrodynamics algorithm[END_REF] This process can be extended to higher (more than second) order conservative reconstruction. (2) The numerical flux (48) with the reconstruction (50), can be seen as an other way to obtain results of Flux Corrected Transport in [START_REF] Schar | A synchronous and iterative fluxcorrection formalism for coupled transport equations[END_REF], [START_REF] Váchal | Sequential flux-corrected remapping for ALE methods[END_REF]. We want to keep high order evaluation everywhere if discrete maximum principle is satisfied and reduce order locally in an iteratively manner until it is obtained.

(3) This approach can be done on arbitrary extensive quantity, for intensive one see next subsection for specific quantities.

Specific quantities remapping

Here, we will give a way to obtain discrete maximum principle for specific energy (and velocity) for "second order" approximation. Let consider the following projections ρ A → ρ N and (ρe) A → (ρe) N so that the specific energy is defined by e N := (ρe) N ρ N . We have at least two ways to consider second order approximation:

(1) a linear representation for (ρe):

(ρe) A,2 = (ρe) A C + ∇(ρe) A C (x -x A C ) (51)
It is well known that there is no discrete maximum principle in e even if it is true for ρ and (ρe).

(2) In [START_REF] Vanderheyden | Compatible fluxes for van leer advection[END_REF], the authors propose to rewrite (51) as (ρe

) 2 = ρ C e C + e C ∇ρ C (x - x C ) + ρ C ∇e C (x -x C ) so that: e A,2,nl = e A C + ρ A C (x -x A C ) ρ A C + ∇ρ A C (x -x A C ) ∇e A C (52)
giving a non-linear representation of e.

In the case of staggered hydrodynamics, we do the same with the velocity (on the dual mesh). The iterative process applied to non-linear representation writes:

e A,2,nl,i = e A C + ρ A C (x -x A C ) ρ A C + ∇ρ A C (x -x A C ) 1 2 i ∇e A C ( 53
)
Remark [START_REF] Hirt | An arbitrary lagrangian-eulerian computing method for all flow speeds[END_REF] We suppose here that discrete maximum principle is obtained for density (for instance using ( 50)). This ensures the reconstruction (53) to be continuous, and the iterative process to converge also.

In the following, we perform a numerical study of the projection for specific energy and velocity. In our test, we use a gradient computed by green formula and we use the limiter in [START_REF] Dukowicz | Accurate conservative remapping (rezoning) for arbitrary lagrangian-eulerian computations[END_REF], and we don't need to apply (53) in our test cases. We want to show that self-intersection and non-linear representation as proposed in this paper are well suited for a remapping step of an ALE scheme.

For the test cases, we consider three options:

I Flux volume computation: (0) Without self-intersection procedure.

(1) With self-intersection procedure.

II Gradients computed by Green formula and limitation based on [START_REF] Dukowicz | Accurate conservative remapping (rezoning) for arbitrary lagrangian-eulerian computations[END_REF]. The initial data of the test cases on density and temperature are based on [START_REF] Vanderheyden | Compatible fluxes for van leer advection[END_REF], with the following velocity (Figure 21) below : For this initial data, we will consider three different "rezoning" test cases. Test 1 is diagonal advection see Figure 22, Test 2 is random cyclic rezoning Figure 23, and Test 3 is a cyclic rezoning Figure 24. For the latter, the advection speed is given by:

if 0 ≤ t ≤ 0.5      0.1 t sin(5πx)cos(5πy) 0.1 t cos(5πx)sin(5πy) if 0.5 ≤ t ≤ 1      0.1 (1 -t) sin(5πx)cos(5πy) 0.1 (1 -t) cos(5πx)sin(5πy) (54)
For all cases, we take cartesian grid with nx=ny=50. The remapping with the following options gives the best results:

(1) Self-intersection for the volumes fluxing.

(2) Linear limitation on ρ (primal cells) and non-linear on specific quantities.

(3) Linear limitation on ρ (dual cells) and non-linear on velocity.

It always ensures a discrete maximum principle on ρ, and specific quantities for first order scheme if and only if the positivity condition (46) is verified. For second order (in difficult cases) or for higher order methods, a local iterative process is sometimes needed on (50) for extensive quantities and (53) for intensive ones. We will focus on higher order schemes in a future paper.

5 Numerical results when coupling with hydrodynamic scheme

In this section, we now want to show how to use the tools developped in the previous section when coupled to a Lagrangian scheme. We recall that we propose:

(1) Rezoning in section 3: Nodal mesh quality tools, mixed Elliptic-Optimization algorithms.

(2) Remapping in section 4: Self-intersection flux volumes, non-linear limitation for specific quantities and velocity. Local iterative gradient reduction for discrete maximum principle enforcement.

In each of the following cases, we use the Wilkins scheme as hydrodynamic scheme. We also use pseudo-viscosity schulz-tvd [START_REF] Benson | A total variation diminishing shock viscosity[END_REF] with linear coefficient 0.5 and quadratic coefficient 0.7. The anti-hourglass [START_REF] Flanaghan | A uniform strains hexaedron and quadrilateral with orthogonal hourglass control[END_REF] coefficients is choosen in such a way:

(1) For the Sedov test case (Figure 25), we take standard hourglass coefficient as 0.25. (2) For the planar shock case (Figure 30) and three materials case of [START_REF] Després | The vofire method for multi-component flows on unstructured meshes[END_REF] see Figure 32, we do not take any hourglass method.

Sedov test case

The Sedov test case is defined by the initial state (ρ = 1, P = 0, U = 0), γ = 5/3, on an initial cartesian grid 50 × 50. A dirac energy ǫ = 5000 is put in the cell in contact with the origin. Remark 16 With the nodal mesh quality control, we can note that:

(1) the Lagrangian mesh is only smoothed wherever it is effectively distorded.

(2) With the adaptation of Q ref , the total effective number of smoothing is very small (because

Q M lag ∞ > Q ref )
, moreover the number of steps between two such effective smoothing is not constant in time, because it can not be known in advance. Here, we decide to relax (to consider a degree of freedom) a node until the end from the moment where an angle around it is less than π 4 . Note that our new algorithms (Bottom) and strategy is superior in both shock arrival times and numerical diffusion with respect to a standard smoothing algorithm (Top). Note here that:

• In this test case we use isothermal/isobar mixing modeling.

• The Lagrangian computation produces invalid grids, causing the time step to be very small and finally to stop. • The ALE where rezoning is done for all nodes at each time step (Top) produce a quasi Eulerian grid, causing more numerical diffusion see oblique shock (Figure 32 Middle and Bottom) . • The nodal mesh quality ALE (Bottom in Figure 32 or Bottom Figure 33) produces smoothed grid only in regions with poorly positioned vertices. This strategy permits to gain naturally both accuracy and small mass fluxing.

Conclusion and future works

In this paper we have proposed a strategy for Arbitrary Lagrangian Eulerian computation. The rezoning step is done under nodal quality control indepen-dently of the smoothing algorithm. This definition is based on some simple rules and permits to compare iso-connectivity meshes. This definition gets several consequences. Mainly, in the case where we want the smoothed mesh to be close to the Lagrangian, it naturally gives a way to control it through a non linear nodal mesh relaxation. Second, it gives two global mesh quality definition that permit to track non convex cell. And thirdly, it gives a priority order for node that allows to mix different smoothing algorithms (depending on the values of the node quality) and also updating method (Jacobi or Gauss Seidel).

A future extension of the present work will be to design such definition when physical weights (scalar or tensorial) need to be taken into account. Moreover, we have proposed to extend the Escobar algorithm [START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF] to heterogenous conformal mesh: it gives naturally different Jacobian matrix for each simplex.

A natural extension would be to adapt it to anisotropic nodal mesh quality when weights or tensor are defined in the smoothing process.

For the remapping step, we propose to take into account a segment-segment intersection (simpler than a full polygon-polygon intersection) and compare result for density, velocity and internal energy. The second order extension for specific quantities use the non linear limitation of [START_REF] Vanderheyden | Compatible fluxes for van leer advection[END_REF]. We show results on the coupling between this new approximate flux volume scheme and this non linear limitation. Moreover, we proposed an iterative process to enforce the discrete maximum principle both on extensive and intensive quantities using this new approximation. Future work based on this strategy will be investigated: it involves higher order remapping which are forced to verify a discrete maximum principle by iterative processes when density and internal energy are computed. Finally, we have shown that the overall tools developped here are well suited when coupled with Lagrangian hydrodynamic schemes to control mesh displacement for rezoning phase and discrete maximum principle (local bound preserving) for remapping phase.

  d dt C(t) ρUdx + ∂C(t) P N dl = 0, d dt C(t) ρEdx + ∂C(t) (P U).N dl = 0. dl = 0 for volume conservation.

  is the set of nodes of cell C, #N (C) is the number of vertices in cell C, n C l : l is the local index of n in the cell C.

Figure 1 .

 1 Figure 1. Connectivity cell/node (counter-clock wise) local index n C l and node/cell local index C n l

Figure 2 .

 2 Figure 2. Neighborhood of vertex n in a cell: short/medium/long range and underlying cutting into simplices (see triangles T above) for three consecutive nodes involving n.

1 2 (Figure 3 .

 23 Figure 3. Extension of functionnal by denominator modification construction of such rational function in IR 2 , let T a simplex in the physical spacex k = (x 1 k , x 2 k ) ∈ IR 2 , k=0,1,2 and T R the reference simplex u 0 = (0, 0), u 1 = (1, 0),u 2 = (0, 1). The affine map between T and T R is given by x = x 0 + Au, where A is the Jacobian matrix of the map referenced to node x 0 .

Figure 4 .

 4 Figure 4. Different transformations between Real, Reference and Ideal triangle.

Figure 5 .

 5 Figure 5. Neighborhood of node, and adaptation of local Escobar functional to local connectivity for short and medium neighborhood after cutting cells by simplices.Here, C i are the i t h local cell linked to node n, T S i the SHORT range sub-simplex inside C i , and T i,1 , T i,2 the MEDIUM range sub-simplices. The angle α i (resp. α i,1 ,α i,2 ) is the angle inside T S i (resp. T i,1 , T i,2 ) from n.

Figure 6 .

 6 Figure 6. Medium range optimal coefficient

Figure 11 .

 11 Figure 11. Evolution of the nodal quality for different position in one cell

Figure 12 .

 12 Figure 12. Patch for the different ranges : short and medium range

120 OFigure 13 .

 12013 Figure[START_REF] Garimella | Untangling of 2d meshes in ale simulations[END_REF]. Quality of central node. Left patch: Q = 1 for all these functions, Middle patch: Q sin = 1, Q sin,rel = 1, Q aire,rel = .25, Right patch: Q sin = .5, Q sin,rel = 1, Q aire,rel = 1. for all these functions.

Figure 16 .

 16 Figure 16. Mesh before (left) and after rezoning (right) : Remapping, quantity Q C (resp. Q n ) stands for centered quantity (resp. for nodal), and the N (resp. A) super script for Q stands for new quantity (resp. old).

Figure 17 .

 17 Figure 17. Edge connectivity (3) N (C) the set of all C' that share at least one vertex with C.

Figure 18 .

 18 Figure 18. First order flux based on area swept by edge displacement δF e of Margolin-Shashkov flux: donor cell interpretation (only one swept region at most)

Figure 19 .Figure 20 .

 1920 Figure 19. No contribution for zero algebraic first order flux volume of scheme (37), far from exact solution if large discontinuity.

  t(e), e) = +1, sign(right(e), e) = -1.

) Proposition 13

 13 If the flux volumes (40) and cell volumes |C A | and |C N | are computed with the same formula, then the following positivity condition:

)

  First we have to prove that each volume |C N | of new cells and the volume balance |C N | := |C A | + e∈C nblmt(e) k=1

  (a) The limitation for specific quantities s C (projection of (ρs) C ) (0) Non linear limitation (first ρ C then s C ) (1) Standard linear limitation (on (ρs) C ) (b) and velocity (-1) First Order (0) Non linear limitation (same as specific quantity) (1) Standard linear limitation (same as specific quantity)

Figure 21 .

 21 Figure 21. Initial data for test cases Planar geometry (density ρ and temperature T on Left, velocity U=(u,v) on Right)

Figure 23 . 1 Figure 24 .

 23124 Figure 23. Random cyclic rezoning and resulting temperature: No global discrete maximum principle on density and temperature if no self-intersection

Figure 25 .Figure 26 .

 2526 Figure 25. Lagrangian mesh and corresponding density at time 1e-5s

Figure 27 .Figure 28 .

 2728 Figure 27. Lagrangian mesh and corresponding density at time 2e-5s

Figure 31 .

 31 Figure 31. ALE computation without nodal mesh quality control using Jun (Top), and with nodal mesh quality controld with Mixed Jun/Escobar (Bottom) based on remark/equation (31) at times 20ns see Figure30for pure Lagrangian computation. Here, we decide to relax (to consider a degree of freedom) a node until the end from the moment where an angle around it is less than π 4 . Note that our new algorithms (Bottom) and strategy is superior in both shock arrival times and numerical diffusion with respect to a standard smoothing algorithm (Top).

Figure 33 .

 33 Figure 33. Resulting mesh/concentration of Figure 32. Top : Without any nodal quality treatment, and Bottom with nodal quality control.

  Comparaison with different variation of Escobar functionnal (20) on a simple initial distorded mesh

			Test 2: Untangling of an unstructured simplicial mesh Test 3: Untangling of an unstructured quadrilateral mesh
			Mesh M1			Mesh M1
	2.0×10 0	Iter no 1		Iter no 2 1.67 1.67	iter no 3
				Mesh M1
	2.0×10 0
	1.0×10 0
					1 1
	-4.4×10 -9				0.5078 0.5078
	-2.1×10 -9		1.0×10 0	2.0×10 0	0.4578	1	1.653
		iter no 4 Figure 9. Untangling of quadrilateral mesh. Tangled mesh (red) and untangled mesh iter no 5 iter no 6
	1.0×10 0	(blue) after 2 iterations, Left Fullview, Right Zoom
				Mesh M1
	2.0×10 0
	-7.1×10 -9
	-1.8×10 -9 1.0×10 0	1.0×10 0	2.0×10 0
			Our connectivity adaptation algorithm is always better than the original Escobar Figure 8. Untangling of static simplicial mesh : (red) Tangled mesh, (blue) Final
			(always in the middle) smoothed mesh
	0.0×10 0
			0.0×10 0	1.0×10 0	2.0×10 0

. The second (resp. third) test Figure

[START_REF] Dukowicz | Accurate conservative remapping (rezoning) for arbitrary lagrangian-eulerian computations[END_REF]

(resp. Figure

[START_REF] Escobar | Simultaneous untangling and smoothing of tetrahedral meshes[END_REF]

) consists in untangling a simplicial mesh (resp. unstructured quadrilateral mesh). Test 1: Figure 7. Comparison of evolution w.r.t global update (Gauss Seidel) for three different variations of Escobar algorithm: Original Escobar (red, all α = π/3), Escobar ortho. (blue, all α = 0), and new Escobar (green, α connectivity adaptation (20)).

Figure 10. Smoothed mesh after 20 global iterations (blue) obtained from the tangled mesh (red)

  stop else (2) Compatible with any kind of nodal statute (depending on criteria such as angles, area, neighborhood). (3) Adaptive choice (in time) of the nodal quality function and adaptation of Q ref . The function that analyses all the meshes can jump between Q sin n and Q area,rel n

  Figure 22. No global discrete maximum principle on temperature for linear limiter.

							Test 1 : Influence on specific quantity
				Test Case : Advection (diagonal direction (0.005,0.005))
		Initial density min/max :0.001/1, initial temperature min/max : 0/1
	nx=ny				Flux				cell limiter	density min/max Temperature min/max
	50		self intersection or not Non-Linear 0.001/0.953157874	0/0.999775443
	50		self intersection or not	Linear	0.001/0.953157874	0/1.023732158
						Temperature										Temperature
												0.9								1
		0.9 0.9998										0.8		1						0.9
		0.8										0.7								0.8
		0.7																		0.7
												0.6
	Z-Axis	0.4 0.5 0.6										0.4 0.5	Z-Axis							0.5 0.6
		0.3																		0.4
												0.3
		0.2																		0.3
		0.1										0.2								0.2
		0	0.3									0.1		0	0.3					0.1
		Z Y -A x is 0.4	0.5	0.6	0.7	0.7	0.6	0.5	X-Ax is 0.4	0.3	0		Z Y -A x is 0.4	0.5	0.6	0.7	0.7	0.6	X-Ax is 0.4 0.5	0.3	0
	X		Y										X		Y
							Test 2 : Random Cyclic Rezoning
	nx=ny				Flux			cell limiter	density min/max		Temperature min/max
	50		Self intersection	Non linear	0.001/0.999994937		0/0.999999724
	50		Self intersection			linear	0.001/0.999994937		0/1.000276494
	50	No Self intersection Non linear 0.000999984/0.999997086 -0.000008013/0.999999854
	50	No Self intersection		linear	0.000999984/0.999997086 -0.000015653/1.000144011

|S||S -1 | , the

n end if end for Evaluation of Nodal quality and Global mesh quality of

We use non linear relaxation [START_REF] Shashkov | The repair paradigm and application to conservation laws[END_REF] with M 1 n , M 2 n beeing two of the following meshes : Lagrangian M lag , Best M BEST , Last Iteration M M axIterSmooth , Previous M prev , or any intermediate (Lagrange/Rezoned or Previous/Rezoned). Moreover, for each of the non linear relaxation, we replace Q opt by a quantity depending on

In short terms, this approach can be viewed as an over layer of some previous versions of an ALE smoothing algorithm. At the end, we can check

for both Lagrangian and final rezoned meshes (named M ale ). This approach is:

(1) Independent of the smoothing algorithm.

The case of Planar shock with non regular grid

The initial gas is at rest, and the density is equal to 1, γ = 1.4, and we impose to the left boundary condition a constant velocity v=(vz,0) so that an initial shock is coming from left and traveling to the right (vz=5.10 7 ). The mesh is depicted in Figure 29. It is composed with unstructured quadrilaterals. We run this test in pure Lagrangian Figure 30 and in two ALE regimes Figure 31.