P Hoch

Nodal extension of Approximate Riemann Solvers and nonlinear high order reconstruction for finite volume method on unstructured polygonal and conical meshes: the homogeneous case

In this article, we consider the nodal flux extension of classical Eulerian edge flux schemes for linear and nonlinear systems of two dimensional hyperbolic conservation laws. In a first part, we present structural properties implying local conservativity and consistency. Such criteria apply to flux based scheme (VFFC, Roe,..) as well as viscosity based scheme (Rusanov, HLL, VF-Roe,..) (invariant domain preservation and entropy stability are still highly scheme dependent). We propose an algebraic framework that straightforwardly applies on arbitrary unstructured polygonal or curved conical meshes. This naturally defines composite edges/nodal fluxes taking into account either straight or curved meshes. We extend it to finite volume methods but it could also be used for a Discontinuous Galerkin discretization or Residual Distribution scheme. In a second part, we deal with the stability (of discrete unknowns) and admissibility (in flux evaluation) of Euler system of gas dynamics. We focus on a nonlinear reconstruction for massic variables and propose a way to circumvent the direct velocity limitation which is not physically relevent. Finally, we show how to deduce from relation of total energy E = + 1 2 |U| 2 a limited velocity reconstruction induced by the direct limitation of internal massic energy .

Introduction

In this article, we consider the following conservative system to be solved on domain Ω ∈ R d (d = 2 here):

∂ t U + div(F (U )) = 0, U (0, x) = U 0 (x).
(1.1)

where U , the conservative unknown is such that U (t, x) ∈ R n and F (.) = (F 1 (.), .., F d (.)) the (linear or nonlinear) flux is such F (.) ∈ R n×d , and div(F (U )) = d i=1 ∂ xi F i (U ) in Cartesian coordinates. Here we have the commonly used hypothesis and notations (see [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]):

1. Hyperbolicity: ∀ξ ∈ R d , with |ξ| = 1, the jacobian matrix J(U , ξ) := ∂F (U ).ξ ∂U = i=1,d

∂F i (U ) ∂U ξ i is diagonalizable.(1.2)
2. We suppose that the n real eigenvalues λ i=1,..,n of (1.2) are ordered:

λ 1 (U , ξ) ≤ .. ≤ λ i (U , ξ) ≤ ... ≤ λ n (U , ξ). (1.3)
3. We note by R(U , ξ) (resp. L(U , ξ)) the right (resp. left) eigenvectors of Jacobian matrix (1.2). We suppose that the following equality holds for all U and ξ : R(U , ξ)L(U , ξ) = I (the identity matrix), (1.4) and we call Diag(λ i ) the diagonal matrix containing eigenvalues λ i (U , ξ):

J(U , ξ) := R(U , ξ)Diag(λ i )L(U , ξ).

(1.5)

4. We suppose that U 0 (.) does live in an admissible (convex) set A adm for the associated flux F (1.1) (ie physically admissible states):

A adm ⊂ R n . (1.6)
We recall that in the multi-dimensional case, stability results are weak and a common approach is to deal with admissible states of (1.1). More precisely, we recall what is an invariant domain (see e.g. [START_REF] Seguin | Méthode de volumes finis pour les fluides compressibles[END_REF][25]):

Definition 1.1. A (convex) invariant set C ⊆ A adm ∈ R n is said to be an invariant domain for (1.1) if for any U 0 (x) ∈ C then solution U (t, x) is also in C for all t > 0.

Finite volume method comes from integrating (1.1) on each cell Ω j : Ωj ∂ t U + div(F (U ))dv = 0, and most edge based fluxes are based on applying the Green-Riemann formula on it:

∂ t U j (t) = - 1 |Ω j | ∂Ωj F (U ).Nds, (1.7)
where the discrete unknown are U j (t) = 1 |Ωj | Ωj U (t, x)dv. In practice, evolving (1.7) is done using a numerical flux defined at edge of ∂Ω j . This function must take into account information from both parts of each edge: two values U j , U i and a (local) normal direction N ji coming from a point located on the edge e ji (Figure 1). Now consider a first order scheme for (1.7) on straight edges, it writes (see e.g. [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]):

∂ t U j (t) = - 1 |Ω j | eji |e ji |Φ(U j (t), U i (t), N ji ). (1.8)
N ji e ji Ω j Ω i

Mr+1 Mr

Mr-1

N ji e ji Ω j Ω i Ω k Ω n
For edge based numerical flux Φ: green cells are involved in exchanges with Ω j red cells are ignored in exchanges with Ω j ⇒ ignored cells may be more numerous than thoose involved The function defining a numerical flux Φ must fullfills some necessary (and not sufficient) requirements: local conservation and consistency, those can be any of the well known schemes: Godunov, Roe, Rusanov, HLL... Our first motivation is to extend the classical eulerian edge finite volume method to eulerian nodal finite volume method. For example, see Figure 1, we want our numerical flux to make it possible to exchange informations contained in Ω j with adjacent edge (green) cells (Ω i , Ω k , ..) but also with corner (red) cells (Ω n , ..).

As mentionned below, pure nodal finite volume method enjoy good properties. In Lagrangian update hydrodynamic [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], such a nodal discretization enable to obtain a natural local update of the mesh which is coherent with Geometrical Conservation Laws, in the context of pure edge finite volume method this is actually not the case [START_REF] Benson | Computational methods in lagrangian and eulerian hydrocodes[END_REF]. Nodal flux also enable to obtain Asymptotic Preserving properties on unstructured meshes [START_REF] Buet | Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes[END_REF][START_REF] Del Pino | An asymptotic preserving multidimensional ale method for a system of two compressible flows coupled with friction[END_REF]. Nevertheless sometimes they exhibit some parasitic modes that may be cured with ad-hoc processes. We can found a first example in [START_REF] Franck | Design and numerical analysis of asymptotic preserving schemes on unstructured meshes. Application to the linear transport and Friedrichs systems[END_REF] where a process permits to clean the cross stencil phenomenom of [START_REF] Buet | Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes[END_REF]. A second example can be found in the seminal work [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF] where the pure nodal GLACE scheme exhibit unphysical hourglass modes. Note also that up to our knowledge, no generic nodal Approximate Riemann Solver do exists for arbitrary hyperbolic systems writted in all frames: Eulerian, updated Lagrangian.

Figure 2: Numerical solution of hyperbolic P1 model on cartesian mesh with Dirac like Cauchy data see [START_REF] Buet | Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes[END_REF][START_REF] Blanc | Asymptotic preserving schemes on conical unstructured 2d meshes[END_REF]. Left: with a pure nodal polygonal scheme. Right: with a (composite) conical degenerate scheme. The pure nodal scheme exhibits some cross stencil unphysical phenomenom (here cured by the composite scheme, see Figure 5 and section below), both are first order in time and space.

Figure 3: Numerical solution at t = 1 for Sedov test case (Euler gaz dynamics) solved with Lagrange update formulation. Left: with a polygonal (pure nodal) GLACE scheme [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF]. Right: with a (composite) GLACE conical degenerate extension [START_REF] Boutin | Extension of ALE methodology to unstructured conical meshes[END_REF][START_REF] Bernard-Champmartin | Extension of centered hydrodynamical schemes to unstructured deforming conical meshes : the case of circles[END_REF]. The pure nodal scheme exhibits numerical instability due to lack of viscosity in diagonal direction (cured here by the composite scheme, see Figure 5 and section below), both are first order in time and space.

In parallel, pure edge finite volumes schemes well behave with respect to problems whose mesh is aligned with the flow, and may benefit most of the 1-dimensionnal numerical flux design theory. Moreover these are mainly used in industrial codes with a lot of publications and benchmarks. Nevertheless, for first order explicit finite volume schemes the propagation is done by cell/edge adjacency, this may cause a severe cfl constrain on some given meshes, specifically for nodes having a lot of support cells (see node at (0.55, 0.55) in Figure 4). A drawback in using pure edge scheme is the lack of Asymptotic Preserving property for arbitrary unstructured meshes see [START_REF] Buet | Design of asymptotic preserving finite volume schemes for the hyperbolic heat equation on unstructured meshes[END_REF] (Delaunay hypothesis on mesh in order to obtain consistency).

As motivated by pros and cons behavior of pure nodal and pure edge schemes, we focuse on composite (or nodal/edge) scheme because it is both defined on edge (quasi one-dimensional) and at node (multi-dimensional) whatever the order of the scheme. With such a mixed localisation in flux computation, we hope (and expect) the benefits of each of the two classes. Thereafter, we will give a precise definition of such new general composite schemes.

From previous works, we also want to deal with conical meshes, because they will give special cases of composite schemes defined both on straight and curved meshes, and we give a generalisation of some known schemes defined for two dimensionnal linear advection system (e.g. [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF]).

More precisely, we summarize our requirements:

1. For a generic cell Ω j , we want to take into account all surrounding cells, that is not only those sharing an edge but also those sharing a node (see Figure 1), thus accessing to the corner (red) cells. Then, we would like to be able to give a qualitative measure of edge's importance versus node neighboring cells inside the discrete balance laws.

2. We want to deal with arbitrary polygonal and conical cells. Such meshes are formed by curved edges obtained by adding to each straight edge a control point and a non-negative weight ω that control the curvature at mid edge point. Extensions of Lagrangian-update hydrodynamic schemes such as GLACE [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF][START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF] and EUCCLHYD [START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF] have been proposed in [START_REF] Boutin | Extension of ALE methodology to unstructured conical meshes[END_REF][START_REF] Bernard-Champmartin | Extension of centered hydrodynamical schemes to unstructured deforming conical meshes : the case of circles[END_REF].

S 0 r+1/2 S 0 r+1/2 = 1 2 (M r + M r+1 ) Ω j Ω i 0 ← ω M r+1 S ω r+1/2 M r Nodal Flux Classical Edge Flux
Composite Nodal/Edge Flux 1) while nodal fluxes involve both green and red ones. Straight polygonal meshes must be a special case of conical meshes for ω = 0 (in this case the mid-edge flux still contribute).

In [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF], it is applied on two-dimensional vectorial advection and we compare classical edge flux with composite nodal/edge flux up to third order on curved conical meshes. We also use it both on the hyperbolic heat equation (P1 system) [START_REF] Blanc | Asymptotic preserving schemes on conical unstructured 2d meshes[END_REF] to construct new asymptotic preserving schemes, and on the limit diffusion equation obtaining new diffusion schemes adapted to both straight or curved meshes.

The common novelty is that in the limit ω → 0 of a conical degenerated mesh (straight edge), the conical schemes give truly new schemes different from pure polygonal ones. Indeed the degree of freedom located at mid-edge (S ω r+1/2 ) does not disappear (see Figure 5). Moreover, it should be noticed that they may show better stability when compared to purely nodal polygonal schemes (cross stencil phenomenom).

3. We want to obtain a continuous (with respect to a scalar weight attached to each edge) extension of classical edge fluxes (Rusanov, HLL, Roe, VFFC, ..) from edges to nodes in such a way that a genuinely curved version does not need much extra effort (for the extension of fully multi-dimensional HLLC solver see [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF]).

4. We want to obtain third order schemes on such meshes under the constraint that both density and internal massic energy remain non-negative during the fluxes evaluation, using same technics of real degree valued reconstruction [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF], suitable for both a priori (see [START_REF] Maire | Positivity-preserving cell-centered Lagrangian schemes for multimaterial compressible flows: From first-order to high-orders. Part II: The two-dimensional case[END_REF]) and a posteriori limitation framework ( [START_REF] Clain | A high-order finite volume method for systems of conservation laws multi-dimensional optimal order detection (mood)[END_REF][START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF]).

The article is organized as follows: in Section 2, we introduce the definition of polygonal and conical meshes, and also the notation of geometrical objects needed in the construction of composite schemes. Especialy, we propose a construction of normal vectors which makes it possible to deal with both edge and nodal finite volume fluxes. In Section 3, we present the abstract form of both flux based and viscosity based numerical fluxes. We extend some known Approximate Riemann Solvers from edges to nodes on such meshes and in such a way that they are consistent and localy conservative, giving new extended finite volumes schemes called CARS (Composite Approximate Riemann Solvers). In Section 4, we introduce some nonlinear reconstructions on massic variables to deal with high order accuracy (second and third order). For example, on the gas dynamic system, using the notion of indirect limitation, we show how to circumvent a direct limitation of the velocity (or momentum) variable. Section 5 is devoted to numerical test-cases using results of sections 3 and 4. We show some comparisons of different numerical composite solvers, both on regular and discontinuous solutions. Variations are made upon the order of the reconstruction, but also on different limitations strategies.

2 Conical meshes and composite geometrical informations

Definition and basic properties

In all the paper, we consider a domain Ω ⊂ R 2 which is Lipschitz and bounded. We assume that a mesh M is given on Ω. By mesh, we mean the following:

1. (Ω j ) 1≤j≤N is a finite family of nonempty connected open disjoint subsets of Ω (the cells) such that

Ω = N j=1 Ω j ;
2. (M r ) 1≤r≤P is a family of P distinct points (the nodes) of Ω such that, for all r, M r ∈ ∂Ω j , for some j;

3. for all j, the boundary ∂Ω j of Ω j is the union of a finite number of edges, each of which is a smooth one-dimensional curve joining two nodes (say, M r and M r+1 ) of Ω j . Moreover, two different edges share at most one node. Finally, any edge is contained either in ∂Ω or in Ω j ∩ Ω k for two distinct cells Ω j and Ω k .

4. we call E the set of all edges of the mesh M.

The edges are not necessarily segments, so we describe below how we parametrize them. Let M 0 and M 2 be two points linked by a regular arc γ arc M0,M2 . It is assumed that this arc can be represented by a closed part zero-level curve of a function Γ: Γ : R 2 -→ R (x, y) -→ Γ(x, y) such that γ arc M0,M2 admits an implicit representation:

γ arc M0,M2 ⊂ {(x, y) ∈ R 2 ; Γ(x, y) = 0}. Definition 2.1. If the implicit function Γ is a bi-variate polynomial of (maximal) degree p, we say that the arc γ arc M0,M2 is of order p + 1. In particular:

Linear form P 1 (x, y) : a x + b y + c = 0, straight lines. (2.1) Quadratic form P 2 (x, y) : a x 2 + b y 2 + 2 c x y + d x + e y + f = 0, conics. (2.2)
We recall that solutions of the quadratic equation (2.2) generate the family of conics of the plane.

Instead of considering only straight edges (2.1), we are interested in conical curves (2.2), in the context of finite volume methods for which fluxes are computed on edges. Proposition 2.2 (Quadratic Rational Bezier curve). Conical arcs may be parameterized by rational Bezier curve of degree 2, charaterized by 2 endpoints M 0 and M 2 , a control point M 1 , and a scalar non negative weight ω (cf Figure 6):

M ω (q) = x(q) y(q) = M 0 (1 -q) 2 + 2ωq(1 -q)M 1 + q 2 M 2 (1 -q) 2 + 2ωq(1 -q) + q 2 , q ∈ [0, 1]. (2.3)
Hence, we have: Each conical edge depends on its two endpoints M 0 , M 2 and 3 (in dimension 2) additional parameters: the control point M 1 and a scalar weight ω. The link with conics comes from the following characterization: Proposition 2.3. The edge is said to be of type:

γ arc M0,M2 = {M ω (q), q ∈ [0, 1]}. M 0 (M 1 , ω) M 2 M ω (q) M 0 (M 1 , ω) M 2 ω = 0 ω = +∞ ω = 0.5 ω = 1 ω = 3
-planar degenerate if ω = 0, -elliptic if ω ∈ ]0, 1[, -parabolic if ω = 1, -hyperbolic if ω ∈ ]1, +∞[.
In addition, we introduce the shoulder point S located at the midpoint of arc parametrization, S := M ω (q = 0.5). It may be noted that the unit normal at this point is orthogonal to

[M 0 , M 2 ] M 0 (M 1 , ω) M 2 S Q 0 Q 2 Figure 7: Shoulder point S.
The shoulder point S satisfies the following property:

S = 1 2 (Q 0 + Q 2 ), with Q 0 = 1 1 + ω (ωM 1 + M 0 ), Q 2 = 1 1 + ω (ωM 1 + M 2 ).

Relationship on geometric objects

We recall that the area of a cell is given by:

|Ω j | = Ωj 1dv = 1 2 ∂Ωj OM(s) • N(s)ds (2.4) = e 1 2 1 0
OM ω,e (q) • N e (q)dq (2.5)

where N(s) is the external unit normal of cell Ω j whose boundary is parametrized by s and N e (q) is the non unit one (N e (q) = Ne (q) ds dq ). The point O is the origin, but it can be any fixed point in the plane. In practice, it is usually set as the centroïd of the cell Ω j .

|Ω j | = e (A(O, M e 0 , M e 2 ) + A(M ω,e (q), M e 0 , M e 1 , M e 2 )), (2.6)
where A(O, M e 0 , M e 2 ) is the area of a simplex defined by origin O, M e 0 and M e 2 vertices of edge e, and A(M ω,e (q), M e 0 , M e 1 , M e 2 )) corresponds to the area between straight edge M e 0 M e 2 and the curved arc defined by (2.3). In addition,

A(M ω,e (q), M e 0 , M e 1 , M e 2 ) = 1 2 1 0
x(q) dy dq (q) -dx dq (q)y(q) dq.

(2.7)

Moreover:

A(M ω,e (q), M e 0 , M e 1 , M e 2 ) = f (ω)A(M e 0 , M e 1 , M e 2 ), (2.8)
with f (ω) a piecewise analytical function of ω of class at least C 1 (R + ) (see [START_REF] Bernard-Champmartin | Extension of centered hydrodynamical schemes to unstructured deforming conical meshes : the case of circles[END_REF]). The exact formula is [START_REF] Guojin | Computing integral values involving nurbs curves[END_REF], [START_REF] Boutin | Extension of ALE methodology to unstructured conical meshes[END_REF]:

           (1) if ω = 0, then f (ω) = 0, ( 2
) if 0 < ω < 1, then f (ω) = ω 1-ω 2 2 √ 1-ω 2 atan( 1-ω 1+ω ) -ω , (3) if ω = 1, then f (ω) = 2 3 , (4) if ω > 1, then f (ω) = ω ω 2 -1 ω -2 √ ω 2 -1 atanh( ω-1 ω+1 ) .
(2.9)

We recall that the area of Ω j , whose boundary is made of an arbitray number of conical edges, may be expressed from vectors, cf [START_REF] Boutin | Extension of ALE methodology to unstructured conical meshes[END_REF] and Figure 8:

|Ω j | = 1 2 dof C dof,ω j .OM dof . (2.10)
Here, dof stands for degrees of freedom, this terminology is used as an abuse of notations and may be assimilate to all control points defining the cell. Hence, the sum runs over all degrees of freedom in the cell (that is to say, the nodes and the control points of the edges). See Figure 8, noting N r-1,r := -N r,r-1 and N r,r+1 := -N r+1,r ,

M r+1 M r+1/2 M r M r-1/2 M r-1 N r + 1 , r N r , r - 1 N r + 1 / 2 ,r + 1 N r, r + 1 / 2 Nr - 1 / 2 ,r N r-1 ,r-1 /2 C r + 1 / 2 ,ω j C r,ω j C r -1 / 2 , ω j Ω j
C r,ω j = 1 2 (1 -f (ω r-1/2 ))N r-1,r + (1 -f (ω r+1/2 ))N r,r+1 + f (ω r-1/2 )N r-1/2,r + f (ω r+1/2 )N r,r+1/2 , (2.11) C r+1/2,ω j = f (ω r+1/2 ) 2 (N r,r+1/2 + N r+1/2,r+1
).

(2.12)

Due to property of (2.9), we recall that the dependency property of such C dof,ω j (2.11) (2.12) based on Computer Aided Design nodes M r+1/2 :

C r+1/2,ω j → w→0 0, C r,ω j → w→0 C r j := 1 2 (N r-1,r + N r,r+1
).

(2.13) so that the conical schemes based on (2.13) tend to their polygonal version [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF][START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF].

Despite the concise writing of (2.10), and property (2.13), a disadvantage is that the CAD control point M r+1/2 does not lie on the arc γ arc M0,M2 (see (2.3)). This is why we will use a description of the arc with its shoulder point S, cf Figure 9 (at the end), as in [START_REF] Bernard-Champmartin | Extension of centered hydrodynamical schemes to unstructured deforming conical meshes : the case of circles[END_REF].

To summarize, (2.8) rewrites:

A(M ω,e (q), M e 0 , M e 1 , M e 2 ) = h(ω)A(M e 0 , S e , M e 2 ), (2.14)
This change in the spatial location of the degree of freedom involves now the function

h(ω) := f (ω) ω + 1 ω (2.15)
which is strictly decreasing on R + with values in ]1, π 2 ]. This allows us to derive the following formulas, involving only the points that belong to the edges of cell Ω j , hence from now dof should be understood as quadrature points.

1. GLACE type:

|Ω j | = 1 2 dof Cdof,ω j .OM ω dof (q), (2.16)
that is,

|Ω j | = 1 2   r Cr,ω j .OM r + r+1/2 Cr+1/2,ω j .OS ω r+1/2   := 1 2 newdof Cnewdof,ω j .OM newdof (2.17)
with

Cr,ω j = 1 2 (1 -h(ω r-1/2 )) Ñr-1,r + (1 -h(ω r+1/2 )) Ñr,r+1 + h(ω r-1/2 ) Ñr-1/2,r + h(ω r+1/2 ) Ñr,r+1/2 ,(2.18) Cr+1/2,ω j = h(ω r+1/2 ) 2 ( Ñr,r+1/2 + Ñr+1/2,r+1 ). (2.19)
Unknowns are located at the same points as Cdof,ω j : points M r , S ω r+1/2 inside cell, moreover they satisfy [10][6]:

For cell Ω j : dof Cdof,ω j = 0 (2.20)
For all degree of freedom (dof) inside the domain (Figure 12)

: j Cdof,ω j = 0 (2.21)
2. EUCCLHYD type:

M r+1 S ω r+1/2 M r S ω r-1/2 M r-1 Ñr + 1 , r Ñr , r -1 Ñr + 1 / 2 , r + 1
Ñr,

r+ 1/ 2 Ñr-1/ 2, r Ñr -1 , r -1 / 2 Cr + 1 / 2 , ω j Cr, ω j C r -1 / 2 , ω j Ω j
             Ñr,ω- j := 1 2 (1 -h(ω r-1/2 )) Ñr-1,r + h(ω r-1/2 ) Ñr-1/2,r , Ñr,ω+ j := 1 2 (1 -h(ω r+1/2 )) Ñr,r+1 + h(ω r+1/2 ) Ñr,r+1/2 , Ñr+1/2,ω- j := h(ω r+1/2 ) 2
Ñr,r+1/2 , Ñr+1/2,ω+

j := h(ω r+1/2 ) 2
Ñr+1/2,r+1 .

(2.22)

We recall that the unknowns are located at the same points as Ñdof,ω± j : points M r , S ω r+1/2 of the cell, morevover they satisfy [10][6]:

For all cells Ω j : dof Ñdof,ω,+ j + Ñdof,ω,- j = 0 (2.23)
For all dof inside the domain (Figure 12) :

j Ñdof,ω,+ j + Ñdof,ω,- j = 0 (2.24) M r+1 S ω r+1/2 M r S ω r-1/2 M r-1 Ñ r , ω - j Ñr ,ω + j Ñr + 1 / 2 ,ω - j Ñ r + 1 / 2 , ω + j Ñr- 1 / 2 ,ω + j Ñr -1 / 2 , ω - j Cr + 1 / 2 , ω j Cr, ω j C r -1 / 2
, ω j Ω j

lim ω→0 h(ω) > 0 (2.27)
As we will see thereafter, for planar degenerate conics, the property (2.27) enables us to design new finite volume schemes (when compared to original polygonal ones) for straight meshes.

More precisely, we look at the behavior of the limit (2.13) with this new nodal vectors (2.18)(2.19) (see Fig 9). With h(0) = π 2 and Ñr,r+1/2 + Ñr+1/2,r+1 = (M r M r+1 ) ⊥ after some computations, we obtain that degenerate conical normals do verify:

Cr+1/2,ω j → w→0 π 4 (M r M r+1 ) ⊥ = π 4 Ñr,r+1 , Cr,ω j → w→0 1 2 (1 -π 4 )(M r-1 M r+1 ) ⊥ = (1 -π 4 )C r j (2.28)
so that now these new normals (2.28) ( Cdof,0 j ) tends at nodes to pure polygonal ones with weight coefficient (1π/4) and at shoulder point to the non unit normal to the edge itself with weight coefficient π/4.

A continuous unification from edge to nodal fluxes: one parameter family composite normal vectors

More generaly than polygonal [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF][START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF][START_REF] Maire | A cell-centered Lagrangian scheme for two-dimensional compressible flow problems[END_REF]) and degenerate conical Cdof,0 j (and associated Ñdof,0,+ j and Ñdof,0,j ) (2.28) both valid on straight meshes, we may define a new general class of one parameter normal vector satisfying the same conditions ((2.17), and relations (2.20) (2.21)) for the geometrical first order (flat) mesh.

C r j (= 1 2 (M r-1 M r+1 ) ⊥ , N r,+ j (= 1 2 (M r M r+1 ) ⊥ ),N r,- j (= 1 2 (M r-1 M r ) ⊥ ) cf
Proposition 2.5. Composite normal on straight edge and link with quadrature formula Each edge may be parametrized by:

M planar (q) = M r + q(M r+1 -M r ), q ∈ [0, 1]. (2.29)
Let {q i } N i=1 ∈]0, 1[ (interior quadrature parameters) and {θ i } N i=1 ≥ 0 (interior quadrature weights), and suppose that:

1 - N i=1 θ i ≥ 0, (2.30)
Let q 0 = 0 and q N +1 = 1 extremity parameters, and associate weight are:

θ 0 = θ N +1 = 1 2 (1 - N i=1 θ i ), ( weight
extremities).

(2.31)

Suppose that full edge quadrature {q i , θ i } N +1 i=0 is exact for linear function:

N +1 i=0 θ i (Aq i + B) = 1 0 (Aq + B)dq, (2.32)
then the following dof normal vectors Ĉdof j : Ĉr+i/(N+1)

j = θ i (M r M r+1 ) ⊥ ,
ponctual vector evaluated at q i in (2.29) i = 1, ., N

Ĉr j = 1 2 (1 - N i=1 θ i )(M r-1 M r+1 ) ⊥ . (2.33)
do verify the volume relation (2.17), and relations (2.20) (2.21).

Proof. Firstable, note that hypothesis on quadrature formula impose the two relations:

-Exact for constant:

N +1 i=0 θ i = 1(= 1 0 1dq), (2.34)
which is verified by construction (2.31) (whatever the explicit values of each θ i ).

-Exact monomial of degree one:

N +1 i=0 θ i q i = 1 2 (= 1 0 qdq). (2.35)
thanks to q 0 = 0, and q N +1 = 1 with (2.31), now (2.35) is equivalent to:

N i=1 θ i (q i - 1 2 ) = 0. (2.36)
1. First, we look at the following volumic formula (does V j equal to volume of cell Ω j ):

V j = 1 2 dof Ĉdof j • OM dof = 1 2   r Ĉr j • OM r + r;i=1,.,N Ĉr+i/(N+1) j OM planar (q i )   (2.37)
In (2.37), we replace the explicit formula of Ĉdof j by those of (2.33).

V j = 1 2 (1 - N i=1 θ i ) r C r j • OM r + 1 2 r;i=1,.,N θ i (M r M r+1 ) ⊥ ((1 -q i )OM r + q i (OM r+1 ))(2.38) = 1 2 (1 - N i=1 θ i ) r C r j • OM r + 1 2 r (M r M r+1 ) ⊥ N i=1 θ i ((1 -q i )OM r + q i (OM r+1 )).(2.39)
Due to (2.36), we have

N i=1 θ i q i = N i=1 θ i (1 -q i ) = 1 2 N i=1 θ i ,
so that with obtain the equality:

N i=1 θ i ((1 -q i )OM r + q i (OM r+1 )) = N i=1 θ i OK [r,r+1] , with K [r,r+1] := M r + M r+1 2 .
and (2.39) rewrites:

V j = 1 2 (1 - N i=1 θ i ) r C r j • OM r + 1 2 N i=1 θ i r (M r M r+1 ) ⊥ • OK [r,r+1] , = (1 - N i=1 θ i )|Ω j | + N i=1 θ i |Ω j |, V j = |Ω j |.

Conservation in cell :

dof

Ĉdof j = r Ĉr j + r;i=1,.,N Ĉr+i/(N+1) j = (1 - N i=1 θ i ) r C r j + r;i=1,.,N θ i (M r M r+1 ) ⊥ = (1 - N i=1 θ i ) r C r j + N i=1 θ i r (M r M r+1 ) ⊥ = 0 + 0.
3. Conservation around an (interior) edge dof:

-a node r j,r∈j

Ĉr j = (1 - N i=1 θ i ) r C r j = 0.
-on an edge r + i/(N + 1) (only two contributions):

j,r+i/(N +1)∈j Ĉr+i/(N+1) j = N i=1 θ i (M r M r+1 ) ⊥ + N i=1 θ i (M r+1 M r ) ⊥ = N i=1 θ i (M r M r+1 ) ⊥ + (M r+1 M r ) ⊥ = 0
Remark 2.6. We could also theoretically deal with a general edge curve parametrization:

M γ (q), q ∈ [0, 1], with (∀q ∈ [0, 1]; | d dq M γ (q)| = 0). (2.40)
with M γ (0) = M r and M γ (1) = M r+1 . In this case, to obtain the same results than for straight meshes, we make the assumption of the existency of dof normal vectors Cdof j :

Cr+i/(N+1) j , ponctual vector evaluated at q i in (2.40) i = 1, ., N Cr j .

(2.41) (2.42)

verifying
A composite normal set is then a set for which normal at nodes extremities and at least at one interior edge point do not vanish. Both nodes and edges do contribute in the Green-Riemann formula approximation. A special (and fruitfull) case is to consider N = 1 in (2.33), because we have an unification of known results: Proposition 2.8. A continuous one parameter family composite normal vectors for straight meshes Let q 1 = 1 2 and θ 1 = θ ∈ [0, 1], we define composite nodal/edge (dof ) normal vectors by:

Ĉr+1/2 j = θ(M r M r+1 ) ⊥ = θ Ñr,r+1 , (2.43)
Ĉr j = 1 2 (1 -θ)(M r-1 M r+1 ) ⊥ = (1 -θ)C r j .
(2.44)

The decomposition (2.43)(2.44) is linked to the quadrature formula that compute the volume of a polygonal cell (2.16). Note that knowns generated normal vector can be recovered for special choices of θ (see Table ( In the following of the paper, when considering straight mesh, we will only deal with the natural limit of degenerate conical mesh ω = 0 thus taking θ = π 4 . We postpone the use of (Simpson rule) θ = 2 3 and other values of θ for an other study. Note that except for extremal values of the parameter θ = 0 and θ = 1, it generates truly composite nodal/edge normal vectors within the meaning of (2.42). We will then be able to construct a θ-scheme like between a pure normal edge (Godunov, Roe, Rusanov, HLL, etc [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]) and a pure nodal normal ( [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]) in view of formula (2.43) and (2.44). Hence, the next section is devoted to nodal extension of pure edge Approximate Riemann Solver.

Approximate Riemann Solvers (ARS) extended from edge to nodal flux

In this section, we propose a discretization of the hyperbolic system of equations (1.1). For a general presentation as well as practical guide we refer the reader to [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]. For an introduction to finite volume schemes in Lagrangian frames (total Lagrangian or updated Lagrangian) see [START_REF] Mazeran | Sur la structure mathématique et l'approximation numérique de l'hydrodynamique Lagrangienne bidimensionnelle[END_REF][START_REF] Després | Lois de Conservations Eulériennes, Lagrangiennes et Méthodes Numériques[END_REF][START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF]. We are able to design a quadrature formula at natural DOF of any cell Ω j . Let G be a scalar, vectorial or tensorial function. We have:

M r+1 G r+1 j S ω r+1/2 G r+1/2 j M r G r j S ω r-1/2 G r-1/2 j M r-1 G r-1 j Cr + 1 / 2 , ω j Cr, ω j C r -1 / 2 , ω j Ω j
∂Ωj G Nds dof G dof j Cdof,ω j = Mr∈Ωj G r j C r j , for polygons, Mr∈Ωj G r j Cr,ω j + S ω r+1/2 ∈Ωj G r+1/2 j
Cr+1/2,ω j , for arbitrary conical cells. (3.1) where G dof j is an approximation of G(x). Indeed if G is:

-a continuous function at M dof , then

G dof j = G dof G(M dof ). (3.2)
-not regular around M dof :

G dof j lim |M dof M |→0;M ∈Ωj G(M ). (3.3) M r+1 S ω r+1/2 M r S ω r-1/2 M r-1 M l M m M k Ω j Ω i Ω m Ω l Ω k
Figure 12: Two types of degree of freedom (dof) on the cell boundary ∂Ω j : the edge end points M r and mid edge points S ω r+1/2 . These are the location where the Riemann problems are solved (exactly or approximately) and the set of cells sharing this degree of freedom with Ω j . For example, from the cell Ω j point of view, at the node M r (resp. at the node S ω r+1/2 ) we need to define a numerical flux

G r j (U Ωm , U Ω l , U Ω k , U Ωj , U Ωi ) (resp. a flux G r+1/2 j (U Ωj , U Ω k )).
Integration of (1.1) in each cell Ω j :

Ωj ∂ t U + div(F (U ))dv = 0, (3.4)
gives by Green-Riemann formula (in Eulerian frame):

∂ t Ωj U (t, x)dv + ∂Ωj F (U ) Nds = 0, (3.5)
Denoting by U j (t) := 1 |Ωj | Ωj U (t, x)dv, the finite volume unknown, we approximate (3.5) with a quadrature formula involving dof of cell boundary using (3.1) (see Figures 11 and12):

d dt U j (t) + 1 |Ω j | dof ∈∂Ωj G dof j (U k1(dof ) (t), .., U km(dof ) (t), Cdof k1 , .., Cdof km ) • Cdof j = 0. (3.6)
A generic dof is linked to: nsupp(dof ) cells noted m in (3.6) for short.

(3.7)

We will use the following connectivity at dof:

k i (dof )
is the global index cell of the i'th local index linked to dof, i = 1, .., nsupp(dof ).

(3.8)

The quantities U ki(dof ) (resp. Cdof ki ) are point-wise approximation of exact solution (resp. normal vector) at M dof inside the local cell k i (where k i := k i (dof ) if no ambiguity). Remark 3.1. An idealistic approach would be to consider multi-dimensional Riemann problem around each dof (M r /S ω r+1/2 ) (see Figure [START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF]):

∂ t U (t, x) + div(F (U (t, x))) = 0, U (0, x) = U i , if x ∈ Ω i . (3.9)
in the spirit of [START_REF] Balsara | A two-dimensional HLLC Riemann solver for conservation laws: Application to Euler and magnetohydrodynamic flows[END_REF] for a genuinely multi-dimensional extension of HLL type schemes. In the present study, we will not try to find directly truly two-dimensional Riemann solvers for (3.9). We only deal with nodal extension of existing edge Approximated Riemann Solver that are consistent and locally conservative around each dof in such a way that our approach reduces to original edge formula when the degree of freedom is located at the mid edge (shoulder) point. To this end, in the next section, we propose an algebraic approach that naturally enable composite (or mixted) node/edge fluxes.

Remark 3.2. Except for vectorial linear advection equations, we will mainly focus on the nonlinear Euler system of gas dynamics, governing the conservation laws of density, momentum and total energy:

U = (ρ,
F (U ) = (F 1 (U ) F 2 (U )) =     ρu 1 ρu 2 ρu 1 u 1 + P ρu 1 u 2 ρu 2 u 1 ρu 2 u 2 + P u 1 (ρE + P ) u 2 (ρE + P )     , U = (u 1 , u 2 ), (3.13)
For our purpose the pressure P is considered as P (ρ, ) and will obey either a perfect gas EOS

P = (γ -1)ρ , (γ > 1 : adiabatic constant) (3.14)
or a Mie-Grüneisen type EOS (see [START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF]):

P = P 0 (τ ) + Γ 0 τ 0 ( -0 (τ )), (τ = 1 ρ ). (3.15)
Here Γ 0 and τ 0 are respectively the Grüneisen coefficient and the covolume of the reference state.

Definition of numerical nodal flux

The nodal flux type function G dof j (with value in R n×d ) expressed at degree of freedom dof (seen as quadrature points) of the cell j (see Figure 12) must fulfill two main standard requirements:

1. Consistency: ∀Ω j , ∀M dof ∈ Ω j G dof j (U , .., U , Cdof k1 , .., Cdof km ) • N = F (U ) • N (3.16) 2. Local conservation around dof: ∀M dof ∈ Ω j;dof ∈Ωj G dof j (U k1(dof ) , .., U km(dof ) , Cdof k1 , .., Cdof km ) • Cdof j = 0.
(3.17)

The relation (3.17) means that all components of G are -a priori-discontinuous at position of a global dof (there are jumps between cells containing the same dof ). This relation comes from:

-the discrete balance laws (3.6), and -the global conservation:

d dt j |Ω j |U j (t) = 0 =⇒ j dof ∈Ωj G dof j (U k1(dof ) , .., U km(dof ) , Cdof k1 , .., Cdof km ) • Cdof j = 0. (3.18)
so that local conservation (3.17) is just obtained as a sufficient condition to obtain (3.18) when swapping the two sums in j and dof . Until the end of the section for ease of reading, we omit the dependance on Cdof k1 , .., Cdof km in G (3.18). 3. Monotonicity (only for scalar case):

Function U j → G dof j (U k1(dof ) , .., U km(dof ) ) • Cdof j is increasing, ∀k = j function U k → G dof j (U k1(dof ) , .., U km(dof ) ) • Cdof j is decreasing. ( 3
G r+1/2 j (U j , U k ) • Cr+1/2,ω j + G r+1/2 k (U k , U j ) • Cr+1/2,ω k = 0. (3.20)
Thanks to geometrical property Cr+1/2,ω j + Cr+1/2,ω k = 0, we recover classical numerical edge flux conservation condition, and the continuity of (normal component) of numerical flux at this point. Let Njk be the unit normal going from j to k, and we choose Φ(U j , U k , Njk ) as a numerical edge flux (Roe, VFFC, Rusanov, HLL type variants, etc). The common properties of local conservation of such standard edge schemes is:

Φ(U j , U k , Njk ) + Φ(U k , U j , -Njk ) = 0. (3.21)
We also notice that this Njk (the unit normal of edge between j and k cells) and Ñr+1/2,ω

G r+1/2 j (U j , U k ) • Cr+1/2,ω j := Φ(U j , U k , Cr+1/2,ω j ) (3.22)
so that we must recover EXACTLY a given standard edge flux for all mid edge points S ω r+1/2 . Therefore, choosing a numerical edge flux solves partially our problem of finding general numerical flux fulfiling condition (3.16)-(3.17).

In the next two sections, our aim is to extend a known edge Approximate Riemann Solver (ARS) to nodal type of dof and fulfilling the ad hoc consistency (3.16) and local conservation (3.17).

Flux based schemes (Lax-Wendroff, VFFC, Roe, Flux Vector Splitting)

We recall the notion of flux schemes [START_REF] Ghidaglia | Flux schemes for solving nonlinear systems of conservation laws[END_REF] for edge based numerical fluxes Φ(U j , U k , Njk ) in (3.21) (3.22). Definition 3.4. Edge Flux schemes [START_REF] Ghidaglia | Flux schemes for solving nonlinear systems of conservation laws[END_REF] Let U ({Ω j , U j }, {Ω k , U k }, Njk ) be a R n×n matrix also called an "upwind matrix", with property:

U ({Ω j , U j }, {Ω k , U k }, Njk ) + U ({Ω k , U k }, {Ω j , U j }, -Njk ) = 0.
(3.23)

An edge flux scheme writes:

Φ(U j , U k , Njk ) = 1 2 (F (U j ) + F (U k )) • Njk -U ({Ω j , U j }, {Ω k , U k }, Njk )(F (U k ) -F (U j )) • Njk , (3.24)
Remark 3.5.

-the balance condition on tensor (3.23) ensure that (3.24) verify local conservation (3.21).

-in practice (see exemple below), an average µ is taken between the two states U j , U k (and eventually with geometrical cell property) to evaluate the "upwind" matrix:

U ({Ω j , U j }, {Ω k , U k }, Njk ) = U (µ({Ω j , U j }, {Ω k , U k }), Njk ) (3.25)
Example of average formula µ(., .):

µ({Ω j , U j }, {Ω k , U k }) = 1 2 (U j + U k ), µ({Ω j , U j }, {Ω k , U k }) = |Ωj |U j +|Ω k |U k |Ωj |+|Ω k | , ..etc.. (3.26)
In the following, for simplicity, we denote by µ(U j , .) even if there is a dependency of cell geometry µ({Ω j , U j }, .).

-Lax-Wendroff (see [START_REF] Lax | Systems of conservation laws[END_REF])

In case of Lax-Wendroff scheme, one possibility:

U (µ(U j , U k ), Njk ) = λJ(µ(U j , U k ), Njk ), λ = ∆t h . (3.27)
where µ(U j , U k ) may be linked either to Roe scheme below or (3.26).

-VFFC [START_REF] Ghidaglia | Une méthode volumes finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF] In case of VFFC ("volumes finis à flux caractéristique"), U (upwind sign matrix) in (3.24) is given by:

U (µ{U j , U k }, Njk ) = sign(J(µ{U j , U k }, Njk )) U (µ{U j , U k }, Njk ) = R(µ{U j , U k }, Njk )sign(Diag(λ i (µ{U j , U k }, Njk )))L(µ{U j , U k }, Njk ), with sign(x) =    1 x > 0 -1 x < 0 0 x = 0 (3.28)
-Roe Scheme (see [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors and Difference Schemes[END_REF]) In case of Roe scheme, µ{U j , U k } is the Roe mean value µ Roe (U j , U k ). In the case of mono-material Euler system of a perfect gas:

µ Roe (U j , U k ) =            ρ Roe = √ ρ j ρ k , U Roe = 1 √ ρj + √ ρ k ( √ ρ j U j + √ ρ k U k ), H Roe = 1 √ ρj + √ ρ k ( √ ρ j H j + √ ρ k H k ) (massic enthalpy: H = E + P/ρ), c Roe = (γ -1)(H Roe -1 2 |U Roe | 2 ) (sound speed).
(3.29)

and the upwind matrix U is defined by the following. If there exists a matrix A Roe (U j , U k ) with real eigenvalues and a corresponding complete set of eigenvectors in R n with

A Roe (U , U ) = J(U , ξ), ξ ∈ R 2 ; |ξ| = 1, A Roe (U , V )(V -U ) = (F (V ) -F (U )).ξ, (3.30)
the upwind sign matrix is finally:

U (U j , U k , Njk ) = sign(J(µ Roe (U j , U k ), Njk )) = R(µ Roe (U j , U k ), Njk )sign(Diag(λ i (µ Roe (U j , U k ), Njk )))L(µ Roe (U j , U k ), Njk ) (3.31)
As it is well known, under this flux form that VFFC scheme has the same structure as Roe scheme (when the latter exists), only the average state may be different.

-Flux vector splitting [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] (Chapter 8, p 285) A flux vector splitting scheme is such that:

Φ(U j , U k , Njk ) = (F + (U j ) + F -(U k )) • Njk . (3.32)
where F ± stands for a splitting method (F = F + + F -) (see Steger and Warming, Van Leer, .. ). F ± are constructed so that the resulting Jacobian J + (resp J -) has only positive (resp negative) eigenvalues (in practice it is used when continuous flux is s-homogeneous F (λU ) = λ s F (U ), for example in the case of linear advection or Euler system of gas dynamics (polytropic)).

Before giving some generalizations of edge flux schemes (3.24) to nodal flux schemes, we wish to emphasize the following: in order to obtain (3.17), the main problem is to give an extension of (3.23) at each dof (a node or a mid edge). The connectivity may be not uniform at each node and the number of linked cell may be arbitrary. Note that as for edge scheme, we will consider upwind flux matrix with an evaluation at some average state (3.25)(3.26)(3.29).

Definition 3.6. Flux based nodal schemes We give thereafter two sub-classes of flux based nodal scheme.

1. An all-state average nodal scheme: (called Version 1)

G dof j (U k1(dof ) , .., U km(dof ) ) • Cdof j = 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) • Cdof j -f ({α dof k } k )Λ f lux dof,k (µ({U ki(dof ) } i ), Cdof k )(F (U j ) -F (U k )) • Cdof j (3.33)
Here, µ({U ki(dof ) } i ) is a mean value of all states around the dof (depending eventually on geometric information of each cell Ω ki(dof ) ).

2. A two-state average nodal scheme: (called Version 2)

G dof j (U k1(dof ) , .., U km(dof ) ) • Cdof j = 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) • Cdof j -f ({α dof k } k )Λ f lux dof,k (µ(U j , U k ), Cdof k )(F (U j ) -F (U k )) • Cdof j (3.34)
Here, µ(U j , U k ) is a mean value of the two states (depending eventually on geometric information in each cell Ω j , Ω k ). In this case, the map µ is supposed to be symmetric, in the sense that:

µ(U j , U k ) = µ(U k , U j ) (3.35)
We give now a common assumption for (3.33) and (3.34):

Definition 3.7.

-The upwind flux matrices Λ f lux dof,k in (3.33) and (3.34) satisfy the tensorial local balance law (TLBL) if for a choice of averaged value : .36) this means that each tensor

µ({U ki(dof ) } ni i=1 ) (with ni = 2 (3.34) or ni = nsupp(dof ) (3.33)): TLBL condition : ∀dof ∈ Ω, k;dof ∈k Λ f lux dof,k (µ({U ki(dof ) } ni i=1 ), Cdof k ) = 0 R n×n ∀ U ki(dof ) ∈ A adm . ( 3
Λ f lux dof,k (µ({U ki(dof ) } ni i=1 ), Cdof k ) is defined at M dof in cell k that contains it.
-∀k, α dof k ≥ 0 and the function f ({α dof k } k ) does not depend on k (but on all α dof k , for ex. a symmetrical mean) and is supposed to be positive (e.g.

α dof k = 1 and f ({α dof k } k ) = 1). Remark 3.8.
-It should be noted that the constraint (3.36) is based both upon the physical problem (data and flux : linear or not) AND the mesh property (connectivity and both straight or curved mesh). Obtaining the (TLBL) is not a trivial task.

-Note that for the shoulder point the two numerical fluxes (3.33) and (3.34) coincide when using a common mean value formula µ(.). Hence, they give exactly the same classical generic edge formula (3.24): indeed the correspondings numerical fluxes writes :

G r+1/2 j (U j , U k ) • Cr+1/2 j = 1 2 (F (U j ) + F (U k )) • Cr+1/2 j -Λ f lux dof,k (µ(U j , U k ), Cr+1/2 k )(F (U j ) -F (U k )) • Cr+1/2 j , = 1 2 (F (U j ) + F (U k )) • Cr+1/2 j + Λ f lux dof,j (µ(U j , U k ), Cr+1/2 j )(F (U j ) -F (U k )) • Cr+1/2 j . using Λ f lux dof,j (µ(U j , U k ), Cr+1/2 j ) = U (µ(U j , U k ), Cr+1/2 j | Cr+1/2 j |
) gives the result (3.24) after normalisation by

| Cr+1/2 j |.
In the following, we first admit that the set of tensors satisfying the condition (3.36) is not empty, then we show how to build some of them. Proof. The consistency is obvious. For the local conservation (3.17), we first check for (3.33) (Version 1):

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) • Cdof j -f ({α dof k } k )Λ f lux dof,k (µ({U ki(dof ) } i ), Cdof k )(F (U j ) -F (U k )) • Cdof j (3.37)
In a first part, we group the same terms F (U k ) in (3.37):

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) • Cdof j + f ({α dof k } k )Λ f lux dof,k (µ({U ki(dof ) } i ), Cdof k )F (U k ) • Cdof j = 1 k;dof ∈k α dof k   k;dof ∈k α dof k F (U k ) • j;dof ∈j Cdof j + k;dof ∈k f ({α dof k } k )Λ f lux dof,k (µ({U ki(dof ) } i ), Cdof k )F (U k ) • j;dof ∈j Cdof j   .
The two sums over j is zero thanks to the geometric relation (2.21).

The second part including only F (U j ) term in (3.37) writes:

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k f ({α dof k } k )Λ f lux dof,k (µ({U ki(dof ) } i ), Cdof k )F (U j ) • Cdof j = f ({α dof k } k ) k;dof ∈k α dof k k;dof ∈k Λ f lux dof,k (µ({U ki(dof ) } i ), Cdof k ) j;dof ∈j F (U j ) • Cdof j
here, the mean µ({U ki(dof ) } i ) is independent of j, and so we use (3.36) to conclude. Now we check for (3.34) (Version 2):

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) -f ({α dof k } k )Λ f lux dof,k (µ(U j , U k ), Cdof k )(F (U j ) -F (U k )) • Cdof j (3.38)
We isolate the first term in (3.38):

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) • Cdof j = 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) • j;dof ∈j Cdof j
which is zero thanks to the geometric relation (2.21). The second part in (3.38):

f ({α dof k } k ) k;dof ∈k α dof k k;dof ∈k j;dof ∈j Λ f lux dof,k (µ(U j , U k ), Cdof k )(F (U j ) -F (U k )) • Cdof j (3.39)
To conclude that (3.39) vanish, we use the following Lemma 3.10 with A l = F (U l ).

Lemma 3.10. Let a set of tensors {Λ k } k=1,..,ni ∈ R n×n verifying TLBL (3.36), and { Cdof j } j=1,..,ni ∈ R d verifying (2.21), then for any matrix field {A l } l=1,..,ni ∈ R n×d , we have:

k j Λ k (A k -A j ) Cdof j = 0. (3.40)
Proof. We slice the difference, so that (3.40) rewrites:

k,j k =j Λ k A k Cdof j - k,j k =j Λ k A j Cdof j k Λ k A k j j =k Cdof j - k k =j Λ k j A j Cdof j (3.41)
now using first (3.36), we then have k k =j Λ k = -Λ j so that (3.41) is also:

k Λ k A k j j =k Cdof j + j Λ j A j Cdof j and finally k Λ k A k ( j j =k Cdof j + Cdof k ).
Hence, thanks to (2.21), we have:

k Λ k A k j Cdof j = 0.
Remark 3.11. We focus on two major dependencies with respect to Cdof k in (3.36): a linear (e.g. Lax-Wendroff ) or nonlinear way (linked to sign function like VFFC). If not otherwise stated, we will take in (3.33) and (3.34):

α dof k = 1 (so that k α dof k = nsupp(dof )) and f ({α dof k } k ) = 1.
(3.42)

Extension of (linear) Lax-Wendroff to nodal flux

Lemma 3.12. In relation (3.36), suppose that

Λ f lux dof,k (., Cdof k ) is linear w.r.t Cdof k .
Then we have the local conservation property around each dof (3.17).

Proof. Obvious using (2.21).

In particular, a possible extension for the Lax-Wendroff edge scheme (3.27) is (3.34):

Λ f lux dof,k (µ(.), Cdof k ) = ∆t |Ω| dof | Cdof k |J(µ(.), Ñdof k ), (3.43)
where |Ω| dof is a measure of some "dual" cell around dof.

Extension of (nonlinear) sign type matrix to nodal flux

Unfortunately, for the sign function like in VFFC or Roe, we can not use linear dependency w.r.t normal corner vectors. The extension of local conservation property is more tricky.

Lemma 3.13. Relationship with virtual opposite direction

We have the following matrix identity: considering a (normal) direction D and his opposite -D: Proof. We use the fact that

sign(J(V , D)) + sign(J(V , -D)) = 0. (3.44) M r Ω j Ω k Ω l
J(V , D) + J(V , -D) = 0. so that sign(J(V , D)) = sign(-J(V , -D)), (3.45)
we recall that taking a diagonlisable matrix: J = R diag(λ i ) L, the associated matricial sign function is:

sign(J) = R diag(sign(λ i )) L, (3.46)
and we have in (3.45):

sign(J(V , D)) = R(V , -D)sign(-Diag(λ i (V , -D)))L(V , -D) = -R(V , -D)sign(Diag(λ i (V , -D)))L(V , -D) sign(J(V , D)) = -sign(J(V , -D)). (3.47)
We give now a constructive way to build a nodal flux scheme (hence showing that TLBL (3.36) is a non empty set whatever the local connectivity is at dof). Proposition 3.14. Conservation using opposite normal vectors and sum rearrangement The following generic matrix defined at each dof in cell k

Λ f lux dof,k (µ{.}, Cdof k ): Λ f lux dof,k (µ{.}, Cdof k ) = 1 2   sign(J(µ{.}, Ñdof k )) - 1 nsupp(dof ) -1 nsupp(dof ) s=1, =k sign(J(µ{.}, Ñdof s ))   (3.48)
satisfies condition (3.36).

Proof. We first sum the relation (3.44) for each corner normal Cdof k coming from cells around: the idea is to get a formula without the index k using a sum rearrangement. We introduce a double index real sequence θ j s :

k;dof ∈k sign(J(µ{.}, Cdof k )) + sign(J(µ{.}, -Cdof k )) = 0. ( 3
Double index real sequence θ j s :

         nsupp(dof ) s=1 θ j s = 1, nsupp(dof ) j=1 θ j s = 1,
θ j s ≥ 0 s = j and θ j j = 0.

(3.54) then we define:

Θ j := nsupp(dof ) s=1 θ j s sign(J(µ{.}, -Cdof ks )) (3.55) using (3.54), it satisfies: nsupp(dof ) j=1 Θ j = k;dof ∈k sign(J(µ{.}, -Cdof k )) (3.56)
and we recover the second part of the initial sum in (3.49). Now, choosing (for example) in (3.54):

θ j s = 1 nsupp(dof )-1 , if s = j, θ j j = 0, else. (3.57)
we can rewrite (3.49):

k;dof ∈k   sign(J(µ{.}, Cdof k )) + 1 nsupp(dof ) -1 s;s =k;dof ∈s sign(J(µ{.}, -Cdof s ))   :=2Λ f lux dof,k (µ{.}, Cdof k ) = 0. (3.58)
Now, we use the fact that:

sign(J(µ{.}, λξ)) = sign(J(µ{.}, ξ)), ∀λ > 0, (3.59)
and

sign(J(µ{.}, -ξ)) = -sign(J(µ{.}, ξ)), (3.60)
we finally deduce the expression with normalized normal corner vectors inserted in jacobian of physical system of conservation laws.

Λ f lux dof,k (µ{.}, Cdof k ) = 1 2   sign(J(µ{.}, Ñdof k )) - 1 nsupp(dof ) -1 s;s =k;dof ∈s sign(J(µ{.}, Ñdof s ))   Remark 3.15.
Equality (3.48) gives exactly classical edge fluxes (eg. VFFC, Roe, ..) at shoulder points. Indeed, in this case, because nsupp(dof ) = 2, we have only two cells j and k so that because s = k, we have s

= j ( Ñr+1/2 s = Ñr+1/2 j ) and using that Ñr+1/2 j = - Ñr+1/2 k and we recover Λ dof,k (µ{.}, Cr+1/2 k ) = sign(J(µ{.}, Ñr+1/2 k )) = -sign(J(µ{.}, Ñr+1/2 j )).
Note that using property (3.60), we finally have an expression involving only the normal Ñr+1/2 j (and no longer their opposite as in (3.58)). Formula (3.48) guaranties that the set of tensors satisfying the (TLBL) constraints (3.36) is not empty (more generaly, see the construction of double index sequence (3.54), note that it may be also extended to a matrix form relation

nsupp(dof ) s=1 Θ s j = I n = nsupp(dof ) j=1
Θ s j (I n identity matrix) and Θ j j = 0).

In view of the definition (3.33)(3.34) and (3.48) we are able to give a constructive extension from edge to nodal (or composite) known flux schemes. -the upwind dof matrices given by (3.48), -µ(U j , U k ) the Roe average (see (3.29)) between any couple (jk) around dof.

Extension of

Extension of VFFC edge scheme to nodal flux (Version 1 or Version 2)

Definition 3.17. An extension of classical VFFC edge flux scheme (3.28) to node is given by -either (3.33) (if we consider all couple of cells (jk) named Version 2) or (3.34) (if we consider all cells simultaneously named Version 1)

-the upwind dof matrices given by (3.48), -µ(.) any symmetrical mean value formula of all cells.

Remark 3.18. This VFFC nodal flux scheme is very general and can be applied to arbitrary hyperbolic problems. However, it is not entropic and it needs an entropic correction like classical edge fluxes.

Extension of arbitrary vectorial advection system

Let's define the "Heaviside" matrix of Jacobian J (it is linked to the definition of "sign" matrix (3.28) of J):

H(J) = R diag(H(λ i )) L, where H(x) =    1 x > 0, 1 2
x = 0, 0

x < 0.

(3.61) then:

2 H(J) = I + sign(J).

(3.62) Proposition 3.19. Let's consider here an "upwind" nodal scheme for (linear thereafter) hyperbolic systems under the form (3.33) or (3.34) and consider the choice Proof.

           Λ up dof,j (µ{}, Cdof j ) = nsupp(dof )H(J(µ{}, Ñdof j )) -I , if k = j, Λ up dof,k (µ{}, Cdof k ) = nsupp(dof ) I -H(J(µ{}, Ñdof j )) B k (J(µ{}, Cdof k )) -I if k = j and with B k ∈ R n×n are such that k =j B k = I.
k =j Λ up dof,k = k =j nsupp(dof ) I -H(J(µ{}, Ñdof j )) B k (J(µ{}, Cdof k )) -I , (3.64) = nsupp(dof ) I -H(J(µ{}, Ñdof j )) k =j B k (J(µ{}, Cdof k )) - k =j I, (3.65) = nsupp(dof ) I -H(J(µ{}, Ñdof j )) -(nsupp(dof ) -1)I, (3.66) = I -nsupp(dof )H(J(µ{}, Ñdof j )), (3.67) = -Λ up dof,j . (3.68)
It is worth noting that formula (3.63) can be considered as an extension of known schemes. Indeed, we have the following Lemma 3.20.

Lemma 3.20. For the advection system associated to the (splitted) Arbitrary Lagrangian Eulerian method (e.g. for Euler system of Gas Dynamics), and for each conservative component U i of vectorial unknown U = (ρ, ρU, ρE) is transported at speed a:

∂ t U i + div(a U i ) = 0, (3.69)
we recover a nodal scheme [START_REF] Paul | Etude d'un schéma numérique volumes finis avec bilan des flux aux coins pour l'ALE[END_REF] (for polygons) and [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF] (for arbitrary conical cells) as special case of (3.63). More precisely, let [START_REF] Paul | Etude d'un schéma numérique volumes finis avec bilan des flux aux coins pour l'ALE[END_REF][START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF] writes:

V + j (dof ) = k ∈ V (j); a dof • Cdof k > 0 , V - j (dof ) = k ∈ V (j); a dof • Cdof k < 0 , the nu- merical nodal scheme in
G dof j (U k1(dof ) , .., U km(dof ) ) • Cdof j =          a dof • Cdof j U j if a dof • Cdof j > 0, a dof • Cdof j k∈V + j (dof ) a dof • Cdof k U k k∈V + j (dof ) a dof • Cdof k else. ( 3
H(a dof • Cdof j ) a dof • Cdof j U j + (1 -H(a dof • Cdof j )) a dof • Cdof j k =j β k U k , where (3.71)
β k =      a dof • Cdof k k∈V + j (dof ) a dof • Cdof k if k ∈ V + j (dof ), 0 else.

Viscosity based schemes (Rusanov, HLL, Roe, VF-Roe)

We briefly recall such edges schemes:

-Rusanov edge flux (see [START_REF] Rusanov | On difference schemes of third order accuracy for nonlinear hyperbolic systems[END_REF])

Φ(U j , U k , Njk ) = 1 2 (F (U j ) + F(U k )) • Njk - λ jk 2 (U k -U j ), where
λ jk ≥ max(ρ(J(U j , Njk )), ρ(J(U k , Nkj ))).
(3.73) and ρ(J(U , Njk )) is the spectral radius of Jacobian matrix J:

ρ(J(U , Njk )) = max i=1..n (|λ i (U , Njk )|) (= max(|λ 1 (U , Njk )|, |λ n (U , Njk )|)). (3.74)
In order to obtain the local conservation Φ(U j , U k , Njk ) + Φ(U k , U j , -Njk ) = 0, note that λ jk must verify:

λ jk = λ kj . (3.75)
-Roe flux(see [START_REF] Roe | Approximate Riemann Solvers, Parameter Vectors and Difference Schemes[END_REF])

Φ(U j , U k , Njk ) = 1 2 (F (U j ) + F (U k )) • Njk -1 2 |A Roe (µ Roe (U j , U k ))|(U k -U j ) (3.76)
where the Roe matrix and the Roe state value are given (3.30) (3.29).

-VF-Roe flux (Gallouet-Masella [START_REF] Gallouet | Un schéma de Godounov approché[END_REF]) It is based on a variable change V = φ(U ), if we note T the associated matrix change of variable ( T = ∇φ):

Φ(V j , V k , Njk ) = 1 2 (F (V j ) + F (V k )) • Njk -1 2 | T (µ(V j , V k )) J(µ(V j , V k ), Njk ) T -1 (µ(V j , V k ))|(V k -V j ), µ(V j , V k ) = Vj +V k 2 . (3.77)
Making the update of V and going back to initial variable U = φ -1 (V).

-HLL edge flux schemes (see [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF][START_REF] Batten | On the choice of wavespeeds for the HLLC Riemann solver[END_REF]) Let S min jk ( Njk ) and S max jk ( Njk ) be some estimation of minimal and maximal wave speed of the system along the normal direction of edge shared by Ω j and Ω k :

           S min jk ( Njk ) := min(λ 1 (U j , Njk ), λ 1 (U k , Njk ), λ 1 (U roe , Njk )), S max jk ( Njk ) := max(λ n (U j , Njk ), λ n (U k , Njk ), λ n (U roe , Njk )),
U roe (U j , U k ) the Roe mean value of U j and U k (used when needed and when available for the underlying physical system).

Φ(U j , U k , Njk ) =        F (U j ) • Njk , if S min jk > 0, 1 S max jk -S min jk (S max jk F (U j ) -S min jk F (U k )) • Njk - S min jk S max jk S max jk -S min jk (U j -U k ), if S min jk ≤ 0 ≤ S max jk , F (U k ) • Njk , if S max jk < 0. (3.79)
We recall that the estimated extremal speeds in (3.78), satisfy:

S min kj (-Njk ) = -S max jk ( Njk ), S max kj (-Njk ) = -S min jk ( Njk ).

(3.80) so that the numerical flux (3.79) is indeed locally conservative:

Φ(U j , U k , Njk ) + Φ(U k , U j , -Njk ) = 0.
Other variants of HLL exist (HLLC, HLLE, ..) and need some slight modification.

Thereafter, we give a definition of nodal extension to previous standard edge viscosity schemes: Definition 3.21. Viscosity based nodal schemes Like the nodal flux schemes (3.2), we consider two sub-classes:

1. An all-state average nodal scheme: (called Version 1)

G dof j (U k1(dof ) , .., U km(dof ) ) • Cdof j = 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) Cdof j + f ({α dof k } k )Λ visc dof ({U ki(dof ) , Cdof ki } i )(U j -U k ) (3.81)
2. A two-state average nodal scheme: (called Version 2)

G dof j (U k1(dof ) , .., U km(dof ) ) • Cdof j = 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) Cdof j + f ({α dof k } k )Λ visc kj (g({U k , Cdof k }, {U j , Cdof j }))(U j -U k ) (3.82)
where the function g : Proof. The consistency is obvious. For (3.81), we prove conservation as follows:

R n+d × R n+d → R n+d × R n+d in (3.82) is such that: g(X, Y ) = g(Y, X) so that Λ visc kj = Λ visc jk ( 3
j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) Cdof j + f ({α dof k } k )Λ visc dof ({U ki(dof ) , Cdof ki } i )(U j -U k ) (3.84)
The first term vanishes by geometric relation (2.21). And the second part:

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k f ({α dof k } k )Λ visc dof ({U ki(dof ) , Cdof ki } i )(U j -U k ) = f ({α dof k } k ) k;dof ∈k α dof k Λ visc dof ({U ki(dof ) , Cdof ki } i ) j;dof ∈j k;dof ∈k (U j -U k ) = 0.
For (3.82), we prove conservation as follows:

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k α dof k F (U k ) Cdof j + f ({α dof k } k )Λ visc kj (g({U k , Cdof k }, {U j , Cdof j }))(U j -U k ) (3.85)
The first part is the same in (3.84), and the second part writes:

j;dof ∈j 1 k;dof ∈k α dof k k;dof ∈k f ({α dof k } k )Λ visc kj (g({U k , Cdof k }, {U j , Cdof j }))(U j -U k ) = f ({α dof k } k ) k;dof ∈k α dof k j;dof ∈j k;dof ∈k Λ visc kj (g({U k , Cdof k }, {U j , Cdof j }))(U j -U k )
By exchanging the sums, we note that with (3.83) we have the relation:

k;dof ∈k j;dof ∈j G kj = k;dof ∈k j;dof ∈j -G kj = 0.
In practice the viscous matrices (3.81) and (3.82) are supposed to be positive (so that discrete approximation is really a viscous approximation of the conservation laws):

∀v ∈ R n , t v Λ visc {.} v ≥ 0. (3.86)
In view of the definition (3.81)(3.82), we are able to give a constructive extension from edge to nodal (or composite) known viscosity schemes.

Extension of Rusanov scheme to nodal flux

Definition 3.23. An extension of classical Rusanov edge scheme (3.73) to node is given by -For (3.81), we can take Λ visc dof = µ dof I, where µ dof is proportional to the local maximum wave speed of the system in all corner normal (and their opposite) directions from dof:

µ dof := max k;dof ∈k max l;dof ∈l ρ(J(U l , Ñdof k )) | Cdof k | (3.87) or µ dof := max k;dof ∈k ρ(J(µ({U l } l ), Ñdof k )) | Cdof k | (3.88)
-For (3.82), we can take Λ visc kj = µ dof kj I where:

µ dof kj := max ρ(J(U k , Ñdof k )) | Cdof k |, ρ(J(U j , Ñdof k )) | Cdof k |, (3.89) ρ(J(U j , Ñdof j )) | Cdof j |, ρ(J(U k , Ñdof j )) | Cdof j | or µ dof kj := max ρ(J(µ(U k , U j ), Ñdof j )) | Cdof j |, ρ(J(µ(U k , U j ), Ñdof k )) | Cdof k | (3.90)
Note that both definition (3.87) and (3.89) (for the edge shoulder point) are equal to formula (3.73) (3.75).

Remark 3.24. In practice, we sometimes need to compute the maximum speed of 1d Riemann problem associated to each two by two Riemann problem on direction Ñdof j , Guermond-Popov [START_REF] Guermond | Second-order invariant domain preserving approximation of the Euler equations using convex limiting[END_REF] do propose a practical computation for perfect gas EOS. -the viscosity matrix given by :

Λ visc kj (g({U k , Cdof k }, {U j , Cdof j })) = 1 2 (| Cdof k ||J(µ Roe (U j , U k ), Ñdof k )| + | Cdof j ||J(µ Roe (U j , U k ), Ñdof j )|) (3.91)
where µ Roe (U j , U k ) is the Roe average (see (3.29)) between any couple (jk) around dof. 1. For Version 1:

Extension of

Λ visc dof ({V ki(dof ) , Cdof ki }) = 1 nsupp(dof ) ( k;dof ∈k | Cdof k || Tk J(µ({V ki(dof ) } i ), Ñdof k ) T -1 k |) (3.92)
where Tk in (3.92) is equal to:

Tk = T (µ({V ki(dof ) } i ), Ñdof k ). (3.93)
2. For Version 2:

Λ visc kj (g({V k , Cdof k }, {V j , Cdof j })) = (3.94) = 1 2 (| Cdof k || Tkj J(µ(V j , V k ), Ñdof k ) T -1 kj | + | Cdof j || Tjk J(µ(V j , V k ), Ñdof j ) T -1 jk |)
where Tjk in (3.94) is equal to:

Tjk = T (µ(V j , V k ), Ñdof j ). (3.95)
Let summarise some remarks:

Remark 3.27.

-Stencil and multiple cell dependencies. It is worth noticing that in view of Figure 12 and the above Definitions (3.6) and (3.21), we need to define different type of neighborood of generic cell Ω j (for example cell sharing a common edge with Ω j also appear in cell sharing a common node). Definition 3.28.

V e j = {k ∈ M (k = j); ∃e * ∈ Ω j ; e * ∈ Ω k } (3.96)

V n j = {k ∈ M (k = j); ∃M r ∈ Ω j ; M r ∈ Ω k } (3.97)
We can define the corner neighbor cell of Ω j :

V c j = V n j \ V e j .
(3.98)

M r+1 M r M r-1 S ω r+1/2 N ji Ω j Ω i Ω k M r+1 M r M r-1 S ω r+1/2 N ji Ω j Ω i Ω k
M r+1 M r S ω r+1/2 Ω j Ω k M r+1 M r S ω r+1/2 Ω j Ω k

Some properties extended to nodal dof flux

First of all, let the first order in space and time (explicit) of equation semi-discrete (3.6):

U n+1 j = U n j - ∆t |Ω j | dof ∈j G dof j (U n k1(dof ) , .., U n km(dof ) ) • Cdof j , (3.99)
where we have omitted the dependency of all geometrical vector in (3.99) for ease of reading.

Note here, that the situation is completely different from the standard edge flux view point. Due to a different number of arguments at each dof (due to the fact the cell support set has a different cardinal), a single numerical flux function Φ(., ., N ) depending on only two states is no more valid. And we cannot extend easily for example the notion of monotonicity.

Whatever, we propose an indexing of the unknowns based on neighborhood kind see remark 3.27 and see Figure 15.

Without loss of generality, we suppose that for each dof we rewrite the dependency of variable according to the type of neighborood and the identity V n j = V e j ⊕ V c j . We would consider a local form of these relation:

V e j (dof ) = {k ∈ V e j ; dof ∈ k} (3.100) V n j (dof ) = {k ∈ V n j ; dof ∈ k} (3.101) V c j (dof ) = {k ∈ V c j ; dof ∈ k} (3.102)
We suppose that card(V e j (dof

)) + card(V c j (dof )) = nsupp(dof ) -1 (3.103)
and now the nodal dof flux of (3.99) may write:

G dof j (U j , {U } l∈V e j (dof ) , {U } l∈V c j (dof ) ) (3.104)
For example, if the degree of freedom dof is a shoulder point V c j (dof ) = ∅ and V e j (dof ) = {k(dof )} so that (3.104) writes:

G dof j (U j , U k ) (3.105)
and for a vertex we have V e j (dof ) = {k 1 (dof ), k 2 (dof )} and V c j (dof ) is not fixed:

G dof j (U j , U k1 , U k2 , {U } l∈V c j (dof ) ) (3.106)
Note that for a mesh containing only nodes linked with exactly 3 cells V c j (dof ) ≡ ∅ (for example Voronoy type mesh).

In the following, we will give some qualitative results based on the scheme (3.99)(3.104) and we make some adaptation of edge flux type schemes [START_REF] Seguin | Méthode de volumes finis pour les fluides compressibles[END_REF][25] [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF] to our composite fluxes.

Because of the high dependency on the mesh, we will make some assumptions on it.

Hypothesis 3.29. The mesh must fullfill the following hypothesis, let h := sup j diam Ω j the maximum cell diameter over the mesh.

-Volumic

κ v |∂Ω j |h ≤ h d ≤ 1 κ v |Ω j |, (3.107)
-Curvilinear bound

κ e h ≤ | Cdof j | ≤ 1 κ e h or κ e |∂Ω j | ≤ dof ∈j | Cdof j | ≤ 1 κ e |∂Ω j |. (3.108)
-Shape regularity:

∀Ω j , ∃M * j ∈ Ωj , such that ∀dof ∈ ∂Ω j , Cdof j • (M dof -M * j ) > 0.
(3.109)

Monotonicity condition on scalar problems

In a more compact form (3.99) (3.104):

U n+1 j = H(U n j , {U n } l e ∈V e j , {U n } l c ∈V c j ) (3.110)
where the set are defined in (3.96) (3.98). In order to obtain monotonic scheme we argue the following: Definition 3.30. A nodal scheme (3.110) is called monotonic: -Individual dof monotonicity:

∀l ∈ j ∪ V e j ∪ V c j , ∂ U n l H(U n j , {U n } l e ∈V e j , {U n } l c ∈V c j ) ≥ 0. ( 3
the function U j → G dof j (U n j , {U n } l e ∈V e j (dof c ) , {U n } l c ∈V c j (dof c ) ) • Cdof j is increasing, ∀l = j ∈ V e j (dof c ); V c j (dof c ) : the function U l → G dof j (U n j , {U n } l e ∈V e j (dof c ) , {U n } l c ∈V c j (dof c ) ) • Cdof j is decreasing. (3.112) -Time step restriction: If L is the maximum of all dof in Ω j of the Lipschitz constant of G dof j (., {U n } l e ∈V e j (dof c ) , {U n } l c ∈V c j (dof c
) ) (and depends on time values at t n ).

∆t ≤ min j |Ω j | L dof ∈j | Cdof j | . ( 3

.113)

Proof. For corner neighbor cell:

∀k c ∈ V c j ; ∃! dof c , k c ∈ V c j (dof c ), ∂ U n k c H(U n j , {U n } l e ∈V e j , {U n } l c ∈V c j ) = -∆t |Ωj | ∂ U n k c G dof c j (U n j , {U n } l e ∈V e j (dof c ) , {U n } l c ∈V c j (dof c ) ) • Cdof c j ≥ 0, (3.114)
which means that the associated numerical flux G dof j • Cdof j is decreasing with respect to corner neighboring states of the generic cell j. For edge neighbor cell: ∀k e ∈ V e j ; ∃! dof e , k e ∈ V e j (dof e ) and ∃ dof c1 , dof c2 ; k e ∈ V c j (dof c1 ), k e ∈ V c j (dof c2 ),

∂ U n k e H(U n j , {U n } l e ∈V e j , {U n } l c ∈V c j ) = -∆t |Ωj | ∂ U n k e G dof c1 j (U n j , U n k e 1 , U n k e , {U n } k c ∈V c j (dof c1 ) ) • Cdof c1 j + ∂ U n k e G dof e j (U n j , U n k e ) • Cdof e j + ∂ U n k e G dof c2 j (U n j , U n k e , U n k e 2 , {U n } k c ∈V c j (dof c2 ) ) • Cdof c2 j ≥ 0.
(3.115) the sum in (3.115) involves three terms (two at the extremities and one at mid edge). With the sufficient condition (but not necessary!) that each of the numerical flux G dof j • Cdof j is decreasing with respect to edge neighboring states of the generic cell j. For cell j: 3.4.2 Invariant domain preserving (scalar or system)

∂ U n j H(U n j , {U n } l e ∈V e j , {U n } l c ∈V c j ) = 1 - ∆t |Ω j | dof ∈j ∂ U n j G dof j (U n j , {U n } l∈V e j , {U n } l∈V c j ) • Cdof
U n+1 j = U n j - ∆t |Ω j | dof ∈j G dof j (U n k1(dof ) , .., U n km(dof ) ) • Cdof j , (3.118)
where the numerical nodal flux G dof j do verify the local conservation constrain:

for all dof ∈ Ω j G dof j (U n k1(dof ) , .., U n km(dof ) ) • Cdof j = 0. (3.119)
In case of systems, we look at preserving physical properties of numerical solution (3.118):

Definition 3.32. Invariant domain preserving Let an invariant domain C ⊆ A adm , a scheme is invariant domain preserving if:

U n j ∈ C ⇒ U n+1 j ∈ C. (3.120)
where for example C = {ρ > 0, > 0} for Euler system of gas dynamics with a polytropic EOS. Using (2.20), we have dof ∈j F (U j ) Cdof j = 0, so that we can rewrite (3.99):

U n+1 j = U n j - ∆t |Ω j | dof ∈j G dof j (U n k1(dof ) , .., U n km(dof ) ) -F (U n j ) • Cdof j . (3.121)
Let some γ dof j > 0 to be precised later such that dof ∈j

γ dof j = 1. U n+1 j = dof ∈j γ dof j U n j - ∆t |Ω j |γ dof j G dof j (U n k1(dof ) , .., U n km(dof ) ) -F (U n j ) • Cdof j . (3.122)
We can define the dof's contributive part in the solution update U n+1 j in terms of individual contributions coming from nodal flux at dof: Definition 3.33. Let a dof around which we fix a unit normal direction N and σ > 0 given (σ = σ dof j (U k1(dof ) , .., U km(dof ) , N)), then we define the following function: | so that σ dof j = σ j (same for all dof in cell Ω j ) finally:

Ū dof j (U k1(dof ) , .., U km(dof ) , N, σ) = U j - 1 σ G dof j (U k1(dof ) , .., U km(dof ) ) -F (U j ) • N (3.123) S r r+1/2 S r+1 r+1/2 S r-1 r-1/2 S r r-1/2 S r-2 r-3/2 S r-1 r-3/2 S r+2 r-5/2 S r-2 r-5/2 S r+1 r+3/2 S r+2 r+3/2 M r+1 M r S ω r-1/2 M r-1 S ω r+1/2 x j S k k+1/2 S k+1 k+1/2 x k S l l+1/2 S l-1 l+1/2 x l S m m+1/2 S m-1 m+1/2 x m x i M l M m M k Ω j Ω i Ω m Ω l Ω k S r r+1/2 S r+1 r+1/2 S r-1 r-1/2 S r r-1/2 S r-2 r-3/2 S r-1 r-3/2 S r+2 r-5/2 S r-2 r-5/2 S r+1 r+3/2 S r+2 r+3/2 M r+1 M r S ω r-1/2 M r-1 S ω r+1/2 x j S k k+1/2 S k+1 k+1/2 x k S l l+1/2 S l-1 l+1/2 x l S m m+1/2 S m-1 m+1/2 x m x i M l M m M k Ω j Ω i Ω m Ω l Ω k
U n+1 j = dof ∈j | Cdof j | dof ∈j | Cdof j | Ū dof j (U n k1(dof ) , .., U n km(dof ) , Ñdof j , |Ω j | ∆t dof ∈j | Cdof j | ) ( 3
U k1(dof ) , .., U km(dof ) ∈ C ⇒ Ū dof j (U k1(dof ) , .., U km(dof ) , Ñdof j , σ dof j (U k1(dof ) , .., U km(dof ) , Ñdof j )) ∈ C. (3.126)
The property is still true for all σdof j ≥ σ dof j . Proposition 3.35. Let C be an invariant domain. If the scheme (3.118) is invariant domain preserving for C, then the associated numerical flux is invariant domain preserving (3.123) with:

σ dof j (U k1(dof ) , .., U km(dof ) , Ñdof j ) = |Ω j | ∆t| Cdof j | . ( 3

.127)

Proof. Consider the generic scheme (3.118) in a cell j, we call dof * a special dof, for all other dof, we consider that ∀dof ∈ Ω j (dof = dof * ), U k l (dof ) = U j , for all l = 1, .., nsupp(dof ).

(3.128)

By doing this, we suppose that we have separated each nodal flux from each other. Considering then each sub-cell around them and defining an associated dual control volume cell. By consistency of the numerical flux G dof j , the scheme rewrites:

U n+1 j = U n j - ∆t |Ω j |   dof =dof * ∈j F (U n j ) • Cdof j + G dof * j (U n k1(dof * ) , .., U n km(dof * ) ) • Cdof * j   (3.129)
using the property

dof =dof * Cdof j = -Cdof * j
, we obtain:

U n+1 j = U n j - ∆t |Ω j | (G dof * j (U n k1(dof * ) , .., U n km(dof * ) ) • Cdof * j -F (U n j ) • Cdof * j ) (3.130)
which is nothing but saying that U n+1 j is equal to:

Ū dof * j (U k1(dof * ) , .., U km(dof * ) , Ñdof * j , σ dof * j (U k1(dof * ) , .., U km(dof * ) , Ñdof * j )). (3.131) with σ dof * j (U k1(dof * ) , .., U km(dof * ) , Ñdof * j )) = |Ωj | ∆t| Cdof * j | . S r r+1/2 S r+1 r+1/2 S r-1 r-1/2 S r r-1/2 S r-2 r-3/2 S r-1 r-3/2 S r+2 r-5/2 S r-2 r-5/2 S r+1 r+3/2 S r+2 r+3/2 M r+1 S ω r+3/2 M r S ω r-1/2 M r-1 S ω r+1/2 M r-2 S ω r-3/2 M r+2 S ω r-5/2 x j Ω j = dof Ω dof j Ω r j Ω r+1/2 j Ω r-1/2 j Ω r+1 j Ω r+3/2 j Ω r-1 j Ω r-3/2 j Ω r-2 j Ω r-5/2 j Ω r+2 j
σ dof j (U k1(dof ) , .., U km(dof ) , Ñdof j )∆t ≤ β dof j |Ω j | | Cdof j | . (3.133)
Proof. Let's define the following approximate solution in virtual sub cell of Ω j linked to dof:

U dof,n+1 j := U n j - ∆t| Cdof j | β dof j |Ω j | G dof j (U n k1(dof ) , .., U n km(dof ) ) -F (U n j ) • Ñdof j ( 3
1 σ = ∆t| Cdof j | β dof j |Ω j | ≥ 1 σ scheme . (3.135)
We deduce that U dof,n+1 j belong to C. Moreover, the scheme rewrites:

U n+1 j = dof ∈j β dof j U dof,n+1 j (3.136)
Due to the fact that C is convex, U n+1 j in (3.136) belong also to C. Pratically, the β dof j are designed to be linked to cell properties:

β dof j = 1 ndof (j) , ndof (j
) is the number of dof in cell j or, (3.137)

β dof j = nsupp(dof ) dof ∈j nsupp(dof )
or, (3.138)

β dof j = | Cdof j | dof ∈j | Cdof j | (3.139)
We can also associate truly physical sub cells.

Proposition 3.37. Admissible geometrical sub cells With the geometric hypothesis on the mesh (3.109) |Ω dof j | the volume of the sub-cell of Ω j attached to dof (see Figure 18):

|Ω dof j | = 1 2 Cdof j • (M * j M dof ) (3.140)
and so the following choice do verify (3.132):

β dof j := |Ω dof j | |Ω j | (3.141)
The construction of geometric boundary of the sub cell defining Ω dof j (see Figure 18) rely on a correct position of intermediate nodes around each dof. More precisely, around M r let S r r-1/2 and S r r+1/2 and around S ω r+1/2 let S r r+1/2 and S r+1 r+1/2 , all these points permits to build the sub-cell Ω dof j verifying volumic equality (3.140). Remark 3.38. Case of straight meshes -For degenerate conical mesh (ω ≡ 0), taking M * j = x j (the centroïd), the intermediate nodes S r r+1/2 , S r+1 r+1/2 are given by:

S r r+1/2 = π 4 M r + (1 - π 4 )S 0 r+1/2 = S 0 r+1/2 + π 4 (M r -S 0 r+1/2 ), (3.142)
S r+1 r+1/2 = π 4 M r+1 + (1 - π 4 )S 0 r+1/2 = S 0 r+1/2 + π 4 (M r+1 -S 0 r+1/2 ). (3.143)
where S 0 r+1/2 = 1 2 (M r + M r+1 ) is also the mid planar edge, it may be rewritten:

S r r+1/2 = 1 2 (1 + π 4 )M r + 1 2 (1 - π 4 )M r+1 , (3.144)
S r+1 r+1/2 = 1 2 (1 + π 4 )M r+1 + 1 2 (1 - π 4 )M r . (3.145)
-More generaly, within the meaning of (2.43) (2.44) with a single weight quadrature parameter θ ∈ [0, 1], we have the following: Mr-1

S r r+1/2 = 1 2 (1 + θ)M r + 1 2 (1 -θ)M r+1 , (3.146)
S r+1 r+1/2 = 1 2 (1 + θ)M r+1 + 1 2 (1 -θ)M r . ( 3
S 0 r+1/2 x j S k k+1/2 S k+1 k+1/2 x k S l l+1/2 S l-1 l+1/2 x l S m m+1/2 S m-1 m+1/2 xm x i M l Mm M k Ω j Ω i Ω m Ω l Ω k S r r+1/2 S r+1 r+1/2 S r-1 r-1/2 S r r-1/2 S r-2 r-3/2 S r-1 r-3/2 S r+2 r-5/2 S r-2 r-5/2 S r+1 r+3/2 S r+2 r+3/2 Mr+1 Mr S 0 r-1/2
Mr-1

S 0 r+1/2 x j S k k+1/2 S k+1 k+1/2 x k S l l+1/2 S l-1 l+1/2 x l S m m+1/2 S m-1 m+1/2 xm x i M l Mm M k Ω j Ω i Ω m Ω l Ω k S r r+1/2 S r+1 r+1/2 S r-1 r-1/2 S r r-1/2 S r-2 r-3/2 S r-1 r-3/2 S r+2 r-5/2 S r-2 r-5/2 S r+1 r+3/2 S r+2 r+3/2 Mr+1 Mr S 0 r-1/2
Mr-1 -Link with DDFV two dimensionnal discretizations Hence, in view of Figure [START_REF] Després | Lois de Conservations Eulériennes, Lagrangiennes et Méthodes Numériques[END_REF] (Right) and Figure [START_REF] Dumbser | ADER discontinuous Galerkin schemes for aeroacoustics[END_REF] (Right), a composite finite volume scheme naturally defines both a dual diamond cell (around each edge) AND a dual nodal cell (around each node). In this way, these sub-cells construction can be used when some unknowns are defined either at edge or at node or both. This is in some sense a possible direct extension of Discrete Duality Finite Volume (DDFV) method (both on hyperbolic and elliptic equations) on straight meshes (but can also be transposed to curved conical meshes).

S 0 r+1/2 x j S k k+1/2 S k+1 k+1/2 x k S l l+1/2 S l-1 l+1/2 x l S m m+1/2 S m-1 m+1/2 xm x i M l Mm M k Ω j Ω i Ω m Ω l Ω k

Case of Rusanov composite scheme

Now, like standard edge Rusanov scheme we have a stability result concerning preservation of convex domain for discrete unknowns. Using Definition 1.1.

Proposition 3.39. The first order in space and time (explicit) Rusanov nodal scheme (3.90) is invariant domain preserving under the time step restriction:

∆t ≤ min j |Ω j | dof ∈j 2 nsupp(dof ) k;dof ∈k k =j µ dof jk (3.148) Proof. U n+1 j = U n j - ∆t |Ω j | dof ∈j 1 nsupp(dof ) k;dof ∈k F (U n k ) Cdof j + µ dof jk (U n j -U n k ). ( 3
U n+1 j = U n j - ∆t |Ω j | dof ∈j 1 nsupp(dof ) k;dof ∈k (F (U n k ) -F (U n j )) Cdof j + µ dof jk (U n j -U n k ). (3.150)
and using U n j = 2U n j -U n j :

U n+1 j =   1 - ∆t |Ω j | dof ∈j 1 nsupp(dof ) k;dof ∈k k =j 2µ dof jk    U n j + (3.151)    ∆t |Ω j | dof ∈j 1 nsupp(dof ) k;dof ∈k k =j (F (U n j ) -F (U n k )) Cdof j + µ dof jk (U n j + U n k )    .
As usual, we introduce the intermediate states:

Ū jk := 1 2 (U n j + U n k ) + 1 2µ dof jk (F (U n j ) -F (U n k )) Cdof j (3.152)
so that (3.151):

U n+1 j =   1 -∆t |Ωj | dof ∈j 1 nsupp(dof ) k;dof ∈k k =j 2µ dof jk    U n j + ∆t |Ωj | dof ∈j 1 nsupp(dof ) k;dof ∈k k =j 2µ dof jk Ū jk (3.153)
in order to obtain a convex combination of old value U n j and intermediate values Ū jk , we need to show two major things:

-Explicit constraint on time step: -Each of intermediate states Ū jk do verify Ū jk ∈ C like U n j . As in [START_REF] Guermond | Second-order invariant domain preserving approximation of the Euler equations using convex limiting[END_REF], in the spirit of original paper of [START_REF] Harten | On upstream differencing and Godunov-type schemes for hyperbolic conservation laws[END_REF]. Let U RS (ξ, U l , U r , N) be the Riemann Solver

1 -∆t |Ωj | dof ∈j 1 nsupp(dof ) k;dof ∈k k =j
-Consistency: U RS (ξ, U , U , N) = U , ∀ξ ∈ R.
-Consistency with integral formula with Riemann data:

∀ U l , U r ∈ A adm , 1 ∆x ∆x/2 -∆x/2 U RS ( x δt , U l , U r , N)dx = 1 2 (U l + U r ) - δt ∆x (F (U r ) -F (U l ))N, (3.155)
with the pseudo time step restricted to CFL 1/2 condition:

δt ≤ ∆x 2 max i=1..n |λ i (U , N)| . (3.156)
Moreover, with (3.156), the Riemann solution: Some other type of discretization such as Finite Element scheme [START_REF] Guermond | Second-order invariant domain preserving approximation of the Euler equations using convex limiting[END_REF] or Residual Distribution scheme [1] also deal with this type of Rusanov approach.

U RS ( x δt , U l , U r , N) is itself invariant domain preserving. ( 3

Numerical example of first order composite nodal/edge flux type schemes and viscosity type schemes

In this section, we check the convergence of two numerical composite conical scheme of flux type (3.33) or (3.34) and of viscosity type (3.81) or (3.82). More precisely, we use composite VFFC and composite Rusanov (with Version 2). For the first, we show that it may capture a non entropic solution (rarefaction shock in Figure 21). Hence, we need to modify it in order to make the scheme entropic. We propose a treatment based on the analogous of classical edge flux and based on a use of Rusanov entropic ones wherever an eigenvalue is changing sign for a genuinely non-linear field. For nodal numerical scheme, in order to detect if a changing sign occurs at a dof for an eigenvalue i associated to a GNL field (CSAGNL), we consider if this event occurs:

∃ i; λ i (U j , Ñdof j ) < 0 < λ i (U k , Ñdof j ) for k; dof ∈ Ω k , called CSAGN L event. (3.160)
and the j's dependancy (at dof) associated set: (3.163) Remark 3.43. Note that none of the correction (3.162) (3.163) do lead to Lipschitz continuity of numerical fluxes so that high order accuracy may be lost in the vicinity of sonic points (ie states U * such that any of eigenvalues do vanish λ i (U * , Ñdof j ) = 0). In the spirit of [START_REF] Helluy | A simple parameter-free entropy correction for approximate Riemann solvers[END_REF], if we want to obtain a continuous correction of non entropic edge flux scheme, we may define a continuous flux based numerical nodal scheme (FBNF) (see (3.33) or (3.34)) :

k CSAGN L dof (j) = { k; dof ∈ Ω k ; verifying (3.
G dof,F BN F -CC j =    G dof,F BN F j + 1 k;dof ∈k α dof k min k∈k CSAGN L dof (j) (|λ i (U k , Ñdof j )|| Cdof j |) k;dof ∈k (U j -U k ) if CSAGN L, G dof,F BN F j else. (3.164) dof,F BN F -CC j =      G dof,F BN F j + 1 k∈k CSAGN L dof (j) α dof k k∈k CSAGN L dof (j) min(|λ i (U j , Ñdof j )|, |λ i (U k , Ñdof j )|)| Cdof j |(U j -U k ) if CSAGN L, G dof,F BN F j else. (3.165
We emphasize that this correction may lead to insufficient amount of numerical viscosity to proove entropy or invariant domain preservation property.

Remark 3.44. Approximation of two dimensional test problems In the following, for genuinely two dimensionnal problems except in some special test cases, the given meshes are generated with some widely used (free) mesh generator. We can cite Gmsh (see https://gmsh.info/ or [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in preand post-processing facilities[END_REF]), or FreeFem++ (see http://www3.freefem.org/ or [START_REF] Hecht | New development in freefem++[END_REF]), or structured mesh (cartesian or polar). In many test cases, cells are only made by straight edges:

-Purely nodal flux require only evaluation of "corner" normal vector at node r C r j in each cell j (see also [START_REF] Després | Lagrangian gas dynamics in two dimensions and lagrangian systems[END_REF][START_REF] Carré | A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension[END_REF]).

-Composite nodal/edge (conical) flux require now "corner" normal vector at node r Cr,ω j but also at mid edge normal vector Cr+1/2,ω j in each cell j. Now, it remains to take into account the prescribed curvature of an edge. If you do it at mid-edge shoulder point S r+1/2 , we recall a result:

κ r+1/2 := κ(S r+1/2 ) = 8ω r+1/2 A T (M r , M r+1/2 , M r+1 ) |M r M r+1 | 3 = 8(ω r+1/2 + 1) A T (M r , S r+1/2 , M r+1 ) |M r M r+1 | 3 ) (3.166)
1. Zero curvature edge : we consider that ω r+1/2 = 0 and M r+1/2 = 1 2 (M r + M r+1 ) = S r+1/2 , so that we have a "degenerate conic edge".

2. κ r+1/2 = 0, there are three unknowns (M r+1/2 , ω r+1/2 ) and one equation, we may complete the system by also imposing the curvature at each end point. In case of circular edge (having O origin and R radius), we recall that ω r+1/2 = cos( θ 2 ) (θ is the angle M r OM r+1 ), and the control point M r+1/2 (and then S r+1/2 ) is on the mediatrice of [M r , M r+1 ], using the constrain R = |OS r+1/2 | we find an explicit position of S r+1/2 (and then thoose of M r+1/2 ).

Approximation of one dimensional test problems

Despite the formalism which only makes sense for two-dimensional problems, we want to use all this 2d-frame for the approximation of one dimensional test cases. For the computation of numerical solution of (1.1) on x 1 ∈ [0, L] (see Figure 20), we fix the second coordinate interval to x 2 ∈ [0, 1] and we consider rectangular cells of height ∆x 2 = 1. We use a symmetric boundary condition to {x 2 = 0} and {x 2 = 1}. Note that in this situation, we are forced to use degenerate conical case edge to take into account these rectangular cells.

O L 1 ∆x 2 = 1 x 1 x 2 Ω j Ω k

Entropic correction of nodal/edge flux VFFC/Roe schemes

We want to approximate the Sod shock tube with a sonic point (Cf. E. Toro [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics[END_REF]). Let's consider the Riemann problem at x 1 = 0.5 with final time tf=0.1: 23)) the numerical convergence of both first order nodal Rusanov and VFFC(-C) schemes and compare qualitatively to reference solution. Both do converge to exact solution: the arrival times associated to each waves seem correct and each level are well approxiamted. However, as expected the nodal VFFC-C flux scheme show better resolution in each waves with respect to nodal Rusanov scheme.

x1 < 0.5 :    ρg = 1 Ug = (0.75, 0) Pg = 1 x1 > 0.5 :    ρ d = 0.125 U d = (0, 0) P d = 0.1 (
4 From low to high order nonlinear composite schemes on straight or curved meshes

In this section, in view of the building of generic first order composite nodal/edge numerical flux : Definition 3.6 with (3.33)(3.34) and Definition 3.21 with (3.81) (3.82). We describe higher order versions and the limitation process for systems. We will consider volumic or massic variables. The target system under consideration is the mono-fluid Euler system either for perfect gas (convex) or Mie-Grüneisen (non-convex) EOS.

∂ t ρ + div(ρU) = 0, (4.1) ∂ t (ρU) + div(ρU ⊗ U + P I 2 ) = 0, (4.2) ∂ t (ρE) + div(ρEU + P U) = 0. (4.3)
where the pressure equation of state P (ρ, ) depends on fluid nature, and his internal energy,

= E - 1 2 |U| 2 . (4.4)
In case of perfect gas EOS, the unknowns remain in a convex set of R 4 given by:

ρ ≥ ρ min > 0, (4.5) U ∈ R 2 , (4.6) ≥ min > 0. (4.7)
In case of Mie-Grüneisen EOS, the unknowns remain in a smaller convex set (see [START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF]).

High order (up to third) spatial nodal polygonal or conical schemes

In order to obtain high order space conical nodal schemes, we first define R j (x dof , U ) a high order reconstruction of conservative variables on each conical cells Ω j , we use a Least Square type method (see [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF]). Typically, if n stands for the degree of the polynomial reconstruction R j (x dof , U ), then the order is n+1. In the second part, we use a composite nodal/edge flux scheme G dof j defined in the previous section 3. For example, we may choose a flux based scheme (3.6) or a viscosity based scheme (3.21):

d dt U j (t) + 1 |Ω j | dof ∈Ωj G dof j (R k1(dof ) (x dof , U ), .., R km(dof ) (x dof , U )) • Cdof j = 0. (4.8)
For the limitation of a quantity Q, we impose some numerical stability criteria (such as thermodynamic admissible states and admissible states compatible with the hyperbolicity of the system):

1. Positivity (scalar variable): Extrapolated values of reconstruction at dof:

min dof R j (x dof , Q) > 0 (4.9)
2. Global bound (scalar variable): Q min , Q max :

min dof R j (x dof , Q) ≥ Q min , (4.10) max dof R j (x dof , Q) ≤ Q max . (4.11)
3. Local bound in space (scalar variable) (here we note V (j) = V n (j) see notation (3.97)):

min V Q V (j) ≤ R j (x dof , Q) ≤ max V Q V (j) .
(4.12)

4. "TVB": there is an extensive literature on how to minimize spatial oscillation during the reconstruction step. A common way is to control extrapolated reconstructed values in terms of some local bound (min/max value around all the neighboring cells), hence a natural limitation is based on extending slope limiting strategy in two dimensional unstructured grids (with regards to the extension of total variation of reconstructed function from one dimensional theory). Here, we can make an adaptation of [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF].

(a) Barth/Jespersen, Dukowicz [START_REF] Barth | The design and application of upwind schemes on unstructured meshes, chapter 6[END_REF], MLP [START_REF] Park | Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids[END_REF], or min-mod for scalar (volumic or massic) variables.

(b) Vip [START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF] for vector (volumic or massic) variables.

For second, third (and possibily higher) order methods, we propose a reconstruction of conservative variables and make possible limitation on variables involved in continuous physical constraint, both volumic and massic. Note also that a TVB algorithm has also been reported [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF] in order to minimize jump of some candidate reconstruction function of finite volume based edge scheme: typical candidate reconstructions are either polynomials and tanh type functions (this latter being oscillation free).

Nonlinear -a priori-strategies for the reconstructions step

We focus on nonlinear strategies for the reconstruction of massic variables using Leibniz formula. If n + 1 is the targeted order, then for any massic quantity S, whenever P n,h j (x, ρ) > 0, we define its nonlinear reconstruction by the rational function R n,h : R n,h j (x, S) :=

P n,h j (x, ρS) P n,h j (x, ρ) . (4.13)
For the reconstruction of P n,h j (x, ρS),

P n,h j (x, ρS) = (ρS) j + (∇(ρS)) j (x -x j ) + 1 2 t (x -x j )(∇ 2 (ρS)) j (x -x j ) (4.14) +... + 1 n! (∇ n (ρS)) j (x -x j , .., x -x j )
where

(ρS) j = (ρS) j - n k=1
Cte k j (ρS) (4.15) and :

Cte k j (Q) := 1 k! 1 |Ω j | Ωj (∇ k Q) j (x -x j , .., x -x j )dv. (4.16)
Instead of dealing with the volumic variable ρS, the idea is to deal with the variations of each variable ρ and S independently and then gather all the information, more precisely, we do the following:

1. Suppose the fields ρ and S are regular, we can use the Leibniz rule on any order differentiation of the product ρS so that:

       ∇(ρS) = ρ∇S + S∇ρ, ∇ 2 ρS = ρ∇ 2 S + S∇ 2 ρ + ∇ρ ⊗ ∇S + ∇S ⊗ ∇ρ, .. = .. ∇ n ρS = ρ∇ n S + S∇ n ρ + ∇ρ∇ n-1 S + ∇S∇ n-1 ρ + .. (4.17)
2. For each variable ρ and S, choose a numerical method to compute high order terms: for example, the specific variable S in each cell Ω j : (∇S) j , (∇ 2 S) j , .. (∇ n S) j (each must be an approximation of (∇ l S)(x j )).

3. Use all the developpement (4.17) to deduce P n,h j (x, ρS) in (4.14).

For example, we apply (4.14) (4.17) to U = u 1 u 2 and E. The gradient of vector U, is noted (∇U) ij = (∂ xj u i ). Here for second order (n = 1) and third order (n = 2) approximation:

∇(ρU) = ρ∇U + U ⊗ ∇ρ, (4.18) ∇(ρE) = ρ∇E + E∇ρ, (4.19) ∇ 2 (ρU) = ρ∇ 2 u 1 u 2 + (∇ 2 ρ) u 1 u 2 + ∇ρ ⊗ ∇ u 1 u 2 + ∇ u 1 u 2 ⊗ ∇ρ, (4.20) ∇ 2 (ρE) = ρ∇ 2 E + E∇ 2 ρ + ∇ρ ⊗ ∇E + ∇E ⊗ ∇ρ. (4.21)

Direct limitation of Density, Velocity and Total massic Energy

Let consider second order reconstructions for density, momentum and volumic total energy:

P 1,h j (x, ρ) = ρ j + (∇ρ) j (x -x j ). (4.22) P 1,h j (x, ρU) = (ρU) j + (∇(ρU)) j (x -x j ), (4.23) = ρ j U j + (ρ j (∇U) j + U j ⊗ (∇ρ) j )(x -x j ).
(4.24)

P 1,h j (x, ρE) = (ρE) j + (∇(ρE)) j (x -x j ), (4.25) = ρ j E j + (ρ j (∇E) j + E j (∇ρ) j )(x -x j ). (4.26)
Using relation (4.13) for two massic variables, velocity U and total massic energy E, after simplification we obtain:

R 1,h j (x, U) = U j + ρ j P 1,h j (x, ρ) (∇U) j (x -x j ), (4.27)
R 1,h j (x, E) = E j + ρ j P 1,h j (x, ρ) (∇E) j (x -x j ). (4.28)
In order to obtain a single limitation in each cell, we define a real scalar cell number attached to the high order term for each independent variable, let α q j be this number (where q is either: density ρ, velocity U, or total massic energy E).

P 1,h j (x, ρ) = ρ j + α ρ j (∇ρ) j (x -x j ), (4.29) R 1,h j (x, U) = U j + ρ j P 1,h j (x, ρ) α U j (∇U) j (x -x j ), (4.30)
R 1,h j (x, E) = E j + ρ j P 1,h j (x, ρ) α E j (∇E) j (x -x j ). (4.31)
For third order reconstruction, we have:

P 2,h j (x, ρ) = ρj + (∇ρ) j (x -x j ) + 1 2 t (x -x j )(∇ 2 ρ) j (x -x j ), (4.32)
P 2,h j (x, ρU) = (ρU) j + (∇(ρU)) j (x -x j ) + 1 2 t (x -x j )(∇ 2 (ρU)) j (x -x j ), (4.33)
P 2,h j (x, ρE) = (ρE) j + (∇(ρE)) j (x -x j ) + 1 2 t (x -x j )(∇ 2 (ρE)) j (x -x j ). (4.34)
Using relation (4.13) and Leibniz formula (4.18)-(4.21) for two massic variables, velocity U and total massic energy E, after simplification we obtain:

R 2,h j (x, U) = U j + -CteRest 2 (ρU) j + ρ j (∇U) j (x -x j ) + 1 2 t (x -x j )(Hess rest (ρU)) j (x -x j ) P 2,h j (x, ρ) ,(4.35) R 2,h j (x, E) = E j + -CteRest 2 (ρE) j + ρ j (∇E) j (x -x j ) + 1 2 t (x -x j )(Hess rest (ρE)) j (x -x j ) P 2,h j (x, ρ) . (4.36)
with notations (4.16) for constant terms Cte k (Q) j :

CteRest 2 (ρU) j := Cte 2 (ρU) j -U j Cte 2 (ρ) j , (4.37) CteRest 2 (ρE) j := Cte 2 (ρE) j -E j Cte 2 (ρ) j . (4.38) and :
(Hess rest (ρU)

) j := ρ j (∇ 2 U) j + (∇ρ) j ⊗ (∇U) j + (∇U) j ⊗ (∇ρ) j , (4.39)
(Hess rest (ρE)) j := ρ j (∇ 2 E) j + (∇ρ) j ⊗ (∇E) j + (∇E) j ⊗ (∇ρ) j . (4.40)
In this case, direct limitation is an adaptation of [START_REF] Hoch | A frame invariant and maximum principle enforcing second-order extension for cell-centered ALE schemes based on local convex hull preservation[END_REF].

1. Barth/Jespersen, Dukowicz [START_REF] Barth | The design and application of upwind schemes on unstructured meshes, chapter 6[END_REF] for scalar variables (for volumic or massic variables):

2. Vip [START_REF] Luttwak | Slope limiting for vectors: A novel vector limiting algorithm[END_REF] for vector variables (for volumic or massic variables).

In order to obtain a stability control for internal massic energy, we need to link the limitation coefficient between all the involved variables U and E.

Direct limitation of Density, Indirect Velocity and Total massic Energy

In this part, we look at direct limitation process dealing with density, velocity U and total massic energy E.

P 1,h j (x, ρ) = ρ j + (∇ρ) j .(x -x j ), (4.41)
P 1,h j (x, ρU) = (ρU) j + (∇(ρU)) j (x -x j ), (4.42)
P 1,h j (x, ρE) = (ρE) j + (∇(ρE)) j .(x -x j ) (4.43)
In this case, we adjust the same limitation coefficient for U and E, because both are linked to massic internal energy , we call it α j .

P 1,h j (x, ρ) = ρ j + α ρ j (∇ρ) j .(x -x j ), (4.44)
R 1,h j (x, U) = U j + ρ j P 1,h j (x, ρ) α j (∇U) j (x -x j ), (4.45)
R 1,h j (x, E) = E j + ρ j P 1,h j (x, ρ) α j (∇E) j .(x -x j ) (4.46)
In this case, using

R 1,h j (x, ) := R 1,h j (x, E) -1 2 |R 1,h j (x, U)| 2
, a possible reconstruction of internal energy in terms of velocity U and total massic energy E:

R 1,h j (x, ) = j + ρ j P 1,h j (x, ρ) α j ((∇E) j -( t ∇U) j U j )(x -x j ) - 1 2 ( ρ j P 1,h j (x, ρ) ) 2 (α j ) 2 |(∇U) j (x -x j )| 2 (4.47)
Now at each dof of cell, the positivity of reconstructed internal energy is (generaly) a quadratic inequality, namely:

P 2 (α j ) := A (α j ) 2 + B (α j ) + C > 0 1 . (4.48)
Indeed except the case where A = 0 corresponding to either

-ρ j = 0 or -(∇U) j = 0 or -(x dof -x j ) ∈ Ker(∇U) j .
this is a quadratic equation and so we have A < 0 and C > 0. So that there are always 2 real roots, and the parabola is always positive between the two roots. Let α dof + the positive one (strictly because P 2 (0) = C > 0). In order to obtain a unique coefficient in cell Ω j , we finally take: We begin with the reconstruction of conservative quantities, each polynomial reconstruction P 1,h writes:

α j = min dof min(α dof + , 1) (4
P 1,h j (x, ρ) = ρ j + (∇ρ) j .(x -x j ), (4.50)
P 1,h j (x, ρU) = (ρU) j + (∇(ρU)) j (x -x j ), (4.51)
P 1,h j (x, ρE) = (ρE) j + (∇(ρE)) j .(x -x j ) (4.52)
with the algebraic relation:

E = + 1 2 |U| 2 , (4.53)
using Leibniz formula on both differential terms:

∇(ρU) = ρ∇U + U ⊗ ∇ρ, (4.54)
∇(ρE) = ρ∇E + E∇ρ, (4.55)
∇(ρ ) = ρ∇ + ∇ρ, (4.56)
∇(ρ |U| 2 2 ) = ρ∇ |U| 2 2 + |U| 2 2 ∇ρ (4.57)
using (4.55),(4.56),(4.57), we rewrite (4.51),(4.52):

P 1,h j (x, ρU) = (ρU) j + (ρ j (∇U) j + U j ⊗ (∇ρ) j )(x -x j ), (4.58)
P 1,h j (x, ρE) = (ρE) j + (∇(ρE)) j .(x -x j ) (4.59) = (ρ( + 1 2 |U| 2 )) j + (∇(ρ( + 1 2 |U| 2 )) j .(x -x j ) (4.60) = (ρ ) j + (ρ |U| 2 2 ) j + (∇(ρ )) j + (∇(ρ |U| 2 2 )) j .(x -x j ) (4.61)
We group all terms involving differential of ρ:

P 1,h j (x, ρE) = ( j + |U j | 2 2 )(ρ j + (∇ρ) j .(x -x j )) + ρ j (∇ ) j + (∇ |U| 2 2 ) j .(x -x j ) (4.62) using ∇ |U| 2 2 = t ∇U.U with ∇U = ∂ 1 u 1 ∂ 2 u 1 ∂ 1 u 2 ∂ 2 u 2 , U = (u 1 , u 2 ), (4.63)
we finally obtain that the reconstruction of ρU in terms of ρ, U (and their derivatives):

P 1,h j (x, ρU) = ρ j U j + (ρ∇U + U ⊗ ∇ρ) j (x -x j ) (4.64)
and the reconstruction of ρE in terms of ρ, U, (and their derivatives):

   P 1,h j (x, ρE) = ( j + |Uj | 2 2 ) (ρ j + (∇ρ) j .(x -x j )) + ρ j ((∇ ) j + ( t (∇U) j .U j )) .(x -x j ), ρE = ρ( + 1 2 |U| 2 ). (4.65)
so that for any massic quantity S, wherever P 1,h j (x, ρ) > 0, we define his nonlinear reconstruction by the rational function R 1,h : R 1,h j (x, S) := P 1,h j (x, ρS) P 1,h j (

We apply (4.66) (4.64) at massic variable U and (4.66) (4.65) for massic variable E:

R 1,h j (x, U) := U j + ρ j P 1,h j (x, ρ) (∇U) j (x -x j ), (4.67)
R 1,h j (x, E = + |U| 2 2 ) := j + |U j | 2 2 + ρ j P 1,h j (x, ρ) (∇ ) j + ( t (∇U) j .U j ) .(x -x j ) (4.68)
going further, we impose the same relation of (4.53) in terms of the reconstruction, namely:

R 1,h j (x, E) = R 1,h j (x, ) + 1 2 |R 1,h j (x, U)| 2 . (4.69)
this gives the reconstruction of variable:

R 1,h j (x, ) = R 1,h j (x, E) - 1 2 |R 1,h j (x, U)| 2 . ( 4
R 1,h j (x, ) = j + ρ j P 1,h j (x, ρ) (∇ ) j .(x -x j ) + |U j | 2 2 + ρ j P 1,h j (x, ρ) ( t (∇U) j .U j ) .(x -x j ) (4.71)
- 1 2 |U j | 2 + 2 ρ j P 1,h j (x, ρ) U j ((∇U) j (x -x j ))) + ( ρ j P 1,h j (x, ρ) ) 2 |(∇U) j (x -x j )| 2 (4.72)
Finally, after some simplifications, we deduce that polynomial or rational (nonlinear) reconstructions for physical data are given by:

P 1,h j (x, ρ) = ρ j + (∇ρ) j .(x -x j ), (4.73)
R 1,h j (x, U) = U j + ρ j P 1,h j (x, ρ) (∇U) j (x -x j ), (4.74)
R 1,h j (x, ) = j + ρ j P 1,h j (x, ρ) (∇ ) j .(x -x j ) - 1 2 | ρ j P 1,h j (x, ρ) (∇U) j (x -x j )| 2 (4.75)
These forms are essentials in making a-priori or a-posteriori limitation using direct control on internal massic energy.

Whatever the limitation strategy, we emphasize that gathered high order polynomial reconstructions of conservative variables (ρ, ρU, ρE) are then given by:

P 1,h j (x, ρ) = ρ j + (∇ρ) j .(x -x j ), (4.76)
P 1,h j (x, ρU) = P 1,h j (x, ρ)R 1,h j (x, U), (4.77)
P 1,h j (x, ρE) = P 1,h j (x, ρ) R 1,h j (x, ) + 1 2 |R 1,h j (x, U)| 2 . (4.78)
In the following, we propose different strategies for limitation of reconstructions (4.73),(4.74),(4.75) that we use to evaluate final polynomial reconstructions of conservative variables in (4.76) (4.77) (4.78).

4.2.3.b Third order

We begin with reconstruction of conservative quantities, each polynomial reconstruction P 2,h writes:

P 2,h j (x, ρ) = ρ j + (∇ρ) j .(x -x j ) + 1 2 t (x -x j )(∇ 2 ρ) j (x -x j ), (4.79)
P 2,h j (x, ρU) = (ρU) j + (∇(ρU)) j (x -x j ) + 1 2 t (x -x j )(∇ 2 (ρU)) j (x -x j ), (4.80)
P 2,h j (x, ρE) = (ρE) j + (∇(ρE)) j .(x -x j ) + 1 2 t (x -x j )(∇ 2 (ρE)) j (x -x j ). (4.81)
where we have set

Q j := Q j , if conservation property of reconstruction is not needed, Q j -1 2 1 |Ωj | Ωj t (x -x j )(∇ 2 Q) j (x -x j )dv := Q j -Cte(Q) j , else. (4.82)
for any conservative quantity Q (=ρ, ρU and ρE). Using (4.54),(4.55),(4.57), we deduce that Leibniz formula writes:

∇ 2 (ρU) = ρ∇ 2 u 1 u 2 + (∇ 2 ρ) u 1 u 2 + ∇ρ ⊗ ∇ u 1 u 2 + ∇ u 1 u 2 ⊗ ∇ρ, (4.83)
∇ 2 (ρE) = ρ∇ 2 E + E∇ 2 ρ + ∇ρ ⊗ ∇E + ∇E ⊗ ∇ρ, (4.84) ∇ 2 (ρ ) = ρ∇ 2 + ∇ 2 ρ + ∇ρ ⊗ ∇ + ∇ ⊗ ∇ρ, (4.85)
∇ 2 (ρ |U| 2 2 ) = ρ∇ 2 |U| 2 2 + |U| 2 2 ∇ 2 ρ + ∇ρ ⊗ ∇ |U| 2 2 + ∇ |U| 2 2 ⊗ ∇ρ, with (4.86)
∇ |U| 2 2 = t ∇UU and ∇ 2 |U| 2 2 = t ∇U∇U + i u i ∇ 2 u i (4.87) using (4.83), (4.84),(4.85),(4.86), we rewrite (4.80),(4.81)
:

P 2,h j (x, ρU) = (ρU) j + (ρ j (∇U) j + U j ⊗ (∇ρ) j ) (x -x j ) + (4.88) 1 2 t (x -x j ) ρ j (∇ 2 u 1 u 2 ) j + (∇ 2 ρ) j u 1 u 2 j + (∇ρ) j ⊗ (∇ u 1 u 2 ) j + (∇ u 1 u 2 ) j ⊗ (∇ρ) j (x -x j ) (4.89)
using (4.82):

(ρU) j = (ρU) j -Cte(ρU) j (4.90

(ρU) j = (ρU) j - 1 2 1 |Ω j | Ωj t (x -x j ) ρ j (∇ 2 u 1 u 2 ) j + (∇ 2 ρ) j u 1 u 2 j + (∇ρ) j ⊗ (∇ u 1 u 2 ) j + (∇ u 1 u 2 ) j ⊗ (∇ρ) j (x -x j )d = (ρU) j - 1 2 1 |Ω j | Ωj t (x -x j ) (∇ 2 ρ) j u 1 u 2 j (x -x j )dv - 1 2 1 |Ω j | Ωj t (x -x j )H rest,ρu j (x -x j )dv
where we have noted:

H rest,ρu j = ρ j (∇ 2 u 1 u 2 ) j + (∇ρ) j ⊗ (∇ u 1 u 2 ) j + (∇ u 1 u 2 ) j ⊗ (∇ρ) j (4.91) Noting: CteRest(ρU) j = 1 2 1 |Ωj | Ωj t (x -x j )H rest,ρu j (x -x j
)dv, we also have:

-CteRest(ρU) j = (ρU) j -U j ρ j , or equivalently (4.92) CteRest(ρU) j = Cte(ρU) j -U j Cte(ρ) j .

(4.93)

Bringing together derivative of ρ, the polynomial reconstruction of ρU:

P 2,h j (x, ρU) = U j (ρ) j + (∇ρ) j (x -x j ) + 1 2 t (x -x j )(∇ 2 ρ) j (x -x j ) -CteRest(ρU) j (4.94) +ρ j (∇U) j (x -x j ) + 1 2 t (x -x j ) ρ j (∇ 2 u 1 u 2 ) j + (∇ρ) j ⊗ (∇ u 1 u 2 ) j + (∇ u 1 u 2 ) j ⊗ (∇ρ) j (x -x j )(4.95)
Now, for the total energy ρE:

P 2,h j (x, ρE) = (ρE) j + (∇(ρE)) j (x -x j ) + 1 2 t (x -x j )(∇ 2 (ρE)) j (x -x j ) (4.96) P 2,h j (x, ρE) = (ρ( + 1 2 |U| 2 )) j + (∇(ρ( + 1 2 |U| 2 ))) j (x -x j ) + 1 2 t (x -x j )(∇ 2 (ρ( + |U| 2 2 ))) j (x -x j ) (4.97) = (ρ ) j + (ρ |U| 2 2 ) j + (∇(ρ )) j + (∇(ρ |U| 2 2 )) j (x -x j ) + 1 2 t (x -x j )((∇ 2 (ρ )) j + (∇ 2 (ρ |U| 2 2 )) j )(x -x j ) with (ρ ) j = (ρ ) j - 1 2 1 |Ω j | Ωj t (x -x j )(ρ∇ 2 + ∇ 2 ρ + ∇ρ ⊗ ∇ + ∇ ⊗ ∇ρ) j (x -x j )dv (4.98) = (ρ ) j - 1 2 1 |Ω j | Ωj t (x -x j )( ∇ 2 ρ) j (x -x j )dv - 1 2 1 |Ω j | Ωj t (x -x j )H rest,ρ j (x -x j )dv (4.99) (ρ |U| 2 2 ) j = (ρ |U| 2 2 ) j - 1 2 1 |Ω j | Ωj t (x -x j )(ρ∇ 2 |U| 2 2 + |U| 2 2 ∇ 2 ρ + ∇ρ ⊗ ∇ |U| 2 2 + ∇ |U| 2 2 ⊗ ∇ρ) j (x -x j )dv (4.100) = (ρ |U| 2 2 ) j - 1 2 1 |Ω j | Ωj t (x -x j )( |U| 2 2 ∇ 2 ρ) j (x -x j )dv - 1 2 1 |Ω j | Ωj t (x -x j )H rest,ρ |U| 2 2 j (x -x j )dv (4.101)
where in (4.99) and (4.101)

H rest,ρ j := (ρ∇ 2 + ∇ρ ⊗ ∇ + ∇ ⊗ ∇ρ) j , (4.102) H rest,ρ |U| 2 2 j
:= (ρ∇ 2 |U| 2 2 + ∇ρ ⊗ ∇ |U| 2 2 + ∇ |U| 2 2 ⊗ ∇ρ) j (4.103)
so that summing (4.98) (4.100):

(ρ ) j + (ρ |U| 2 2 ) j = (ρE) j = E j ρ j -CteRest(ρE) j (4.104)
We group all terms involving differential of ρ:

P 2,h j (x, ρE) = ( j + |U j | 2 2 ) ρ j + (∇ρ) j (x -x j ) + 1 2 t (x -x j )(∇ 2 ρ) j (x -x j ) (4.105) -CteRest(ρE) j + ρ j (∇ ) j + (∇ |U| 2 2 ) j (x -x j ) (4.106) + 1 2 t (x -x j ) ρ j ((∇ 2 ) j + (∇ 2 |U| 2 2 ) j ) + (∇ρ) j ⊗ ((∇ ) j + (∇ |U| 2 2 ) j ) + ((∇ ) j + (∇ |U| 2 2
)) j ⊗ (∇ρ) j (xx j ) (4.107) using (4.87), we have:

P 2,h j (x, ρE) = ( j + |U j | 2 2 )P 2,h j (x, ρ) -CteRest(ρE) j + ρ j (∇ ) j + ( t ∇U) j U j (x -x j ) (4.108) + 1 2 t (x -x j )H rest,ρE j (x -x j ) (4.109)
where

H rest,ρE j = ρ j ((∇ 2 ) j + ( t ∇U) j (∇U) j + i u i j (∇ 2 u i ) j ) + (∇ρ) j ⊗ ((∇ ) j + ( t ∇U) j U j ) + ((∇ ) j + ( t ∇U) j U j ) ⊗ (∇ρ) j
and we have also:

CteRest(ρE) j = 1 2 1 |Ωj | Ωj t (x -x j )H rest,ρE j (x -x j )dv.
Now, using the algebraic relation between total energy, internal and kinetic energy in their counterparts reconstruction:

R 2,h j (x, ) = R 2,h j (x, E) - 1 2 |R 2,h j (x, U)| 2 . (4.110)
To obtain a reconstruction for internal massic energy, we first use the relation (4.94)

R 2,h j (x, U) := P 2,h j (x, ρU) P 2,h j (x, ρ) (4.111)
and the relation (4.108)

R 2,h j (x, E) := P 2,h j (x, ρE) P 2,h j (x, ρ) (4.112)
The reconstruction of physical variables are finally given by:

P 2,h j (x, ρ) = ρ j + (∇ρ) j (x -x j ) + 1 2 t (x -x j )(∇ 2 ρ) j (x -x j ), (4.113)
R 2,h j (x, U) = U j + -CteRest(ρU) j + ρ j (∇U) j (x -x j ) + 1 2 t (x -x j )(Hess rest (ρU)) j (x -x j ) P 2,h (x, ρ) (4.114) R 2,h j (x, ) = j + -CteRest(ρ ) j + ρ j (∇ ) j (x -x j ) + 1 2 t (x -x j )(Hess rest (ρ )) j (x -x j ) P 2,h (x, ρ) (4.115)
+ -M M T 2(ρ( t ∇U)(∇U)) j + 1 2 t (x -x j )(ρ j ( t ∇U) j (∇U) j )(x -x j ) P 2,h (x, ρ) - 1 2 -CteRest(ρU) j + ρ j (∇U) j (x -x j ) + 1 2 t (x -x j )(Hess rest (ρU)) j (x -x j )
P 2,h (x, ρ) 2 using (4.82):

(Hess rest (ρU)) j := ρ j (∇ 2 U) j + (∇ρ) j ⊗ (∇U) j + (∇U) j ⊗ (∇ρ) j -CteRest(ρU) j := (ρU) j -U j ρ j = -1 2

1 |Ω j | Ωj t (x -x j )(Hess rest (ρU)) j (x -x j )dv, M M T 2(ρ( t ∇U)(∇U)) j := 1 2 1 |Ω j | Ωj t (x -x j )(ρ j ( t ∇U) j (∇U) j )(x -x j )dv (4.116) (Hess rest (ρ )) j := ρ j (∇ 2 ) j + (∇ρ) j ⊗ (∇ ) j + (∇ ) j ⊗ (∇ρ) j , (4.117)
-CteRest(ρ

) j := (ρ ) j -j ρ j = - 1 2 1 |Ω j | Ωj t (x -x j )(Hess rest (ρ )) j (x -x j )dv,
where

(∇ 2 U) j = (∇ 2 u 1 u 2 ) j ,
(∇u i is written by rows).

here also as for second order reconstruction, whatever the limitation strategy used, the final polynomial reconstruction of conservative variable, for arbitrary n order:

P n,h j (x, ρ) = ρ j + (∇ρ) j (x -x j ) + 1 2 t (x -x j )(∇ 2 ρ) j (x -x j ) + .., (4.118)
P n,h j (x, ρU) = P n,h j (x, ρ)R n,h j (x, U), (4.119)

P n,h j (x, ρE) = P n,h j (x, ρ) R n,h j (x, ) + 1 2 |R n,h j (x, U)| 2 .
(4.120) -a hierarchical strategy λ i ρ for i th order term ∇ i ρ:

P n,lim j (x, ρ) = ρlim_hi j + λ 1 ρ (∇ρ) j (x -x j ) + 1 2 t (x -x j )(λ 2 ρ )(∇ 2 ρ) j (x -x j ) + .. ( 4
(Q) := 1 k! 1 |Ω j | Ωj λ k Q (∇ k Q) j (x -x j , .., x -x j )dv. (4.127)
-the same function λ ρ for all high order terms (in the same way than (4.123) (4.124)):

P n,lim j (x, ρ) = ρlim_all j + λ ρ (∇ρ) j (x -x j ) + 1 2 t (x -x j )(∇ 2 ρ) j (x -x j ) + .. ( 4
ρlim_all j = ρ j -λ ρ n k=1 Cte k j (ρ) (4.129)
Here, λ ρ in (4.128) or each λ i ρ in (4.125) are designed to obtain at least the positivity criterium P n j (x, ρ) > 0 for all x dof points of Ω j (but also in practice min

V ρ V (j) ≤ P n j (x, ρ) ≤ max V ρ V (j) ).
Remark 4.2. In (4.123), (4.124), we have the following expression for P n polynomials:

-For n=1:

P 1 ≡ 0 -For n=2: P 2 = -M M T 2(ρ( t ∇U)(∇U)) j + 1 2 t (x -x j )(ρ j ( t ∇U) j (∇U) j )(x -x j )
where MMT2 is defined in (4.116).

Choices on limitor λ ρ , λ

In view of (4.8), we need to make sense of flux evaluated with high order reconstructed value at dof of each cells. Many limitor for density (4.128) and for internal energy (4.121) (that induces the one of velocity (4.122)) can be investigated, here we use a minmod/Dukowicz like limitor constructed in the spirit of (4.128), for each cell we define a scalar cell limitor λ ρ j :

-P n,lim j (x, ρ), the density limitation with strict high order part:

P n,lim j (x, ρ) = ρ j + λ ρ j (P n j (x, ρ) -ρ j ). (4.130)
Now, we explain the construction of λ ρ j :

1. Evaluation of unlimited P n j (x, ρ) at x dof in boundary cell Ω j : the associated values P n j (x dof , ρ) are noted P n j,dof , we also note:    M p j := max dof P n j,dof , mp j := min dof P n j,dof .

(4.131) 2. Evaluation of unlimited P n j (x, ρ) at x j centroïd of neighboring cell Ω j , j ∈ V j (:=V n (j) here see notation (3.97)): the associated values P n j (x j , ρ) are noted P n j,j , we also note:    M pnc j := max j ∈Vj P n j,j , mpnc j := min j ∈Vj P n j,j .

(4.132)

Computation of extremal values of previous reconstructed values:

M R j := max(M p j , M pnc j ), mR j := min(mp j , mpnc j ).

           s ρ,max j := 1, if |M R j -ρ j | << 1, M Nj -ρj M Rj -ρj , else. s ρ,min j := 1, if |mR j -ρ j | << 1, mNj -ρj mRj -ρj .
else.

(4.135)

Finally, we take : R n,lim j (x, U) = U j + λ U j Rj (x, U), and using (λ U j = λ j ) (where λ j is defined (4.136)), we get then in (4.124):

λ ρ j = max(0,
R n,lim j (x, U) = U j + λ j Rj (x, U).

Remark 4.3.

We emphasize here that other a priori limitation strategies are possible, and will be investigated in other studies.

In case of the weaker positivity condition on a target quantity f (= ρ (density) or (massic internal energy) ), we just put s f,max j = 1 and mN j = 0 in (4.135).

Time order

In our case, to obtain high order time integration, we use RK3-TVD of [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF] for both second or third order spatial reconstruction (4.118), (4.119), (4.120).

RK3 -TVD

           f (1) h = f n h + ∆tRHS(f n h ), f (2) h = f n h + ∆tRHS(f (1) h ), f (3) h = ( 3f n h +f (2) h 4 ) + ∆tRHS( 3f n h +f (2) h 4
),

f n+1 h = f n h +2f (3) h 3 (4.137)
It is a convex combination of first order Euler time steping:

f update h = f ref h + ∆tRHS(f ref h ) (4.138)
More precisely, for next numerical simulation section, we will use (4.137) for both second and third order space reconstruction. For the first order space reconstruction, one use Euler time steping (4.138). We apply the previous a posteriori limitation to composite nodal/edge numerical schemes for the Euler equations with

Adm = {U := (ρ, ρU , ρE); ρ ≥ 0; = E - 1 2 |U | 2 ≥ 0}. (4.148)
In case of massic quantities , we recall that the approach is slightly different, and we need to consider the reconstruction of non conservative variable U and . In view of (4.124) and (4.123), which are not polynomial, we can not use directly degree limitation (4.143). But if we consider (4.124) with a constant density field, we obtain a polynomial:

R n,lim j (x, U) ρ≡Cte = P n,lim (x, U) = α U,0 j + α U,1 j (∇U) j (x -x j ) + 1 2 (x -x j ) t α U,2 j (∇ 2 U) j (x -x j ) +... + 1 n! (α U,n j )(∇ n U) j (x -x j , .., x -x j ) (4.149)
using same procedure in [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF], we deduce a polynomial reconstruction for momentum:

= α ,0 j + α ,1 j (∇ ) j (x -x j ) + 1 2 (x -x j ) t α ,2 j (∇ 2 ) j (x -x j ) +... + 1 n! (α ,n j )(∇ n ) j (x -x j , .., x -x j ) (4.151)
and then deduce a limited polynomial reconstruction for total energy:

R n,lim j (x, ρE) = R n,lim j (x, ρ( + 1 2 |U| 2 )) = P n,lim (x, ρ( + 1 2 |U| 2 ), (D(d ρ j ), D(d U j ), D(d j ))) (4.152)
Remark 4.6. In view of a priori limitation (4.130), (4.123), (4.136). It is possible to obtain a more compact limitation for both velocity U and internal massic energy . Indeed, a more natural way is to apply a unique coefficient

(α ,l j ) n l=1 ≡ (α j ), ∀l, (4.153)
so that, now we have a non hierarchical limitation (we then deduce that of U by the constraint (α U j = α j )). In this case a more natural reconstructed polynomial is: R n,lim j (x, ρU) = P n,lim (x, (ρU), (D(d ρ j ), α j )) (4.154)

R n,lim j (x, ρE) = R n,lim j (x, ρ( + 1 2 |U| 2 )) = P n,lim (x, ρ( + 1 2 |U| 2 ), (D(d ρ j ), α j , α j )) (4.155)
so that only ONE limitor dictates ALL massic quantities U and .

In the following, we are able to give two version of APITALI algorithms (a hierarchical and another more compact) to obtain a valid numerical solution U ∈ Adm(4.148) if the underlying first order reconstruction scheme verifies this stability.

d i C * ,U C * → d i C * +1,U C * d i C * , C * → d i C * +1, C * (4.162) Eventually      For C ∈ V (C * )( / ∈ BadCellSet) : d i C ,U C → d i C +1,U C d i C , C → d i C +1, C
) ii. Else Bad 0 = ∅ then M T = ∅ : For C * ∈ BadCellSet : α i C * , C * → α i C * +1, C * (4.171) Eventually For C ∈ V (C * )( / ∈ BadCellSet) : α i C , C → α i C +1, C (4.172) M T ← {C ∈ V(C) with C ∈ M T such that C ∈ (4.

Numerical results of high order composite nodal/edge schemes

We present some test cases with shocks and non smooth solutions (Sod, mach3 step and double mach reflection), we want to asses the robustness of third order non-linear reconstructions either on degenerate conical or on truly curved meshes. We begin with one-dimensional test cases and continue with truly two-dimensional ones.

Smooth solutions

In this part, we will focus on obtaining the third order (numericaly) of the schemes (thus on reconstruction, quadrature points and time order) without any spatial limitation process.

Taylor-Green Vortex

We compute the stationnary Taylor-Green vortex on domain :

[-5, 5] × [-5, 5].
The solution ρ(t, x) = (T ∞ + δT ) 1/(γ-1) , U(t, x) = U ∞ + δU, p(t, x) = ρ γ , with x = (x 1 , x 2 ) and:

δU(t, x) = β 2π exp( 1-r 2 2 )(-x 2 , x 1 ) δT (t, x) = -(γ-1)β 2 8γπ 2 exp(1 -r 2 ) ( 5

Regular Periodic Solution

The example is taken in [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF] to check schemes accuracy.

ρ(t, x) = 1 + 0.2sin(π(x 1 + x 2 -t(u 1 + u 2 ))) (5.2)
where u 1 = 1, u 2 = -1/2 (the components of velocity), and P = 1 (the pressure), γ = 1.

Woodward-Colella problem: Blast Waves

This benchmark problem involves strong shocks that reflect on the extremities of the tube and interact further. For this test, final time is t f = 0.038, and the initial data is given by: (ρ 0 , u 0 , p 0 ) =    (1, 0, 1000), 0 ≤ x < 0.1 (1, 0, 0.01), 0.1 < x < 0.9 (1, 0, 100), 0.9 < x ≤ 1.

The Osher-Shu problem

This test case is used to demonstrate the capability of high order schemes to capture small structure in smooth regions of Euler flows. More precisely, a moving mach3 shock interacts with a smooth periodic perturbation in density, final time is t f = 1.8, and the initial data is given by: (ρ 0 , u 0 , p 0 ) = (3.857143, 2.629369, 10.333333), -5 ≤ x < -4, (1 + 0.2 sin(5x), 0, 1.), -4 < x ≤ 5. This two dimensional test case, is an extension of 1D shock tube. In this case, the unit square is divided into four quadrants by lines x 1 = 1/2, x 2 = 1/2. The 2D Riemann problem is defined by initial constant states in each quadrant which are given :

(ρ 0 , u 1,0 , u 2,0 , p 0 ) = (0.5197, -0.6259, -0.3, 0. We want now to solve a radial sod shock tube, the polar circular mesh (see Figure 39) is made of radial circular edges while orthoradial are made with degenerate conical edges. The nodal scheme of Rusanov (Version 2) is used here, and third order limited reconstruction.

(ρ 0 , u 1,0 , u 2,0 , p 0 ) = (1, 0, 0, 1), 0 ≤ r ≤ 0.5 (0.125, 0, 0, 0.1), 0.5 ≤ r ≤ 1. A well known benchmark is the Double Mach reflection problem described in [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. This test case is a extension of a two dimensionnal sod shock tube. On a rectangular domain [0, 2]×[0, 1] containing two curved circular holes, from each part of {x 1 = 1} the two following states are imposed, and symmetry is imposed on exterior rectangle (North,East,West,South), U.N = 0 is imposed at each node/edge of circular holes.

(ρ 0 , u 1,0 , u 2,0 , p 0 ) = (1, 0, 0, 1), 0 ≤ x 1 ≤ 1 (0.125, 0, 0, 0.1), 1 ≤ x 1 ≤ 2. This test case presented in [START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF] shows the loss of viscosity due to higher order (three or more) reconstruction fails to predict a good numerical solution. Indeed the authors demonstrates that only a second order minmod limiter permits to capture the entropic solution, but never the third order extension. As shown in our experiments, we have the same conclusion see Figure 50.

Conclusion and prospects

We have proposed an extension of some well known Approximate Riemann Solvers (ARS) defined on edge to the nodes of arbitrary straight or conical meshes. This defines composite nodal/edge flux schemes that we call CARS (Composite Approximated Riemann Solver) for which conical schemes are a special case.

We have formulated algebraic relations to obtain consistent and localy conservative schemes. These are based on making a difference between two families: flux type and viscosity type schemes, keeping them separated makes an easier understanding and a better implementation. Moreover, the design of conical composite flux schemes permits to reach more than second (third or fourth) order contrary to pure polygonal scheme involving only normal at nodes (as in seminal work dedicated to Lagrange update frame [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] for some extensions see [START_REF] Vilar | A discontinuous Galerkin discretization for solving the twodimensional gas dynamics equations written under total lagrangian formulation on general unstructured grids[END_REF][START_REF] Cheng | A third order conservative Lagrangian type scheme on curvilinear meshes for the compressible Euler equations[END_REF]).

This new class of composite finite volume schemes make a continuous variation between classical pure edge ARS scheme to pure nodal ARS scheme.

This formalism share some of the hypothesis in the design of Distributive Residual schemes, Proof. Let κ 0 (q) :=

α ω 1/2 =                  2 √ 1-ω 2 atan( 1-ω 1+ω ) (1-ω) 2 + ω+1 ω-1 0 ≤ ω < 1,
(1-q)(1-2q)

(1-q) 2 +2ωq(1-q)+q 2 , κ 1/2 (q) := 2(1+ω)q(1-q)

(1-q) 2 +2ωq(1-q)+q 2 , κ 1 (q) := q(2q-1)

(1-q) 2 +2ωq(1-q)+q 2 . we juste compute α ω i := 1 0 κ i (q)dq, using single element decomposition of these functions. -The degenerate conical case corresponding to ω = 0 gives: α 0 = α 1 = 1 -π 4 ( .2146018) and α 1/2 = π 2 -1( .5707963). The limit ω → ∞ gives α 0 = α 1 = 0 and α 1/2 = 1.

-For ω = 1, it corresponds to a parabolic Lagrange polynomial interpolation, giving the Simpson rule.

-It exists a link with quadrature formula α i and h(ω) in (2. The two others α ω 0 , α ω 1 are deduced by (B.7).

Figure 1 :

 1 Figure 1: Mesh dependencies for first order finite volume schemes with edge based numerical flux Φ (1.8): only cells sharing a common edge (green cells) are involved , corner cells (red ones) do not appear. 3

Figure 4 :

 4 Figure 4: Numerical travel time: first order explicit finite volume pure edge scheme depends highly on cell/edge repartition. Consider a fluid initialy at rest and one cell (center (0.575, 0.575)) density field, the numerical arrival time in the cell at left/bottom corner (center (0.475, 0.475)) highly depends on the NUMBER of edges which separates them.

Figure 5 :

 5 Figure 5: Composite Nodal/Edge fluxes for arbitrary conical meshes (with straight or curved edges), edge fluxes involve only green neighbor cells (see Figure 1) while nodal fluxes involve both green and red ones. Straight polygonal meshes must be a special case of conical meshes for ω = 0 (in this case the mid-edge flux still contribute).

Figure 6 :

 6 Figure 6: Conic parametrized by quadratic rational Bezier curves.

Figure 8 :

 8 Figure 8: Normal nodal vector expressed with control point

Figure 9 :

 9 Figure 9: Normal vector from a degree of freedom defined on boundary cell Ω j . Two types: endpoints M r denoted by Cr,ω j or shoulder point S ω r+1/2 denoted by Cr+1/2,ω j

Figure 10 :

 10 Figure 10: Vector decomposition at dof (nodes/edge) (of GLACE type vector) Cdof,ω j into two half-vectors to left Ñdof,ω+ j

Figure 11 :

 11 Figure 11: Flux G dof j located at each dof of a cell Ω j : a quadrature formula point of view

  at shoulder point are the same. Hence by identification of (3.20) and (3.21), we enforce a functional constraint on numerical nodal schemes:

Proposition 3 . 9 .

 39 With the assumption (TLBL) (3.36), the nodal numerical fluxes (3.33) and (3.34) are consistent and locally conservative around each dof verifying (2.21) (interior dof ).

Figure 13 :

 13 Figure 13: Opposite normal corners (at dof) vectors need to be taken into account for local conservation of nodal extension for VFFC or Roe flux schemes

  Roe edge scheme (in flux form) to nodal flux (Version 2) Definition 3.16. A possible extension of classical Roe edge scheme (3.30) to nodes (when it exist at the edge) is to consider the following -the nodal flux (3.34) (named Version 2).

( 3 .

 3 63) then tensors Λ up dof,k in (3.63) satisfy TLBL (3.36).

  .70) Proof. With the Heaviside H function (3.61), the schemes (3.70) rewrites :

( 3 .

 3 72) because of the fact k β k = 1, choosing B k = (β k )I in (3.63) gives the result.

  (3.78) by omitting the dependency of S min / max jk on ( Njk ) in (3.78):

  .83) Proposition 3.22. The numerical fluxes (3.81) and (3.82) are consistent and locally conservative around each dof verifying (2.21) (interior dof ).

3. 3 . 2

 32 Extension of Roe scheme (in viscous form) to nodal flux Definition 3.25. A possible extension of classical Roe edge scheme (3.76) to nodes (when it exist at the edge) is to consider the following -the nodal flux (3.82) (named Version 2).

  VFRoe scheme to nodal flux Definition 3.26. A possible extension of classical VFRoe edge scheme (3.77) to nodes (when it exists at the edge) is to consider the following -the nodal flux (3.81) (Version 1) or (3.82) (Version 2) and -The viscosity matrix is an average of absolute value of the Jacobian matrix (evaluated with mean state values see (3.26) and (3.77) for the notations) :

Figure 14 :

 14 Figure 14: Composite nodal/edges schemes make the set V e j multiple exchanging cell because of V e j ⊂ V n j . Red hashed cells are those belonging to nodal flux at both extremities of edge (Definition 3.97), green cells are those belonging to (shoulder) edge flux (Definition 3.96).

Figure 15 :Figure 16 :

 1516 Figure 15: From any cell Ω j , any edge neighbor cells in V e j (3.96) is counted three times (one from edge itself and the others from his two end points), any of corner neighbor cells in V c j (3.98) is counted once.

.111) Proposition 3 . 31 .

 331 With this following (sufficient) conditions:

j ≥ 0 ,

 0 (3.116) we obtain that if the numerical flux G dof j • Cdof j is increasing w.r.t state in j, a condition on the time step: gives (3.113).

Figure 17 :

 17 Figure 17: Left: Generic conical unstructured cells with all his subcells attached to his own dofs . Right: Dual cells control volume around dof (nodes M r in green and mid edges points S ω r+1/2 in blue): partial dual diamond cell (blue) for mid edge points and partial dual nodal cell (green) for nodes.

  .125) we recover a convex combination. So that in view of (3.123), the aim is to obtain that each Ū dof j in (3.125) must belong to C, which means that under properties on numerical nodal fluxes G dof j , a time step restriction on ∆t, and the existence of a strictly positive σ. An identical definition of invariant domain preservation is transposed to nodal flux at dof : Definition 3.34. Invariant domain preserving nodal flux at dof Let an invariant domain C ⊆ A adm , the nodal flux G dof j is invariant domain preserving if it exist σ dof j > 0 in (3.123) such that:

Figure 18 : 2 .

 182 Figure 18: A generic cell Ω j and the sub-cell Ω dof j attached to each of his composite dof set.

  due to the CFL (3.133) on σ dof j in (3.126), we have:

Figure 19 :

 19 Figure 19: A one parameter family of dual cells depending on the weight value of θ in composite formula (2.43) (2.44) : Left θ = 0 (Nodal Dual Cell), Middle θ = 1 (Diamond Dual Cell), Right 0 < θ < 1 (Composite Dual Diamond/Nodal Cells).

  the time step restriction (3.154), the new value unknown U n+1 j in (3.153) is indeed a convex combination of data.

Lemma 3 . 40 .

 340 .157) Now, taking ∆x = 1 (ie a unit volume cell), and choosing δt = | Cdof j | 2µ dof jk , we can identify Ū jk in (3.152) as the averaged solution of Riemann problem with data U j , U k in the direction Ñdof j In view of (3.121) (3.122), (3.124) (3.125), and CFL (3.154), we can choose:

Definition 3 . 41 .Definition 3 . 42 .

 341342 As for edge flux VFFC scheme, we define an entropic correction of nodal flux scheme called nodal VFFC-C by replacing the VFFC nodal flux (see Definition 3.17) with those of nodal Rusanov flux (see Definition 3.23) when an eigenvalue is changing sign for a genuinely non-linear field at any dof (this definition could be used either for Version 1 or Version 2). As for edge flux Roe scheme, we define an entropic correction of nodal flux scheme called nodal Roe-C by replacing the Roe nodal flux (see Definition 3.16) with those of nodal Rusanov flux (see Definition 3.23) when an eigenvalues is changing sign for a genuinely non-linear (GNL) field at any dof (this definition is valid for Version 2).

Figure 20 :

 20 Figure 20: Case of one dimensional planar test problem approximated with two dimensional degenerate conical schemes. Symmetric boundary conditions are imposed to up and down boundaries and those of one dimensional physical problem are assigned to left ({x 1 = 0}) and right ({x 1 = L}).

Figure 22 :

 22 Figure 22: Sod shock tube with sonic point: zoom of Figure 21 around sonic point location, Right: continuous correction of VFFC-CC :

  First called VFFC-C is nothing but VFFC except when it exists a CSAGNL event (3.160), where the numerical flux is replaced by nodal Rusanov (3.41). However, this correction is non regular, and we observe some inflexion point in the rarefaction wave at the sonic point location. Using a continuous (Lipschitz) correction called VFFC-CC (3.164), we obtain a better qualitative behavior.

3. 5 . 3

 53 Classical Sod shock tube at final time tf=0.16 Now, as a single test case, we compare the results of previous flux based nodal VFFC-C (Version 2) (3.162) and viscosity based nodal Rusanov (Version 2) of the Sod shock tube at final time t f = 0.16.

.49) 4 . 2 . 3

 423 Direct limitation of Density and Internal massic Energy, Indirect for Velocity 4.2.3.a Second order

  .70) using both equation in (4.67)(4.68), (4.70) writes:

4 .

 4 Computation of extremal bounds of solution in first neighborood:    mN j := min(ρ j , min j ∈Vj ρ j ), M N j := max(ρ j , max j ∈Vj ρ j ). (4.134) 5. Evaluation of min and max gap to obtain spatial local bound preservation:

Proposition 4 . 5 .

 45 For any multi-index sequence of coefficients of limitation (α l,k j ) l,k , we can associate a unique sequence of degrees of limitation (d i j ) 0≤i≤N d where This real valued sequence is strictly decreasing, with d 0 j = n and d N d j = 0, and it includes all the integer between n and 0.

  ρU) = P n,lim (x, (ρU), (D(d ρ j ), D(d U j ))) (4.150) In the same way for internal massic energy (4.123) when considering now a constant density field AND a constant velocity field: R n,lim j (x, ) ρ≡Cte U≡Cte = P n,lim (x, )

Algorithm 1 3 .

 13 APITALI method over a full time step: an uncorrelated hierarchical approach In order to pass from U ref C to U update C (cf (4.138), ∀C ∈ M) for an arbitrary degree's limitation n 1. Init: computation of (∇ l f ) C , l = 1, .., n ⇒ for f = (ρ, U, ) unlimited reconstruction (linear and non-linear) and with or without use of a priori limitation (4.123) and (4.124) (4.122) (4.136). 2. Data: ρ update C computed with volumic APITALI [7] from ρ ref C and then D(d dof ρ ) with R C (x, ρ, D(d dof ρ )): M T = M (testing set where U update C in Adm (4.148) need to checked) (4.156) While M T = ∅: Computation of degrees limitation d dof at quadrature points for U and : (a) ∀dof ∈ M T : applied to R C (x, U, D(d dof U )) ρ=cte one use D(d dof U ), and R C (x, , D(d dof )) ρ=cte U=cte one use D(d dof ). ∀C/dof ∈ C: Flux evaluation at (x = M dof ) with gathered high order terms in reconstruction (4.150) (4.152): R C (x, (ρU), (D(d dof ρ ), D(d dof U ))), R C (x, ρ( + 1 2 |U| 2 ), (D(d dof ρ ), D(d dof U ), D(d dof ))). (4.159) (b) For C ∈ M T , update with scheme : predicted conservative values qupdate C , q = (ρU, ρE). (4.160) (c) Let BadCellSet the set of cells that do not fulfill (4.148), we apply APITALI (4.146) decreasing of degree limitation d U C * and d C * a in cells C * belonging to BadCellSet (and eventually V (C * )). (d) Computation Bad 0 ={C ∈ BadCellSet; d N d ,Bad 0 = ∅ then we re-apply APITALI algorithm [7] on density ρ with a zero degree d ρ C = 0 (locally first order), with M T = Bad 0 then (4.160) to obtain a new ρ update C then go back to (4.158) with M T = M T + Bad 0 + N eighbors and update of U update C and update C (so (ρU) update C and (ρE) update C are locally first order like ρ update C ) ii. Else Bad 0 = ∅ then M T = ∅ :    For C * ∈ BadCellSet :

( 4 ..Algorithm 2 3 .

 423 163) M T ← {C ∈ V(C) with C ∈ M T such that C ∈ (4.162)(and(4.163))} (4.164) iii. Loop to 3. until every cells do verify (4.148) and then U update C APITALI method over a full time step: an correlated compact approach In order to pass from U ref C to U update C (cf (4.138), ∀C ∈ M) 1. Init: computation of (∇ l f ) C , l = 1, .., n ⇒ for f = (ρ, U, ) unlimited reconstruction (linear and non-linear) and with or without use of a priori limitation (4.123) and (4.124) (4.122) (4.136). 2. Data: ρ update C computed with volumic APITALI [7] (or compact α dof ρ coefficient) from ρ ref C : M T = M (testing set where U update C in Adm (4.148) need to checked) (4.165) While M T = ∅: Computation of coefficient limitation α dof at quadrature points for (and induced for U): (a) ∀dof ∈ M T :α dof = min C;dof ∈C α i C , C so that α dof U = α dof . (4.167) ∀C/dof ∈ C: Flux evaluation at (x = M dof ) with gathered high order terms in reconstruction (4.150) (4.152): R C (x, (ρU), (D(d dof ρ), α dof )), R C (x, ρ( + 1 2 |U| 2 ), (D(d dof ρ ), α dof , α dof )). (4.168) (b) For C ∈ M T , update with scheme : predicted conservative values qupdate C , q = (ρU, ρE). (4.169) (c) Let BadCellSet the set of cells that do not fulfill (4.148), we apply APITALI (4.146) limitation α C * in cells C * belonging to BadCellSet (and eventually V (C * )). (d) Computation Bad 0 ={C ∈ BadCellSet; α N d , C = 0}. (4.170) i. If Bad 0 = ∅ then we re-apply APITALI algorithm on density ρ with a local first order ( d ρ C = 0 or α ρ C = 0) et M T = Bad 0 then (4.160) to obtain a new ρ update C then go back (4.158) with M T + = Bad 0 + N eighbors and update U update C and update C (so (ρU) update C and (ρE) update C are locally first order like ρ update C

  162)(and(4.163))} (4.173) iii. Loop to 2. until every cells do verify (4.148) and then U update C .

. 1 )

 1 and ρ ∞ = 1, T ∞ = 1, β = 5, γ = 1.4 and r 2 = |x| 2 . Here, we use the nodal scheme of Rusanov (version 2) with third order reconstruction without any limitor.

Figure 24 :

 24 Figure 24: Error study (ρ, u1, E): Polygonal version of nodal Rusanov scheme (version 2) is only second order in space in L1 norm (left) and L∞ norm (right)

Figure 25 :

 25 Figure 25: Error study (ρ, u1, E): degenerate conic version of nodal Rusanov scheme (version 2) is nearly third order in L1 norm (left) and L∞ norm (right)

Figure 30 :

 30 Figure 30: Top: Comparison of nodal Roe Viscosity form (Version 2) scheme (left) and Rusanov (Version 2) scheme (right) with third order reconstruction and limitation (4.123) and (4.124) (4.122) (4.136), Center : zoom around the contact discontinuity, Bottom : comparison from first to third order on 200 cells of each of the two schemes.

Figure 31 :

 31 Figure 31: Third order composite nodal schemes. Left: VFFC-C (Version 2) and Right: Rusanov (Version 2).

Figure 32 :

 32 Figure 32: Third order composite nodal schemes. Left: nodal Roe Viscosity Form (Version 2) and Right: nodal Roe Flux Form (Version 2).

Figure 33 :

 33 Figure 33: Zoom around x 1 = .7 and superposition of the two Figures 32, results are close. 69

Figure 34 :Figure 35 :

 3435 Figure 34: Numerical convergence of third order composite nodal schemes: VFFC-C (Version 2) on Left and Rusanov (Version 2) on Right for the Osher-Shu test case.

  4) (1.0, 0.1, -0.3, 1.0) (0.8, 0.1, -0.3, 0.4) (0.5313, 0.1, 0.4276, 0.4)where each quadrant is described by the two by two matrix above (his first row mean for the two up left and right quadrant, second row mean for the two down left and right quadrant), and the final time is taken as T= 0.2.

Figure 36 :

 36 Figure 36: Density on 400x400 cells: Top (reference) and Bottom: composite nodal Roe Flux Form (Version 2) (Left) / composite nodal Rusanov (Version 2) (Right).

5. 3 . 4

 34 Radial Sod shock tube : Nodal flux on circular mesh (36x200) at t=0.[START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF]

Figure 39 :Figure 40 :Figure 41 :

 394041 Figure 39: Circular mesh and numerical representation of density and velocity .

Figure 42 :

 42 Figure 42: Density with third order nodal Roe (viscosity form) on degenerate conic on 39616 triangular cells (Top) and associated dual Voronoi mesh (Bottom) at final time t=4.

Figure 43 :

 43 Figure 43: Density with third order nodal Rusanov (version 2) on degenerate conic on cartesian (300x100) cells at final time t=4.

Figure 44 :

 44 Figure 44: Density: third order nodal Rusanov flux on degenerate conic triangular cells (Left) and conical edges (randomized ω ∈ (0, 10): Color code for edge type: Green for ellipses, Red for parabolas, Blue for hyperbolas) (Right)

Figure 45 :

 45 Figure 45: Density: third order nodal Rusanov flux on degenerate conic unstructured triangular cells (Left 19224 cells and Right 40201 cells), at final time t=0.2. 77

Figure 46 :

 46 Figure 46: Density with third order nodal Rusanov Version 2 (Left) and nodal Roe (flux form) (Right) on degenerate conic unstructured triangular (40201) cells

Figure 47 :

 47 Figure 47: Density with third order nodal Rusanov Version 2 (Left) and nodal Roe (flux form) (Right) on degenerate conic unstructured polygonal Voronoi like cells (dual cells of triangular mesh in Figure 46)

Figure 48 :

 48 Figure 48: Density with third order nodal Rusanov flux (Version 2) on left and nodal Roe (flux form) on right. Degenerate conic unstructured quadrilateral cells are mainly used, only boundary holes edges are truly curved.

Figure 49 :

 49 Figure 49: Magnitude Velocity (top) and Pressure (bottom) with third order nodal Rusanov flux (Version 2) on left and nodal Roe (flux form) on right. Degenerate conic unstructured quadrilateral cells are mainly used, only boundary holes edges are truly curved.

Figure 51 :

 51 Figure 51: Comparison of second and third order (Left Above) and different zooms (Right Above, and full Bottom): Limited third order nodal scheme lake of dissipation to capture entropic solution (as mentionned also in citeBizar)

6 Figure 52 :

 652 Figure 52: Left: Evolution of coefficients α i with respect to value of parameter ω (α 1/2 is the top curve). Right zoom around ω = 1 for the Simpson's rule. In view of Figure 52, we have the following: Remark B.2.-The degenerate conical case corresponding to ω = 0 gives: α 0 = α 1 = 1 -π 4 ( .2146018) and α 1/2 = π 2 -1( .5707963). The limit ω → ∞ gives α 0 = α 1 = 0 and α 1/2 = 1.

1 / 2 =

 12 (h(ω) -1)(ω + 1) (B.8)

  3. Relationship between "GLACE type vectors" (2.18) (2.19) cf Figure (9) and "EUCCLHYD type vectors" defined by (2.22) cf Figure (10):

	Cr,ω Cr+1/2,ω j	=	j Ñr+1/2,ω-j	, +	Ñr+1/2,ω+ j	.	(2.25)
	and also:
	Cdof,ω				j	.		(2.26)
	Remark 2.4. Property for planar degenerate conic case.
	By a direct computation on (2.15), we have the following:

j = Ñr,ωj + Ñr,ω+ j = Ñdof,ωj + Ñdof,ω+

Table 1 :

 1 Special values of θ and local weights associated to Ĉr

	1)):

  Figure 21: Sod shock tube with sonic point. Left: Entropic correction of first order nodal VFFC is done using first order nodal Rusanov it calls VFFC-C (see Definition 3.41). Right: Comparison with the Rusanov nodal scheme (Version 2).

	1.1										1.1
								dens_vffc_two_800									dens_vffc_two_800
								dens_vffc-c_two_800									dens_vffc-c_two_800
	1								density_ref	1								dens_rusanov_two_800
																				density_ref
	0.9										0.9
	0.8										0.8
	0.7										0.7
	0.6										0.6
	0.5										0.5
	0.4										0.4
	0.3										0.3
	0.2										0.2
	0.1										0.1
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
								dens_vffc_two_800									dens_vffc_two_800
								dens_vffc-c_two_800									dens_vffc-c_two_800
								dens_rusanov_two_800									dens_rusanov_two_800
	0.85								density_ref	0.85							dens_vffc-cc_two_800 density_ref
	0.8										0.8
	0.75										0.75									3.167)
	0.7										0.7
	0.65										0.65
	0.6										0.6
		0.46	0.48		0.5		0.52		0.54	0.56			0.46		0.48	0.5		0.52		0.54	0.56

  4, we have taken the rectangular domain [-1, 1] × [-1, 1] (periodic boundary condition) and a final time is T = 4.

	1.1 1 0.9		0.8 0.9 1 1.1 0.8 0.9 1 1.1 0.8 0.9 1 1.1					rho_ref rho_o3_rusanov_two_800 rho_o3_rusanov_two_400 rho_o3_rusanov_two_200 rho_o3_rusanov_two_100 rho_ref rho_o3_vffc_two_800 rho_o3_vffc_two_400 rho_o3_vffc_two_200 rho_ref rho_o3_vffc_800 rho_o3_vffc_400 rho_o3_vffc_200 rho_o3_vffc_100 rho_ref rho_o3_roe_two_800 rho_o3_roe_two_400 rho_o3_roe_two_200 rho_o3_roe_two_100 rho_o2_roe_two_800 rho_o2_roe_two_400		0.8 0.9 1 1.1 0.8 0.9 1 1.1 1 1.1 1.1 1 0.9 0.9 0.8					rho_ref rho_o3_roeflux_two_800 rho_o3_roeflux_two_400 rho_o3_roeflux_two_200 rho_o3_roeflux_two_100 rho_ref rho_o3_roe_two_800 rho_o3_roe_two_400 rho_o3_roe_two_200 rho_o3_roe_two_100 rho_ref rho_o3_rusanov_800 rho_o3_rusanov_two_800 rho_ref rho_o3_rusanov_400 rho_o3_rusanov_200 rho_o3_rusanov_two_400 rho_o3_rusanov_100 rho_o3_rusanov_two_200 rho_o3_rusanov_two_100 rho_o2_rusanov_two_800 rho_o2_rusanov_two_400
	0.8		0.7 0.7 0.7						rho_o2_roe_two_200 rho_o2_roe_two_100		0.7 0.7 0.7 0.8					rho_o2_rusanov_two_200 rho_o2_rusanov_two_100
	0.7		0.6 0.6 0.6								0.6 0.6 0.6 0.7
			0.5 0.5 0.5								0.5 0.5 0.5
	0.6										0.6
			0.4 0.4 0.4								0.4 0.4 0.4
	0.5		0.3 0.3 0.3								0.3 0.3 0.3 0.5
	0.4		0.2 0.2 0.2		0.100000						0.100000 0.2 0.2 0.2 0.4
	0.3		0.1 0.1 0.1	0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1 1	0.1 0.1 0.1 0.3	0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1
	0.2		0.9 1 0.9 1 0.9 1		0.010000 vit_o3_rusanov_two_800 vit_ref vit_o3_rusanov_two_400 vit_ref vit_o3_vffc_two_800 vit_o3_vffc_two_400 vit_ref vit_o3_vffc_800 vit_o3_vffc_400					0.010000 0.9 1 0.9 1 0.9 1 0.2	vit_ref vit_o3_roeflux_two_800 vit_o3_roeflux_two_400 vit_ref vit_o3_roe_two_800 vit_o3_roe_two_400 vit_o3_rusanov_400 vit_ref vit_o3_rusanov_800
	0.1	0	0.8 0.8 0.8	0.001000 vit_o3_rusanov_two_200 vit_o3_rusanov_two_100 vit_o3_vffc_two_200 vit_o3_vffc_two_100 vit_o3_vffc_200 vit_o3_vffc_100 0.2 0.4	0.6		0.8	1	0.001000 0.8 0.8 0.8 0.1 0	vit_o3_roeflux_two_200 vit_o3_roeflux_two_100 vit_o3_roe_two_200 vit_o3_roe_two_100 vit_o3_rusanov_200 vit_o3_rusanov_100 0.2	0.4	0.6	0.8	1
			0.7 0.7 0.7								0.7 0.7 0.7
	0.5	0.5 0.6 0.5 0.6 0.5 0.6		0.000100				rho_ref rho_o3_roe_two_800 rho_o3_roe_two_400 rho_o3_roe_two_200		0.000100 0.5 0.6 0.5 0.6 0.6 0.5 0.5				rho_ref rho_o3_rusanov_two_800 rho_o3_rusanov_two_400 rho_o3_rusanov_two_200
	0.45	0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2		0.000010 0.01 0.100000				order 2 order 3 error rho_o3_roe_two_100 rho_o2_roe_two_800 rho_o2_roe_two_400 rho_o2_roe_two_200 rho_o2_roe_two_100	0.000010 0.4 0.4 0.4 0.3 0.3 0.3 0.45 0.01 0.2 0.2 0.2 0.100000			order 2 order 3 error	rho_o3_rusanov_two_100 rho_o2_rusanov_two_800 rho_o2_rusanov_two_400 rho_o2_rusanov_two_200 rho_o2_rusanov_two_100
	0.4	0.1 0.1 0.1								0.1 0.1 0.1 0.4
			0 0 0		0.010000						0.010000 0 0 0
	0.35	-0.1 -0.1 -0.1							-0.1 -0.1 -0.1 0.35
				0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1 1		0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1
	0.3	3 3 3		0.001000 eint_ref eint_o3_rusanov_two_800 eint_ref eint_o3_vffc_two_800 eint_ref eint_o3_vffc_800					3 3 3 0.001000 0.3	eint_ref eint_o3_roeflux_two_800 eint_ref eint_o3_roe_two_800 eint_ref eint_o3_rusanov_800
					eint_o3_rusanov_two_400 eint_o3_vffc_two_400 eint_o3_vffc_400							eint_o3_roeflux_two_400 eint_o3_roe_two_400 eint_o3_rusanov_400
			2.8 2.8 2.8		eint_o3_rusanov_two_200 eint_o3_rusanov_two_100 eint_o3_vffc_two_200 eint_o3_vffc_two_100 eint_o3_vffc_200 eint_o3_vffc_100					2.8 2.8 2.8		eint_o3_roeflux_two_200 eint_o3_roeflux_two_100 eint_o3_roe_two_200 eint_o3_roe_two_100 eint_o3_rusanov_200 eint_o3_rusanov_100
	0.25	2.6 2.6 2.6		0.000100					0.000100 2.6 2.6 2.6 0.25
									order 2							order 2
									order 3							order 3
	0.2	2.4 2.4 2.4		0.000010 0.01				error		0.000010 0.01 2.4 2.4 2.4 0.2			error
			0.5		0.100000 0.55	0.6	0.65	0.7	0.75	0.8	0.100000	0.5	0.55	0.6	0.65	0.7	0.75	0.8
	1.1		2.2 2.2 2.2								2.2 2.2 2.2 1.1
									rho_ref							rho_ref
									rho_o3_roe_two_200							rho_o3_rusanov_two_200
	1		2 2 2		0.010000				rho_o2_roe_two_200 rho_o1_roe_two_200		2 2 2 0.010000 1				rho_o2_rusanov_two_200 rho_o1_rusanov_two_200
	0.9		1.8 1.8 1.8								1.8 1.8 1.8 0.9
					0.001000						0.001000
	0.8										0.8
			1.6 1.6 1.6								1.6 1.6 1.6
				0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1 1		0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1
	0.7		1.1 1.1 1.1		0.000100				pres_ref pres_ref pres_ref		0.000100 1.1 1.1 1.1 0.7				pres_ref pres_ref pres_ref
	0.5 0.6		0.8 0.8 0.8 0.9 0.9 0.9 1 1 1		0.000010 0.01			order 2 order 3 error pres_o3_rusanov_two_100 pres_o3_vffc_two_100 pres_o3_vffc_100 pres_o3_rusanov_two_800 pres_o3_rusanov_two_400 pres_o3_rusanov_two_200 pres_o3_vffc_two_800 pres_o3_vffc_two_400 pres_o3_vffc_two_200 pres_o3_vffc_800 pres_o3_vffc_400 pres_o3_vffc_200		0.000010 0.8 0.8 0.8 0.01 0.5 0.9 0.9 0.9 1 1 1 0.6			order 3 error pres_o3_roeflux_two_100 pres_o3_roe_two_100 pres_o3_rusanov_100 order 2 pres_o3_roeflux_two_800 pres_o3_roeflux_two_400 pres_o3_roeflux_two_200 pres_o3_roe_two_800 pres_o3_roe_two_400 pres_o3_roe_two_200 pres_o3_rusanov_800 pres_o3_rusanov_400 pres_o3_rusanov_200
	0.4		0.7 0.7 0.7								0.7 0.7 0.7 0.4
			Figure 26: Error study on degenerate conic schemes for regular periodic solution (5.2) (L1 norm on Left
	0.3		and L∞ on Right). Top : Nodal Rusanov (version 2), Center : Nodal Roe (Viscosity Form), Bottom : Nodal 0.6 0.6 0.6 0.6 0.6 0.6 0.3
	0.2		0.5 0.5 0.5 Roe (Flux Form). 0.4 0.4 0.4					0.5 0.5 0.5 0.4 0.4 0.4 0.2
	0.1	0	Remark 5.1. 0.3 0.3 0.3 0.2	0.3 0.3 0.3 0.1 1. The errors for both problems (25) (26) seems to indicate that conic degenerate schemes are 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8	1
			0.2 0.2 0.2		quasi third order.					0.2 0.2 0.2
			0.1 0.1 0.1								0.1 0.1 0.1
			0 0 0								0 0 0
				0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1 1		0 0 0	0.2 0.2 0.2	0.4 0.4 0.4		0.6 0.6 0.6	0.8 0.8 0.8	1 1
			Figure 27: Comparison of nodal Rusanov (Version 2) scheme (left) and nodal Roe Flux Form (Version 2) scheme Figure 28: Comparison of nodal Rusanov (Version 2) scheme (left) and nodal Roe Viscosity Form (Version 2) Figure 29: Comparison of nodal VFFC (Version 1) scheme (left) and Rusanov (Version 1) scheme (right) with
			(right) with third order reconstruction and limitation (4.123) and (4.124) (4.122) (4.136) scheme (right) with third order reconstruction and limitation (4.123) and (4.124) (4.122) (4.136) third order reconstruction and limitation (4.123) and (4.124) (4.122) (4.136)
										65 66 67

  Figure50: Bizzarium test case[START_REF] Heuzé | Dissipative issue of high-order shock capturing schemes with non-convex equations of state[END_REF] with second order (Left) and third order (Right) nodal Rusanov (Version 2) scheme.

						dens_rusanov_two_o3_16000 dens_rusanov_two_o2_16000	14000						dens_rusanov_two_o3_16000 dens_rusanov_two_o2_16000
								dens_ref									dens_ref
										13900
										13800
										13700
										13600
										13500
										13400
										13300
	0		0.2		0.4	0.6		0.8		1		0.2	0.22		0.24	0.26	0.28	0.3
						dens_rusanov_two_o3_16000	12800						dens_rusanov_two_o3_16000
						dens_rusanov_two_o2_16000							dens_rusanov_two_o2_16000
								dens_ref									dens_ref
										12600
										12400
										12200
										12000
										11800
	0.23	0.24	0.25	0.26	0.27	0.28	0.29	0.3	0.31	0.32	0.74	0.76	0.78	0.8	0.82	0.84	0.86	0.88	0.9
										1.2x10 11
																pressure_rusanov_two_o3_1000
																pressure_rusanov_two_o3_4000
																pressure_rusanov_two_o3_16000
																	pressure_ref
											0		0.2	0.4		0.6	0.8		1

we could also deal with a more general inequality K 1 ≤ P

(α j ) ≤ K 2

The pure nodal polygonal scheme do not permits to obtain third order only composite nodal (degenerate conic) scheme does.

Indirect limitor for velocity reconstruction for second and third order

In view of physical constrain for the continuous system ρ ≥ 0, ≥ 0, we can decide to apply a direct limitation for ρ, for both second and third order reconstruction. We give thereafter some hypothesis:

-Suppose that any conservative quantity U = (ρ, ρU, ρE) is approximated by a same type of reconstructed function (here polynomials of degree n : P n j (x, U )).

-Suppose there exists two polynomials of degree n, noted P n and P n (see previous section for order 2 (n=1, (4.74) (4.75)) or for order 3 (n=2, (4.114) (4.115))) such that:

, homogeneous of degree one w.r.t second variable

homogeneous of degree two w.r.t second variable.

The notation [q] j i means that the terms describing q involves derivatives only from degree i to j.

Let take λ , λ U limitors for and U, we can write a limited version of reconstructed massic internal energy (at least for n = 1, 2):

Now, by choosing in (4.121)

we obtain an induced limitor for velocity (build from that of internal massic energy). This choice naturally define an equal amount of limitation between kinetic energy e c (e c = 1 2 |U| 2 ) and internal energy at any point of the domain mimeticaly to relation E = e c + . Proposition 4.1. We can obtain a limited reconstruction for internal energy (directly) and velocity U (indirectly) :

, with λ U given by (4.122). (4.124)

-Each of (4.123) and (4.124) may be written as R n,lim j

where Rj (x, Q) is the strict "high order" extension. Taking λ j = 1 gives the unlimited version whereas λ j = 0 reduce to first order R n,lim j

-Classical scalar limitors can be used for λ , Dukowicz, Barth-Jespersen [START_REF] Barth | The design and application of upwind schemes on unstructured meshes, chapter 6[END_REF], MLP [START_REF] Park | Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids[END_REF] or min-mod (at least obtaining positivity R n j (x, ) ≥ 0 for all x dof points of Ω j in (4.123) but also in practice min 4.124), we can easily adapt an APITALI algorithm considering sequence of real valued degree of reconstruction [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF]. In practice, using the a priori limitation (4.123) and (4.124) (4.122) (4.136) whatever the evaluation of reconstructed functions at dof in each cell make the numerical composite flux well defined. Nevertheless, for a high order flux (two or three here), the updated unknowns U n+1 j in (4.137) may not be in the admissible set (ρ n+1 j ≥ 0, n+1 j ≥ 0). We recall the definition of hierarchical degree limitation for a polynomial reconstruction (see [START_REF] Bernard-Champmartin | Stabilité locale et montée en ordre pour la reconstruction de quantités volumes finis sur maillages coniques non-structurés en dimension 2[END_REF]) of degree n (here for exemple the density ρ) :

where each differential operator is evaluated at x j , the (α l j ) n l=1 are cell scalar values belonging to [0, 1]. Let (α l j ) n l=1 coefficients that must belong to the following set:

-

(4.140)

First, we define the application:

.141)

Let A his inverse function going from A n in [0, n], for all (α 1 j , . . . , α n j ) in A n may be defined:

We recall some ingredients discussed and introduced in ([7]):

Definition 4.4. Let a polynomial P n j defined by (4.139) whose coefficients (α l j ) n l=1 are in A n . One call degree of limitation (hierarchical) of P n j , the real:

The hierarchical limitation begin on highest degree and proceed with decreasing degree. For the active index of limitation l, we search for the coefficient α l j ∈]0, 1] nearest from 1 while ensuring stability. We introduce a sampling of size N * (l) ≥ 1 of [0, 1] for each term of degree l:

Non smooth solutions on one-dimensional test cases

In this section, for non smooth solution, we will essentially use the third order reconstruction with limitation (4.123) and (4.124) (4.122) (4.136) .

Sod's shock tube problem

This famous test problem is to test the capability of numerical methods to resolve discontinuous flows for the Euler equations. The initial data is given by: (ρ 0 , u 0 , p 0 ) = (1, 0, 1), 0 ≤ x < 0.5 (0.125, 0, 0.1), 0.5 < x ≤ 1.

The solution consists in a left rarefaction, a contact and a right shock wave. In the following, we compare different nodal composite schemes.

In this test case, the initial data is given by: (ρ 0 , u 1,0 , u 2,0 , p 0 ) = (1, 0, 0.7276, 1) (0.5313, 0, 0, 0.4) (0.8, 0, 0, 1) (1, 0, 0.7276, 1)

Lax-Liu (Case 3)

In this test case, the initial data is given by: (ρ 0 , u 1,0 , u 2,0 , p 0 ) = (0.138, 1.206, 1.206, 0.029) (1.5, 0, 0, 1.5) (0.5323, 1.206, 0, 0.3) (0.5323, 0, 1.206, 0.3) hence in some sens, our composite finite volume solvers make a constructive bridge between edge finite volume scheme (one dimensionnal) to nodal finite volume scheme (trully two dimensionnal) both on straight and curved meshes.

We have also introduced some overall third order in space and time schemes based on nonlinear reconstruction (rationnal function) of massic quantities ( and U). We intensively use Leibniz formula on momentum ρU and total energy ρE to deduce them for velocity and massic internal energy. We have shown that it is possible to use a priori-limitation strategy on that indirectly gives the limitation of U. An immediate consequence is that limitation on massic internal energy and massic kinetic energy are of same scale, so that it mimics the algebraic relation E = + 1 2 |U| 2 . We have compared our new composite schemes on different test-cases from first to third order. We obtain quite similar results as those of litterature on smooth and non regular solutions for both perfect gas and Mie-Grüneisen (non convex) EOS.

Many possible extensions are investigated:

1. Running more challenging test cases (with vacuum), both on straight and curved boundaries (and also straight/curved meshes), and more comparisons (with other numerical nodal fluxes such as HLL, and some variants). Compare our composite scheme with pure edge schemes and pure nodal schemes (polygonal).

2. Adaptation to the multi-material compressible flows.

3. Coupling the TVB non-linear reconstruction [START_REF] Deng | Limiter-free discontinuity-capturing scheme for compressible gas dynamics with reactive fronts[END_REF] with our proposed third order non-polynomial (rational) one.

4. An easy adaptation to the formalism of Ñdof,ω 5. We plan to use this CARS flux to deal with well-balanced discretization of source terms [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources[END_REF]. We also want to use them in the context of Asymptotic Preserving schemes and compare with other AP conical schemes like in [START_REF] Blanc | Asymptotic preserving schemes on conical unstructured 2d meshes[END_REF].

6. Use them for Discontinuous Galerkin discretization [START_REF] Dumbser | ADER discontinuous Galerkin schemes for aeroacoustics[END_REF], or Residual Distributive schemes [START_REF] Abgrall | High-order residual distribution scheme for the time-dependent Euler equations of fluid dynamics[END_REF].

7. Adaptation to Lagrange Update hydrodynamic equations.

Appendices

A Local behavior of numerical nodal flux

A more local formulation (in link to local residual):

(A.1) More specifically, the dependency is neglected in the case of Version 2 of both nodal flux (3.34) and nodal viscosity flux (3.82), we recall that in these two cases, only two state average between U n j and U n k are considered. In this case, we note:

In contrast that for Version 1 of these schemes ((3.33) and (3.81)), all states around dof have an impact to the exchange of the two cells j and k. In this case, we note:

Let's define the contributive part of the solution update U n+1 in terms coming from nodal flux at dof. We need to define two versions corresponding to Version 1 or Version 2 in both nodal flux schemes (see Definition For each dof,

For each dof, it is linked to two states interactions along a normal:

The schemes associated to Version 1 or Version 2 rewrite like this:

Definition A.2. Invariant domain preserving nodal flux (Version 1) Let an invariant domain C ⊆ A adm , a numerical flux is invariant domain preserving if:

)) ∈ C.

B Newton Cotes quadrature formula

Let f (q) a given function (at least continuous) defined on [0, 1] with value in normed vectorial space V (ex. R p ):

Let be P 2,ω a quadratic rational (Bezier) interpolant of function f at endpoints (f 0 , f 1 ) and middle point (f 1/2 ): P 2,ω (q) = f 0 (1q)(1 -2q) + 2(ω + 1)q(1q)f 1/2 + q(2q -1)f 1 (1q) 2 + 2ωq(1q) + q 2 , q ∈ [0, 1].

(B.2)

We want to give a quadrature formula in the sense of Newton-Cotes quadrature formula (note that P 2,ω (q) is really an interpolant at q = {0, 1 2 , 1} of f ). Hence, we replace the exact problem (B.1) by the integral of interpolant (B.2): 1 0 P 2,ω (q)dq (B.3)

Proposition B.1. The Newton-Cotes formula corresponding to a quadratic rationnal interpolation of a continous function f on [0, 1] (with value in V) is given by:

where α ω i are such that: α ω i ≥ 0, i α ω i = 1 and verifies: