S Diot

S Clain
email: clain@math.uminho.pt

R Loubère

Improved Detection Criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials

Keywords: Finite Volume, high-order, conservation law, polynomial reconstruction, limitation, polygonal, non-conformal, unstructured, Euler, MOOD, positivity-preserving

This paper extends the MOOD method proposed by the authors in ["A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)", J. Comput. Phys. 230, pp 4028-4050, (2011)], along two complementary axes: extension to very high-order polynomial reconstruction on nonconformal unstructured meshes and new Detection Criteria. The former is a natural extension of the previous cited work which confirms the good behavior of the MOOD method. The latter is a necessary brick to overcome limitations of the Discrete Maximum Principle used in the previous work. Numerical results on advection problems and hydrodynamics Euler equations are presented to show that the MOOD method is effectively high-order (up to sixth-order), intrinsically positivity-preserving on hydrodynamics test cases and computationally efficient.

Introduction

In a recent paper [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF], an original high-order method, namely the Multidimensional Optimal Order Detection (MOOD) method, has been introduced to provide up to third-order approximations to hyperbolic scalar or vectorial solutions for two-dimensional geometry. The present article deals with new extensions of the method to general unstructured 2D meshes and to sixth-order convergence in space. Classical high-order reconstructions such as MUSCL or ENO/WENO methods are based on an a priori limiting procedure to achieve stability property. The MOOD method follows a fundamentally different way since the limiting procedure (polynomial degree reduction for instance) is achieved a posteriori and provides the optimal local polynomial reconstruction which satisfies given stability criteria.

The quest [START_REF] Van Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF] of the (very) high-order schemes starts in the early 70's with the pioneer works of Van-Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme II. Monotonicity and conservation combined in a second-order scheme[END_REF] and Kolgan [START_REF] Kolgan | Application of the minimum-derivative principle in the construction of finite-difference schemes for numerical analysis of discontinuous solutions in gas dynamics[END_REF][START_REF] Kolgan | Finite-difference schemes for computation of three dimensional solutions of gas dynamics and calculation of a flow over a body under an angle of attack[END_REF][START_REF] Kolgan | Application of the principle of minimizing the derivative to the construction of finite-difference schemes for computing discontinuous solutions of gas dynamics[END_REF]. Since this date, a large literature was dedicated to the limited reconstruction methods for structured and unstructured meshes. Several strategies became very popular due to their intrinsic simplicity such that the MUSCL method [START_REF] Barth | Numerical methods for conservation laws on structured and unstructured meshes[END_REF][START_REF] Buffard | Monoslope and Multislope MUSCL Methods for unstructured meshes[END_REF][START_REF] Clain | L ∞ stability of the MUSCL methods[END_REF][START_REF] Hubbard | Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids[END_REF][START_REF] Leveque | High-Resolution Conservative Algorithms for Advection in Incompressible Flow[END_REF][START_REF] Park | Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids[END_REF] or their efficiency to achieve very high-order accuracy such that the ENO/WENO method [START_REF] Abgrall | On Essentially Non-oscillatory Schemes on Unstructured Meshes: Analysis and Implementation[END_REF][START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Friedrich | Weighted Essentially Non-Oscillatory Schemes for the Interpolation of Mean Values on Unstructured Grids[END_REF][START_REF] Ollivier-Gooch | Quasi-ENO Schemes for Unstructured Meshes Based on Unlimited Data-Dependent Least-Squares Reconstruction[END_REF][START_REF] Shu | Essentially non-oscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws[END_REF][START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing scheme[END_REF][START_REF] Wolf | High-order ENO and WENO schemes for unstructured grids[END_REF][START_REF] Harten | Uniformly highorder accurate nonoscillatory schemes III[END_REF][START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF][START_REF] Shi | A technique of treating negative weights in WENO schemes[END_REF][START_REF] Zhang | Third-order WENO scheme on three dimensional tetrahedral meshes[END_REF][START_REF] Tsoutsanis | WENO schemes on arbitrary mixedelement unstructured meshes in three space dimensions[END_REF], the Discontinuous Galerkin method [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: General framework[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One dimensional systems[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case[END_REF][START_REF] Cockburn | The Runge-Kutta Discontinuous Galerkin Method for Conservation Laws V: Multidimensional Systems[END_REF], the ADER method [START_REF] Dumbser | ADER schemes on unstructured meshes for nonconservative hyperbolic systems: Applications to geophysical flows[END_REF][START_REF] Titarev | ADER schemes for three-dimensional non-linear hyperbolic systems[END_REF][START_REF] Toro | ADER finite volume schemes for nonlinear reactiondiffusion equations[END_REF][START_REF] Dumbser | Arbitrary High Order Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Linear Hyperbolic Systems[END_REF][START_REF] Dumbser | Quadrature-Free Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Nonlinear Hyperbolic Systems[END_REF], the Residual Distribution Scheme [START_REF] Abgrall | Essentially non-oscillatory Residual Distribution schemes for hyperbolic problems[END_REF][START_REF] Csík | A Conservative Formulation of the Multidimensional Upwind Residual Distribution Schemes for General Nonlinear Conservation Laws[END_REF][START_REF] Ricchiuto | Stabilized residual distribution for shallow water simulations[END_REF] and the spectral method [START_REF] Harris | Efficient quadrature-free high-order spectral volume method on unstructured grids: Theory and 2D implementation[END_REF][START_REF] Wang | Spectral (finite) volume method for conservation laws on unstructured grids: basic formulation[END_REF][START_REF] Wang | Spectral (finite) volume method for conservation laws on unstructured grids: extention to two dimensional scalar equation[END_REF].

While second-order methods do not require particular cautions, dealing with higher-order methods leads to at least three specific difficulties which, up to our knowledge, are not always clearly identified. First point one should not consider the mean value of a function equivalent to the cell centroid value as it is often done in the MUSCL community. The point is straightforward to overcome but important to notice for newcomers in the field of higher-order numerical schemes. Second point, for vectorial problems the reconstruction process must be done on mean values of the conservative variables and not on non-linear combinations of them. This point is often implied in the classical ENO/WENO papers but is rarely clearly stated and this may mislead newcomers in the high-order community because the order of accuracy discrepancy can be missed depending on the numerical tests used. Contrarily one proposes the isentropic vortex in motion test case to numerically prove that if the primitive variables are used for the reconstruction process then very high-order of accuracy cannot be reached. Third point the Discrete Maximum Principle property on mean values should not be used anymore as a guide line for limitation. We propose in this paper to overcome this difficulty by a new limiting criteria (or Detection Criteria in the MOOD jargon) adapted to provide a full high-order method still maintaining robust stability. Simple examples are introduced within the text when some difficulties related to these points are to be expected.

The basic idea of the MOOD method consists of determining the higher polynomial degree of each local cell still satisfying some stability restrictions. To this end, an iterative process is developed. We perform a local polynomial reconstruction of degree d i for each cell K i at the current time t n and compute a candidate solution at time t n+1 without any limiting features. Then a detecting procedure is carried out to check the cells which do not respect the stability criteria and we reduce the local polynomial degree to obtain a better stability. We state that the method is a posteriori since the limiting procedure (namely the polynomial degree reduction) is performed after the candidate solution computation. Such a situation is very useful to test the admissibility of the solution. Furthermore, one has to carry out the limiting algorithm if, and only if, it is necessary while the traditional a priori method performs unnecessary limitation.

In this paper we propose extensions of the MOOD method which take into account the three difficulties mentioned above. More precisely different detection processes both for the advection and hydrodynamics equations are developed. We numerically prove that these detection processes provide the effective higher-order of accuracy on smooth profiles (up to sixth-order). Moreover we show that for the hydrodynamics equations the method is positivity-preserving by construction and we numerically observe this behavior. The test case have been carried out on non-regular, polygonal and non-conformal meshes and the last test case of the paper show the ability of the MOOD method to simulate complex physics from an experimental set-up of the impact of a shock wave on a cylindrical cavity.

The paper is organized as follows. Section 2 is dedicated to the generic framework used to describe the MOOD method where the high-order finite volume scheme is presented. Several obstacles to achieve high-order reconstruction are pointed out and the polynomial reconstruction based on the mean value approximation is detailed. In Section 3, we introduce new criteria to obtain very high-order accurate schemes still preserving local stability. To show the MOOD method efficiency, numerical tests both for the scalar and the vectorial case are carried out in Section 4. We mainly focus on the method accuracy and its robustness. We draw some remarks and future developments in the last section.

The MOOD method

General concept

The MOOD method is a generic procedure that solves multidimensional hyperbolic system of equations on an unstructured grid in the Eulerian framework. Given different numerical finite volume schemes the MOOD method provides an optimal choice for each computational cell by mitigating accuracy vs robustness. From an abstract point of view the MOOD algorithm involves two main ingredients: An ordered list of numerical schemes and a set of constraints with detection criteria which defines the desirable properties the numerical solution should have. The over-topping numerical scheme represents the best scheme one would like to employ. Usually this scheme is the most accurate but less robust one. At the very end of the list lays the least accurate but more robust scheme which is assumed to be satisfactory in all possible situations due to the stabilization effect generated by its intrinsic numerical dissipation. In this paper the list is composed of a robust first-order scheme (an upwind or a Rusanov, HLL, HLLC scheme as instance) while several second or higher-order schemes using polynomial reconstructions compose an ordered list of desirable schemes (see Fig. 1 for instance). The second ingredient is the detecting procedure of a set of constraints which determine the local eligibility of the solution for each cell.

We recall that discontinuous solutions may not be handled with high-order reconstructions since local spurious and unphysical oscillations may take place. The low-order numerical scheme should be used to prevent the numerical approximations from oscillating and force to respect some constraints or mathematical properties that depend on equations under consideration. The numerical solution is considered as eligible if it fulfills given properties. As instance the positivity of certain variables such as density or pressure in hydrodynamics equations or the Discrete Maximum Principle for advection equation shall be considered.

In this paper the k-th numerical scheme of the list is a finite volume scheme using unlimited piecewise polynomial reconstruction of degree k. Ultimately this scheme has a k+1th-order of accuracy for smooth solutions. Consequently the LO scheme is the generic firs-order finite volume scheme and the HO-1 scheme corresponds to an unlimited MUSCL method.

The core of the MOOD method is a loop over the cells to determine the optimal polynomial degree one can safely use to produce an eligible numerical solution. It amounts to select a numerical method in the ordered list of Fig. 1.

Low-Order

HO -1 . . .

HO -k

Scheme #k Scheme #(k-1)
HO -k-1 . . .

Increasing scheme robustness

Increasing scheme accuracy Scheme #0

High-Order

Scheme #d max Scheme #1 "Most robust, least accurate" "Least robust, most accurate" Fig. 1. Schematic representation of an ordered list of numerical schemes used in the MOOD method. The bottom scheme is the most robust but least accurate one denoted "Low-Order". All over-topping schemes are successively more accurate but less robust. The MOOD method is designed to choose the more adapted scheme for each cell of the computational domain.

To this end, given a generic cell K i and its neighbor cells K j having edge e ij in common, we first recall two definitions introduced in [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF] and then give a new one to extend the MOOD concept:

• d i is the Cell Polynomial Degree (CellPD) which represents the degree of the polynomial reconstruction on K i .

• d ij = d ji = min(d i , d j ) are the Edge Polynomial Degrees (EdgePD) corre-
sponding to the degrees of the polynomial reconstructions used to compute approximations of the solution on edge e ij . • A is a set of prescribed physical and/or stability constraints. If for each cell K i the mean values of the numerical solution fulfill the constraints then the numerical solution is said to be A-eligible.

The last item concerns the detecting procedure to distinguish if a candidate solution is eligible according to a set of constraints. In practice we decrement the d i for any cell K i which does not respect all the constraints. Such a cell is called problematic. Moreover since neighbor cells fluxes may be affected by this process, the decrementing is spread over the direct neighborhood. Such a polynomial degree decrementing for a problematic cell is repeated up to a d i > 0 for which the set of constraints is fulfilled or to d i = 0. At that ultimate step the robust and diffusive LO scheme is employed and its first-order solution is always taken as valid. In other words unlike traditional high-order schemes (using a priori limiting procedure), we introduce an a posteriori detecting procedure where the decision to alter the polynomial degree is carried out after computing the candidate solution.

We finally highlight that such a procedure may be interpreted as a try and fail algorithm. Such a generic strategy might be adapted to other classes of method such as the Discontinuous Galerkin method and detect the best polynomial degree in each cell or Finite Element method and detect the most appropriate finite element one can employ in a cell.

Framework

Let us consider a generic autonomous hyperbolic equation defined on a domain Ω ⊂ R 2 , t > 0 which casts in the conservative form

∂ t U + ∇ • F (U ) = 0, (1a) U (•, 0) = U 0 , ( 1b
)
where U = U (x, t) is the vector of unknown functions, x = (x, y) denotes a point of Ω, t is the time, F is the physical flux function and U 0 is the initial condition. Boundary conditions shall be prescribed in the following. We assume that the computational domain Ω is a polygonal bounded set of R 2 divided into convex polygonal cells K i , i ∈ E el , c i being the cell centroid and E el the cell index set. For each boundary edge, K i ∩ ∂Ω, we introduce a virtual cell K j with j / ∈ E el which represents the exterior side of Ω and denote by E bd the index set of all virtual cells. E el = E el ∪E bd is the index set of all cells. This notation avoids a special treatment for boundary edges in the scheme, and provides a natural notation for ghost cells should they exist or not.

For each cell K i , one denotes by e ij the common edge between K i and K j , with j ∈ ν(i) ⊂ E el , ν(i) being the index set of all the elements which share an edge with K i . The extended neighborhood is represented by the index set ν(i) ⊂ E el of all K j such that K i ∩ K j = ∅ (see Fig. 2). Moreover |K i | and |e ij | measure the surface of K i and the length of e ij respectively while n ij is the unit outward normal vector to e ij pointing from K i to K j . At last, q r ij , r = 1, ..., R represent the Gaussian quadrature points employed for numerical integration on edge e ij .

The generic first-order explicit finite volume scheme is given by where F(U n i , U n j , n ij ) is a numerical flux which satisfies the classical properties of consistency and monotonicity. To provide higher-order accuracy, we substitute in equation ( 2) the first-order approximation U n i and U n j with better approximations of U at the quadrature points of edge e ij leading to the generic spatial high-order finite volume scheme

U n+1 i = U n i -∆t j∈ν(i) |e ij | |K i | F(U n i , U n j , n ij ), (2)
U n+1 i = U n i -∆t j∈ν(i) |e ij | |K i | R r=1 ξ r F(U n ij,r , U n ji,r , n ij ), (3)
where U n ij,r and U n ji,r , r = 1, ..., R are high-order approximations of U at quadrature points q r ij ∈ e ij , r = 1, ..., R on both sides of edge e ij and ξ r denote the quadrature weights.

For the sake of simplicity, let us write the scheme under the compact form

U n+1 h = U n h + ∆t H R (U n h ), (4)
with

U n h = i∈E el U n i 1I K i the constant piecewise approximation of function U
and operator H R being defined as

H R (U n h ) := - i∈E el   j∈ν(i) |e ij | |K i | R r=1 ξ r F(U n ij,r , U n ji,r , n ij )   1I K i . ( 5
)
Finally to provide a high-order method in time, we use the third-order TVD Runge-Kutta method (RK3, see [START_REF] Shu | Efficient implementation of essentially non-oscillatory shock-capturing scheme[END_REF]) which corresponds to a convex combination of three explicit steps (2) h = U (1) h +∆t H R (U (1) h ) U (3) h = U (2) h +∆t H R ( U (2) h ) [START_REF] Buffard | Monoslope and Multislope MUSCL Methods for unstructured meshes[END_REF] where U (2) h is the convex combination (3U n h + U (2) h )/4.

U n+1 h = U n h + 2U (3) h 3 with              U (1) h = U n h +∆t H R (U n h ) U

Arbitrary degree polynomial reconstruction

In the introduction we have reminded one classical obstacle to reach higherorder of accuracy when polynomial reconstruction is to be used. It is wellknown that the mean value U i of a regular function U on K i is approximated by the value of the solution at the cell centroid, U (c i ), with an error of O(h 2 ) where h represents the characteristic length of the cell. It results that any reconstruction based on geometrical arguments using U (c i ) in place of U i can only provide second-order approximation.

Therefore as classical higher-order finite volume methods the MOOD method is based on polynomial reconstructions from mean values on cells. Let us consider a generic reconstructed polynomial of degree d, given mean values U on a generic cell K, under the form

U (x; d) = U + 1≤|α|≤d R α   (x -c) α - 1 |K| K (x -c) α dx   , ( 7
)
where c is the centroid of K, x a generic point in K and R α are the unknowns polynomial coefficients where α = (α x , α y ) ∈ N 2 is a multi-index with |α| = α x + α y . Note that by construction, the mean value on K of the polynomial function is equal to U since the integral over K of the term between parenthesis in [START_REF] Clain | L ∞ stability of the MUSCL methods[END_REF] vanishes. It thus fulfills the conservation property on K.

There exist several techniques [START_REF] Abgrall | On Essentially Non-oscillatory Schemes on Unstructured Meshes: Analysis and Implementation[END_REF][START_REF] Friedrich | Weighted Essentially Non-Oscillatory Schemes for the Interpolation of Mean Values on Unstructured Grids[END_REF] to determine the coefficients R α . Here, we consider a least square approximation of neighbor mean values U j where K j belongs to a compact stencil S(K). It amounts to minimizing the functional

E( U ) = j∈S(K) ω k 1 |K j | K j U dx -U j 2 , ( 8
)
where ω k are positive weights used to provide a better condition number. In particular, the condition number of the associated linear system dependents on the spatial characteristic length thus we use the solution proposed in [START_REF] Friedrich | Weighted Essentially Non-Oscillatory Schemes for the Interpolation of Mean Values on Unstructured Grids[END_REF] to overcome this problem.

In practice, we do not directly solve the symmetric linear system associated with the minimization problem. Instead we use the technique from [START_REF] Barth | Higher-Order Solution of the Euler Equations on Unstructured Grids Using Quadratic Reconstruction[END_REF][START_REF] Ollivier-Gooch | Quasi-ENO Schemes for Unstructured Meshes Based on Unlimited Data-Dependent Least-Squares Reconstruction[END_REF] where an over-determined linear system is solved in a least-squares sense with a QR decomposition using Householder transformations. The reconstructed polynomial U is thus exact for any polynomial function of degree lower than d which provides the consistency of the reconstruction method and further the status of a (d + 1) th -order numerical method.

Remark 2

In the introduction we have stated that in the general case one should not identify the mean value of a non-linear combination with the nonlinear combination of mean values. Let ρ and φ be two regular functions on cell K i and ρ i , φ i , (ρφ) i , denote their respective exact mean values. A Taylor expansion with respect to the centroid of the cell gives (ρφ

) i = ρ i φ i + O(h 2 ).
For instance let us consider the one-dimensional variables ρ, φ and (ρφ) and their mean values on cell

K 1 = [0, h] ρ(x) = 1 + x, φ(x) = 1 -x, (ρφ)(x) = 1 -x 2 , ρ 1 = 1 + h 2 , φ 1 = 1 - h 2 , (ρφ) 1 = 1 - h 2 3 .
We then deduce that

|(ρφ) 1 -ρ 1 φ 1 | = h 2 /12 leading to a second-order error.
As instance it is well known that for Euler system of equations the non-linear transformation of the conservative mean values into primitive ones introduces a second-order error in the general case.

Algorithm

Let us assume that we have access to a given sequence U n h = (U n i ) i∈E el of mean value approximations at time t n , the goal is to build an eligible sequence

U n+1 h = (U n+1 i
) i∈E el at time t n+1 = t n + ∆t n in the sense that each approximation U n+1 i respects a set of constraints A. We only consider here a forward Euler time step without loss of generality. The MOOD method algorithm is the following: = U h else the solution is recomputed following steps 2., 3. and 4.

Remark 3

Only cells K i where CellPD has been decremented and their neighbors in the compact stencil ν(i) have to be re-updated. Consequently only these cells will have to be checked for the next iterations of the MOOD procedure within the current time step. This dramatically reduces computational cost. [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF] to truncate U i (•; d max ) to obtain lower-order polynomials. However we found that for d max > 2 this technique implies non desirable behavior on discontinuous profiles as the reconstruction stencil remains large. Moreover numerical experiments show that a one-by-one degree decrementing leads to avoidable computational effort since the decrementing procedure is usually performed around discontinuities. We thus slightly modify the decrementing algorithm by jumping from d = d max to d = 2 and then from d = 2 to d = 0 if needed. This also reduces the computational effort while providing equivalent results on a wide range of test cases compared to a one-by-one decrementing.

Remark 4 Since polynomial reconstruction is costly in CPU time and memory, we proposed in

Remark 5 Polynomial reconstruction on boundary cells are treated using ghost cells in order to to be consistent with the prescribed boundary conditions.

The major difficulty remains to determine a list of constraints which both provides a very high accurate solution while avoids numerical artifacts such as spurious oscillations in the vicinity of discontinuity. This is the purpose of the next section.

Detection process

The list of constraints A corresponds to eligible criteria that the numerical approximation has to fulfill. To this end, detection process is necessary to list where the candidate numerical solution fails to respect the constraints. Such process must be very carefully designed to preserve high accuracy for regular solutions whereas discontinuities should be treated with the lower order scheme to avoid non-physical oscillations. The first subsection deals with the advection problem and a new detection process called u2 and based on a smoothness detector. In the second subsection the Euler system is considered: Two detection processes are proposed and we show the positivity-preserving property of the MOOD method.

Advection problem: The u2 detection process

Solutions of autonomous scalar hyperbolic problems satisfy the Maximum Principle property. Such a property is also valid for advection problem with divergence free velocity. Therefore the Discrete Maximum Principle (DMP) seems to be a good candidate to detect problematic cells. Unfortunately, as mentioned in the introduction, the strict DMP applied to mean values reduces the order of accuracy to two (see the appendix for an example), and thus can not be used alone. Classical studies show that the accuracy discrepancy only occurs at extrema [START_REF] Sander | A third-order accurate variation non-expansive difference scheme for single nonlinear conservation law[END_REF][START_REF] Osher | High resolution schemes and the entropy condition[END_REF][START_REF] Jiang | Non-oscillatory central schemes for multidimensional hyperbolic conservative laws[END_REF]. We will then mainly focus on extrema since the DMP detection process is still relevant where the solution is locally monotone. We propose the relaxation of the strict DMP at smooth extrema in order to avoid accuracy discrepancy. This leads to the introduction of an additional procedure to detect smooth extrema. Notice that in (W)ENO type of methods the DMP is not strictly enforced which implies that extrema are well approximated and consequently arbitrary high-order of accuracy is achieved. The first detection criteria is the DMP: No polynomial degree decrementing is performed for cells where the DMP is satisfied. Let us now consider a cell K i where U i does not fulfill the DMP. Two situations may arise whether we deal with a discontinuity or a smooth extrema. The major difficulty is to give a concrete definition of the concept of a smooth extrema from a numerical point of view. Actually a function may be considered irregular for a coarse mesh but regular with a finer one. We try to overcome this difficulty by introducing the following definition. Definition 6 Let K i be a cell and U i = U i (.; 2) a polynomial reconstruction of degree 2 for an underlying function U. We define the second derivatives in x and y directions by

X i = ∂ xx U i ∈ R and Y i = ∂ yy U i ∈ R.
We will refer to these second derivatives as "curvatures". For all cell K j , j ∈ ν(i), we define the maximal and minimal curvatures as

X min i = min j∈ν(i) (X i , X j ) and X max i = max j∈ν(i) (X i , X j ) , Y min i = min j∈ν(i) (Y i , Y j ) and Y max i = max j∈ν(i) (Y i , Y j ) .
We now introduce the new detection criterion to select smooth extrema.

Definition 7 A numerical solution U i in cell K i which violates the DMP is nonetheless eligible if X max i X min i > 0 and Y max i Y min i > 0, ( 9
)
|X min i | |X max i | ≥ 1 -ε i and |Y min i | |Y max i | ≥ 1 -ε i , ( 10
)
where ε i is a cell dependent parameter defined by

ε i = (∆x i ) 1 2m , with ∆x i = |K i | 1 m , m being the spatial dimension (m = 2 here).
Such a detection criterion is motivated by the following considerations. For a given mesh, the solution is locally considered as non-oscillating if condition ( 9) is fulfilled meaning that, at the numerical level, the "curvatures" of the P 2 approximation have the same sign. Moreover for a given mesh, the solution is considered locally C 2 from a numerical point of view if condition ( 10) is fulfilled. The parameter ε is a mesh dependent coefficient which prescribes the tolerance. Such criteria verifies if the "curvatures" are almost identical in the vicinity of cell K i with respect to the local characteristic space length ∆x i . The choice of ε derives from numerous tests. In fact our numerical experiments have shown that ε scales like a cell dependent characteristics length to a power depending on the dimension of space (tests have been carried out in 1D and 2D). It seems to the authors to be the best compromise to gain a very high-order of convergence while maintaining reasonable monotonicity. Finally we remark that at the limit ε = 0 we recover the DMP.

The set of constraints A for advection equation is thus constituted by the DMP relaxed by the smooth extrema detector described above. The detection process is called u2 detection in reference to the second-order derivatives and is summarized in the sequel.

Being given a sequence U h = (U i ) i∈E el , the u2 detection procedure in the case of the advection problem is given by the following algorithm.

1. The DMP criterion is first checked on each cell 9) and [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One dimensional systems[END_REF]. If cell i is not a smooth extrema then d i is decremented, else U i is eligible.

K i min j∈ν(i) (U n i , U n j ) ≤ U i ≤ max j∈ν(i) (U n i , U n j ). ( 11
) 2. If U i does not satisfy (11) then a-Compute X k , Y k for k ∈ ν(i) {i} and coefficient ε i , b-Check criteria (

Euler system: Two detection processes and positivity-preserving

The compressible hydrodynamics Euler system of equations is the following hyperbolic unsteady non-linear system involving conservation of mass, momentum and total energy

∂ t           ρ ρu ρv E           + ∂ x           ρu ρu 2 + p ρuv u(E + p)           + ∂ y           ρv ρuv ρv 2 + p v(E + p)           = 0. ( 12
)
The primitive variables are the density ρ, the velocity U = (u, v) and the pressure p. The pressure is linked to two thermodynamical variables such as density and specific internal energy ε through an Equation Of State (EOS) p = p(ρ, ε). As instance the classical ideal gas law states that p = (γ -1)ρε where γ is the ratio of specifics heats. Moreover the total energy E is such that

E = ρ(ε + 1/2 U 2 ).
Even if the DMP property is used in most of limiting procedures (MUSCL technique as instance), the DMP property does not make sense in the case of the Euler system, for the density or the total energy for instance, since the velocity is not divergence free. Consequently we can not rely only on DMP.

We propose here two detecting procedures which we have been widely experimented and present in the next sections the pros and cons of such procedures.

Physical Admissible Detection (PAD)

The first and minimal detection criteria consists of ensuring the physical meaningfulness of the primitive variables, namely positivity of density and pressure.

Then the set of constraints A are used to test if the candidate solution satisfies ρ i > 0 and p i > 0. Note that p i is not a conservative variable and derives from nonlinear combinations of conservative ones. The PAD algorithm is the following.

1. The Physical Admissibility criterion is first checked on each cell K i

ρ i > 0, p i > 0. ( 13
)
2. If the PAD criterion is not satisfied then d i is decremented, else U * i is eligible.

The PAD procedure only consists of maintaining the physical meaningfulness of the numerical approximation. In other words, the high-order MOOD method coupled with the PAD Detection Process is positivity-preserving for density and pressure. This point is further discussed in section 3.2.3.

Extension of the u2 detection process

Physical admissibility of the solution is not enough to prevent oscillations in the vicinity of discontinuities. It is a precondition but we require an supplementary detection criterion to decide whether the numerical solution is locally smooth or not. To this end, we adapt the u2 criterion to the density variable using local P 2 polynomial reconstruction ρ i = ρ i (.; 2) to evaluate X i = ∂ xx ρ i and Y i = ∂ yy ρ i . The u2 detection algorithm for the Euler system is thus the following.

1. The PAD criterion is first checked on each cell K i . If it is not satisfied then d i is decremented and Steps 2. and 3. are skipped. 9) and [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One dimensional systems[END_REF]. If cell i is not a smooth extrema then d i is decremented for any conservative variable, else U * i is eligible.

The DMP criterion of the density function is checked on each cell

K i min j∈ν(i) (ρ n i , ρ n j ) ≤ ρ i ≤ max j∈ν(i) (ρ n i , ρ n j ). ( 14
) 3. If ρ i does not satisfy (14) then a-Compute X k , Y k for k ∈ ν(i) {i} and coefficient ε i , b-Check criteria (
The set of constraints A consists of the PAD, and the u2 detection process on the density. Note that the density is thus the variable onto which the detection is performed. However there is a large number of possible choices of detection variables and decrementing procedures.

Positivity-preserving property

One important property a scheme must fulfill is to be positivity-preserving, that is given a set of physically admissible mean values the scheme provides another set of physically admissible ones. It is absolutely mandatory for the simulation to continue. In the case of the Euler equations density and pressure must be positive but this is not straightforwardly ensured by most of classical MUSCL or ENO/WENO schemes and most of simulation codes need a special treatment when the positivity is violated. Indeed designing a positivitypreserving scheme may be a difficult task and often leads to a more complex scheme because of the classical a priori limitation philosophy. This classical difficulty is stated by the authors in [START_REF] Zhang | Maximum-principle-satisfying and positivitypreserving high-order schemes for conservation laws: survey and new developments[END_REF] page 2754 as " It is very difficult to design a conservative high-order accurate scheme preserving the positivity". However the a posteriori treatment implies that the MOOD method is intrinsically positivity-preserving assuming the three following points:

1. The lowest order scheme is positivity-preserving, in our case it is the first-order finite volume one. 2. The positivity of density and pressure are parts of the set of constraints A. 3. The EdgePD strategy is upper-limiting see [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF] definition 9 page 4033. This implies that if the CellPD of a given cell is 0 then this cell is fully updated with the first-order scheme.

The proof that the MOOD method is positivity-preserving is analogous to the one in theorem 10 page 4033 of [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF]. In short, given a candidate solution one checks the positivity of density and pressure. If a cell is problematic that is to say density or pressure is negative then the CellPD is decremented. The next candidate solution is computed and checked again: Either this next candidate is positive or the decrementing process carries on until the CellPD is zero. In this latter case points 1. and 3. necessarily imply the positivity of the candidate solution. As this process is the same for any cell it leads to a positivity-preserving solution in a finite number of MOOD iterations.

In the numerical section we propose the Noh test case for which our implementation of the classical MUSCL scheme generates a negative pressure and fails to complete the simulation whereas the MOOD method always gives a physical meaningful solution.

Numerical tests

MOOD has been implemented into a 2D unstructured (polygonal) code which can deal with advection equation and hydrodynamics equations. The polynomial reconstruction ranges from piecewise constant up to piecewise polynomial of fifth degree. Following remark 4 one uses two decrementing sequences: P 5 -P 2 -P 0 and P 3 -P 2 -P 0 . It implies that only two precomputed matrices for the reconstruction step per cell are only stored in memory for d = d max and d = 2. The flux computation involves integrals which are approximated using Gaussian numerical integration. We use two Gaussian points on edges for P 2 and P 3 reconstructions and three for P 5 to reach the expected order of accuracy for numerical integrations. Time integration is performed with the RK3-TVD method given by system [START_REF] Buffard | Monoslope and Multislope MUSCL Methods for unstructured meshes[END_REF]. We apply the MOOD procedure detailed in section 2 to each sub-step of the RK3-TVD. The CellPD are thus reinitialized to d max at the beginning of each time sub-step. By default we use classical time step control with CFL=0.6. In the case of convergence study we use a fixed time step ∆t = ∆x r/3 to reach r th -order of accuracy. Given a variable ϕ the relative L 1 and L ∞ errors are measured by:

err 1 = i∈E el |ϕ N i -ϕ 0 i ||K i | i∈E el |ϕ 0 i ||K i | and err ∞ = max i∈E el |ϕ N i -ϕ 0 i | max i∈E el |ϕ 0 i | ,
where (ϕ 0 i ) i and (ϕ N i ) i are respectively the cell mean values at initial time t = 0 and final time t = t final = N ∆t. The unstructured meshes used in this paper are of different kinds, logically rectangular, Delaunay triangulation, Voronoi tessellation and non-conformal polygonal mesh. Contrarily to what was done in [START_REF] Clain | A high-order finite volume method for hyperbolic systems: Multi-dimensional Optimal Order Detection (MOOD)[END_REF] the whole detection is made a posteriori, namely we do not check if the reconstructed values at Gauss points are physically admissible or not. If they are not, the flux and the cell mean values are usually undefined therefore the cell is flagged as problematic.

Advection equation

Let us consider the scalar linear advection of a quantity u with velocity V (x)

     ∂ t u + ∇.(V u) = 0, u(., t = 0) = u 0 , ( 15
)
where V (x) is a continuous function on Ω ∈ R 2 and u 0 is the initial condition. Boundary conditions are prescribed as periodic ones on ∂Ω. The Double Sine Translation (DST) is first tested on Delaunay triangulations and Voronoi tessellations in order to prove that on smooth solution MOOD can actually maintain very high-order of accuracy with the u2 detection criteria. Only second-order of accuracy is reached when DMP detection criterion is used. The second test is the Solid Body Rotation (SBR) that is used to prove that MOOD-u2 can preserve smooth extrema but can still limit discontinuous profiles. This problem is further used to show the improvement obtained when polynomial reconstruction degree is increased, in other word when high-(P 1 ) and very high-order (P 3 , P 5 ) numerical schemes are used.

Double Sine Translation (DST)

Let Ω be the unit square. We consider a constant velocity V = (2, 1) and the C ∞ initial condition u 0 (x, y) = sin(2πx) sin(2πy).

The final time is t final = 2.0. Periodic boundary conditions imply that the exact final solution coincides with the initial one. The solution is therefore always smooth during the computation. The computations are carried out on series of successively refined Delaunay triangulations (from 456 up to 29184 cells, see an example in Fig. 3 left panel) and polygonal Voronoi tessellations (from 300 up to 19200 cells, see Fig. 3 right panel). Note that the meshes are far from being regular, see right panel of Fig. 3 fro instance. We plot in Fig. 4 the convergence curves obtained on the series of Delaunay triangulations and Voronoi tessellations. The MOOD method with the DMP detection process is displayed on top panels whereas the u2 Detection Process is on bottom panels. It clearly shows the strong limitation implied by the DMP since only 3 rd -order and 2 nd -order are reached in L 1 and L ∞ norms respectively independently of the polynomial degree. On the contrary the proposed u2 Detection Process reaches the expected order of convergence. This is actually explained by the fact that only polynomials of maximal degree are used during the whole computation, i.e. no CellPD decrementing is ever recorded. L 1 and L ∞ errors and rates are given in Table 1 for the DMP and the u2 detection criteria. One observes that the optimal order of convergence is reached for the u2 detection criterion whereas only second-order accurate results are obtained when the DMP is used. This accuracy test on smooth functions is passed by the MOOD method with u2 Detection Process, the next section is thus dedicated to the study of its behavior on non-smooth profiles.

Solid Body Rotation (SBR)

First introduced by R.J. Leveque in [START_REF] Leveque | High-Resolution Conservative Algorithms for Advection in Incompressible Flow[END_REF], the Solid Body Rotation test on the unit domain consists of one rotation of three shapes: a hump, a cone and a slotted cylinder. Each shape is located within a circle of radius r 0 = 0.15 Hump centered at (x 0 , y 0 ) = (0.25, 0.5)

u 0 (x, y) = 1 4
(1 + cos(π min(r(x, y), 1))).

Cone centered at (x 0 , y 0 ) = (0.5, 0.25)

u 0 (x, y) = 1 -r(x, y).
Slotted cylinder centered at (x 0 , y 0 ) = (0.5, 0.75)

u 0 (x, y) =      1 if |x -0.5| < 0.25, or y > 0.85, 0 elsewhere, where r(x, y) = 1 r 0 (x -x 0 ) 2 + (y -y 0 ) 2 .
To perform the rotation, we use the velocity V (x) = (-y + 0.5, x -0.5) and the final time t final = 2π corresponds to one full rotation. For this test case we use a genuinely unstructured and non-uniform mesh made of 5190 triangles see Fig. 5 where we also display the initial data in isolines view, see also Fig. 6 top-left panel where a side view of the initial data is provided. This mesh is refined around the slotted disk, the ratio between the largest and smallest edge length is approximately 7. The three shapes while rotating move across the refined and coarse zones. The purpose is to emphasize the effects on the numerical results of using a truly non-regular mesh. We plot in Fig. 6 profile views of the solution obtained from three methods but all with a P 5 polynomial reconstruction. First the MOOD method with the DMP Detection Process, then the MOOD method with the u2 Detection Process, and finally the unlimited version of the FV scheme. These results show on one hand that the solution with u2 Detection Process on the nonsmooth slotted cylinder is almost the same as for the DMP. On the other hand it shows that the u2 solution on the two smooth profiles are exactly the ones obtained by the unlimited scheme. In other words, the u2 Detection Process maintains the same accuracy as an unlimited scheme on smooth profiles and almost the monotonicity of a limited scheme on non-smooth ones. The same conclusion applies for any other polynomial degrees tested hence we have skipped these figures. In Fig. 7 are displayed a zoom on the slotted disk at the final time for the initial/final, the limited MUSCL scheme (MLP [START_REF] Park | Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids[END_REF]), MOOD-P 1 , MOOD-P 3 and MOOD-P 5 with u2 detection process. In Table 2 are gathered the errors for P 3 and P 5 in order to show that the u2 Detection Process provides a slightly better accuracy than the DMP detection process. Finally we display in Table 3 the min/max values of the final numer- ical solution for the limited MUSCL method (MLP), MOOD-P 1 , MOOD-P 3 and MOOD-P 5 all with DMP detection or u2 Detection Process. This table shows that the u2 detection process permits slight undershoots which is one of the reasons MOOD can reach high-order accuracy.

L 1 Error DMP u2 UNLIMITED P 3 3.219E-1 3.171E-1 3.734E-1 P 5 2.690E-1 2.621E-1 3.223E-1

Euler system

In this section we test the MOOD method on unstructured meshes for hydrodynamics problems governed by the Euler system. First we need to assess the effective numerical accuracy of the method on a smooth problem for which an exact solution exists. We choose an isentropic vortex which presents a smooth profile during the entire simulation and, as such, permits the estimation of errors and convergence orders. In a second test we run the Lax shock tube to assess the essentially non-oscillatory behavior of MOOD compared to classical WENO results. Then we run the Double Mach reflection problem to highlight the good capacity of the MOOD method to capture strong shocks and contact discontinuities. Moreover we provide CPU cost and memory storage tables.

Next the Noh problem is used to assess the positivity-preserving property of the MOOD method. Last we propose a genuine physical problem extracted from [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF] for which experimental results are available.

Isentropic vortex

The isentropic vortex problem is detailed in [START_REF] Shu | Essentially non-oscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws[END_REF] and [START_REF] Yee | Entropy Splitting and Numerical Dissipation[END_REF], therefore we only mention the basic data for the sake of consistency. The simulation domain Ω is the square [-5, 5] × [-5, 5] and we consider an initial gas flow given by the following condition (ambient gas

) ρ ∞ = 1.0, u ∞ = 1.0, v ∞ = 1.0, p ∞ = 1.0,
with a normalized ambient temperature T * ∞ = 1.0 computed with the perfect gas equation of state and γ = 1.4.

A vortex centered at x vortex = (x vortex , y vortex ) = (0, 0) is added to the ambient gas at the initial time t = 0 with the following conditions u = u ∞ + δu, v = v ∞ + δv, and

T * = T * ∞ + δT * δu = -y β 2π exp 1 -r 2 2 , δv = x β 2π exp 1 -r 2 2 , δT * = - (γ -1)β 8γπ 2 exp 1 -r 2 .
with r = x 2 + y 2 , (x = x -x vortex , y = y -y vortex ) and vortex strength is given by β = 5.0. Consequently, the initial density is given by

ρ = ρ ∞ T * T * ∞ 1 γ-1 = 1 - (γ -1)β 8γπ 2 exp 1 -r 2 1 γ-1 (16)
We assume periodic condition on the boundary and the exact solution at any time t is the same vortex but translated.

The goal of the present test is to highlight the stagnation of the rate of accuracy when primitive variables are used for the polynomial reconstructions instead of conservative ones. As pointed out in the previous section, nonlinear operations on means values reduces the method order up to at most a second-order one. We have performed the numerical simulations of the isentropic vortex problem with the same mesh using the less restrictive Physical Admissible Detection (PAD) procedure to provide effective very high-order. A series of refined meshes (from 200 up to 51200 cells) are successively used to compute the numerical solution.

In Table 4 are gathered the L 1 and L ∞ errors and rates of convergence for MOOD-P 2 , MOOD-P 3 , MOOD-P 5 using the Physical Admissible Detection Clearly, conservative variable reconstructions provide the optimal convergence rate whereas the reconstruction with primitive variables is systematically reduced to a second-order one. We also display in Fig. 8 corresponding to the errors of Table 4. Finally we also mention that when the vortex is not in motion, i.e. (u ∞ , v ∞ ) = (0, 0), then the reconstruction using primitive variables does produce the correct order of convergence.

Lax shock tube

The 1D Lax shock tube consists of two states (ρ L , u L , p L ) = (0.445, 0.698, 3.528) and (ρ R , u R , p R ) = (0.5, 0, 0.571) separated by the interface x = 0. In order to compare with the finite volume multi-dimensional WENO results of [START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF], we run the problem on the domain Ω = [-1; 1] × [0, 0.2] until final time t = 0.26 using a mesh made of 100 × 10 quadrangles split into two triangles through the same diagonal for all cells (see Fig. 5.5 of [START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF]). The goal of this test is to compare the essentially non-oscillatory behavior of the MOOD method using the u2 detection process with the classical genuinely multi-dimensional finite volume WENO results based on conservative or characteristic variables.

On the top panels of Fig. 9 we reproduce the density profiles from [START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF] obtained with the 4 th -order WENO method based on conservative (left) and characteristic (right) variables. We recall that this method uses combinations of P 2 (3 rd -order) polynomials to reach 4 th -order. The bottom panel presents the MOOD-P 3 density profiles on conservative variables with the u2+PAD detection process where only one line of triangles is displayed to comply with Fig. 5.6 of [START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF]. On the left we plot the result obtained with the Lax-Friedrichs numerical flux (which is used for the WENO results) and on the right the result using the HLLC flux. We observe that the non-oscillatory behavior of the MOOD method with u2+PAD detection process is equivalent to the WENO on characteristic variables while it is clearly better than the WENO on conservative variables. Moreover we see that the use of HLLC gives better a result for a negligible additional cost, with only three points in the contact discontinuity instead of five and remains essentially non-oscillatory.

Double Mach reflection of a strong shock

The double mach reflection of a strong shock was first proposed in [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. This test problem involves a Mach 10 shock in a perfect gas with γ = 1.4, which is initially positioned at x = 1/6, y = 0 and makes a 60 • angle with the x-axis.

The gas ahead of the shock is at rest and has uniform initial density ρ 0 = 1.4 and pressure p 0 = 1. The reflecting wall lies along the bottom of the domain, beginning at x = 1/6. The region from x = 0 to x = 1/6 along the bottom boundary at y = 0 is always assigned values for the initial post-shock flow.

Inflow boundary condition on the left side and outflow condition on the right side are also set. At the top, the boundary conditions are set to describe the exact motion of the Mach 10 flow (see [START_REF] Cockburn | The Runge-Kutta Discontinuous Galerkin Method for Conservation Laws V: Multidimensional Systems[END_REF]).

The goal of the test is, on one hand, to quantitatively show the effect of the polynomial degree reconstruction when dealing with strong shock and, on the other hand, to observe the capacity of the method to reproduce the complex structure due to the contact discontinuities in the right part of the shock. The mesh has been obtained using the free mesher Gmsh by a refinement of a coarser Delaunay ones, it is constituted of 102720 triangles (see Fig. 10 top.). Moreover for all figures 30 isolines between 1.39 and 23 have been drawn.

We depict in Fig. 10 the impact of the polynomial degree of the reconstruction on the numerical solution using the same mesh. The u2+PAD Detection Processes has been employed to control the oscillations in the vicinity of the shock. Clearly the degree of the reconstruction has a strong impact on the solution accuracy and improve the shock capture. Most relevant parts are the contact discontinuities in the right zone x ∈ [2.3, 2.7] which show the capacity of the scheme to reduce numerical viscosity when employing higher-order reconstructions.

Figure 11 is a comparison between the Physical Admissible Detection (PAD) and the coupling u2+PAD. The u2 Detection Process reduces the oscillations but increases the numerical viscosity close to contact discontinuities. It is worth noting that even with a weak Detection Process, namely the PAD procedure, the MOOD method is still very robust and provides a solution resembling the classical ones from the literature [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. The choice of the detecting procedure depends of the simulation goal: Less oscillations with the u2+PAD or less diffusive with the PAD alone.

To conclude with this test case, we provide in Tables 5 and6 the cost of the MOOD method running on a single core of the three following machines (using -O3 flag for gfortran compiler) M1 : A laptop with Intel Core2Duo P7550 (2 cores) @ 2.26GHz, 3MB of L2 Cache, 8GB of RAM. M2 : A server with two Intel Xeon E5335 (4 cores) @ 2.00Ghz, 8MB of L2 This comparison is done on two different meshes, one made of 57600 uniform quadrilaterals and one Delaunay triangulation with 17624 cells. We compare MOOD-P2, MOOD-P3 and MOOD-P5 for both the PAD and u2+PAD detection processes. We give in Table 5 the memory cost (in left column) and the total number of iterations (in right column) for all simulations, while we provide in Table 6 the total CPU time (in left column) and the time in microseconds needed for one complete time step of a single cell (in right column) including reconstruction, flux computation and time integration (RK3) of all variables. It is fairly difficult to compare the cost of two methods running on different machines, for instance the method is faster on triangles with M 2 compared to M 1 but it is the opposite for quadrilaterals. However according to reference [START_REF] Dumbser | Quadrature-Free Non-Oscillatory Finite Volume Schemes on Unstructured Meshes for Nonlinear Hyperbolic Systems[END_REF] the MOOD method is competitive when compared to truly unstructured methods of the same order.

Noh problem as a positivity-preserving test case.

The goal of the Noh problem in Cartesian geometry is to numerically prove that the MOOD method is positivity-preserving, see section 3.2.3 for a discussion on this point. It is a difficult problem well-known in the Lagrangian community, see as instance [START_REF] Maire | A high-order cell-centered Lagrangian scheme for twodimensional compressible fluid flows on unstructured meshes[END_REF][START_REF] Loubère | Staggered Lagrangian discretization based on cell-centered Riemann solver and associated hydrodynamics scheme[END_REF]. It is noticeable that our implementation of the classical MUSCL scheme is not able to simulate this problem without creating negative pressures. The problem is run in the disk of radius 1.2 centered at (0, 0). We initialize a perfect gas with γ = 5/3, density ρ 0 = 1, pressure p 0 = 10 -10 and velocity U 0 (x, y) = -x/ √ x 2 + y 2 , -y/ √ x 2 + y 2 such that ||U 0 (x, y)|| = 1. A cylindrical shock wave generated at the origin further diverges until final time t final = 2.0. The exact solution at t final is thus given by

{ρ, p, u r } =      16, 16 3 , 0 if r < r s , (1 + 2 r ), 10 -10 , -1 if r > r s , ( 17
)
where r is the radius, u r the radial velocity and r s = 2/3 the shock wave position. This problem is simulated on a polygonal mesh made of 19756 cells with about 100 cells in the radial direction. Notice that the mesh is made of seven layers of quadrangles separated with degenerated polygons, see Fig. 12.

We display the MOOD-P3 results for the density maps (left panels) and the density as a function of cell radius (right panels) in Fig. 12. The top panels correspond to the PAD detection process whereas the bottom ones correspond to the u2+PAD process. One obverves that the symmetry is almost perfectly reproduces. Notice that the PAD detection process is only intended to ensure the physical meaningfulness of the solution but does not prevent oscillations to occur. Independently of the order of the scheme the PAD always provides a meaningful solution. As a consequence the u2+PAD not only provides a valid solution without negative pressure but also removes the oscillation after the shock wave.

Impact of a shock on a cylindrical cavity

We finally test the ability of MOOD method to capture physics in realistic conditions by simulating the experiment proposed in [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF] where a planar shock impacts a cylindrical cavity. We consider the case of a nominal incident shock Mach number of 1.33 in ambient air (with γ = 1.4) at 0.95 bar pressure. Moreover we use the domain configuration A (following notation of [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF]) we detail in Figure 13. The variables initialization is split in two parts, the pre-shock values (ρ, u, v, p) = (1.1175, 0.0, 0.0, 95000.0), and the post-shock ones (ρ, u, v, p) = (1.7522, 166.3435, 0.0, 180219.75), leading to conditions of [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF] at temperature T = 296.15K.

The simulation is only preformed in the lower half part of the domain for symmetry argument, namely from y = 0 mm to y = 75 mm. The 193615 cells mesh is composed of triangles, quadrangles but also more general polygons with non-conformal elements (see Figure 14) to better suit with the complex geometry of the set-up. Notice that non-conformity is simply handled using polygons, i.e. no special treatment is used. We also deliberately use a heterogeneous mesh to highlight that the MOOD method is not much affected by the quality of the mesh. The simulation are carried out with the MOOD-P3 method (fourth-order) using the PAD and the u2 Detection Process. Pictures are rendered as a full mesh by symmetry even if the computation was done on a half-domain to easier compare with physical results of [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF]. Figure 15 represents the density gradient magnitude at six different times to embrace the global behavior of the solution. In details, top right is chosen to be compared to Fig. 7 (a) of [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF], bottom center to Fig. 8 (d) and bottom right to Fig. 9 (c) of same paper. Our results are clearly in agreement with physical results. In Figure 16 different zooms on solution at several times are plotted. On the top part, density gradient magnitude at a late time is given and is to be compared to Fig. 14 (b) in [START_REF] Skews | Flow features resulting from shock wave impact on a cylindrical cavity[END_REF] while we superpose, in the bottom figure, the velocity vectors on the density magnitude gradient to show the created vortices at the entrance of the cavity (left) and highlight the instabilities lying along the wall.

Conclusion and perspectives

The paper presents important new extensions of the MOOD method for unsteady advection and hydrodynamics equations, that ensure high-order approximations (up to the sixth-order) on unstructured meshes. We introduced new efficient detection processes and proved that the MOOD method is intrinsically positivity-preserving for the hydrodynamics system of equations assuming that the first-order scheme is. This has been numerically assessed on the Noh problem for which our implementation of the MUSCL scheme fails due to negative pressures. Then both for the advection equation and the hydrodynamics Euler system, we proposed numerical tests to confirm that the MOOD method provides very high-order of accuracy on unstructured meshes for smooth solutions (e.g. isentropic vortex in motion) and non-oscillatory behaviour on discontinuous ones (e.g. Lax shock tube). Moreover the memory storage and CPU time have also been reported for the double Mach problem, proving that the MOOD method is competitive. The last numerical test showed that the MOOD method, on a relatively coarse and non-conformal polygonal mesh, is able to simulate complex physics from an experimental set-up of the impact of a shock wave on a cylindrical cavity.

Finally we plan to improve the detection procedure, especially for vectorial problems to achieve a very low diffusion but still preventing the oscillations from appearing. Application to full three-dimensional problem is also an attractive task since performing an efficient computational solution is always a challenging problem. The extension of the MOOD method to deal with steady-state solution needs also more investigations. Overall an important perspective is the polynomial reconstruction itself. We have observed that the main computational cost comes from the reconstruction stage and that the reconstruction quality strongly depends on the stencil employed. Such a point is of crucial importance from a computational point of view to obtain tractable complex numerical simulations. Appendix: The Discrete Maximum Principle on mean values provides at most a second-order scheme.

We recall that a time explicit scheme preserves the Discrete Maximum Principle (DMP) if for all cell K i min j∈ν(i)

(U n i , U n j ) ≤ U n+1 i ≤ max j∈ν(i) (U n i , U n j ). ( 18
)
It has been shown in [START_REF] Sander | A third-order accurate variation non-expansive difference scheme for single nonlinear conservation law[END_REF][START_REF] Osher | High resolution schemes and the entropy condition[END_REF][START_REF] Jiang | Non-oscillatory central schemes for multidimensional hyperbolic conservative laws[END_REF] that any scheme based on the DMP property reduces the accuracy to second-order for regular functions due to inaccurate approximation at extrema. Indeed following [START_REF] Zhang | Maximum-principle-satisfying and positivitypreserving high-order schemes for conservation laws: survey and new developments[END_REF], let us consider the advection problem in R to avoid boundary condition issues

     ∂ t U + ∂ x U = 0, U (x, t = 0) = cos(x). ( 19
)
We consider a uniform discretization x i = ih, i ∈ Z and h > 0 being the cell size and initialize the mean value on cell K 0 = [0, h] as

U t=0 0 = 1 h h 0 cos(x) dx = sin(h) h . ( 20
)
Now, let us perform one time step with ∆t = h/2 of a finite volume scheme which respects the DMP property. The exact solution at time t = h/2 is U ex (x, h/2) = cos(x -h/2) and accordingly the exact mean value on K 0 is

U ex,t
U t=0 0 = sin(h) h = 1 - h 2 6 + O(h 4 ).
Clearly, the exact mean value U ex,t 0 on cell K 0 is greater than the maximum mean values over all cells at time t = 0 with an error of h 2 /8 as

|U ex,t 0 -U t=0 0 | ≤ | h 2 24 - h 2 6 + O(h 4 )| = h 2 8 + O(h 4 ).
Therefore a scheme which fulfills the DMP property necessarily provides a solution lower than sin(h)/h, hence after the first cycle the numerical solution verifies U t 0 ≤ U t=0 0 = 1 -h 2 6 + O(h 4 ). It follows that the approximation of the mean value has an error of order O(h 2 ) compared to the exact mean value on cell K 0 . Consequently the scheme is at most second-order accurate and DMPtype of criteria cannot be used strictly for higher than second-order schemes and has to be relaxed.

Fig. 2 .

 2 Fig. 2. Mesh notation. Index set ν(i) corresponds to blue cells with dots, ν(i) corresponds to non-white cells.

Remark 1

 1 In 2D, at least N (d) = (d + 1)(d + 2)/2 -1 neighbors are needed to provide the minimal number of equations. However for the sake of robustness more cells are involved. In details, we use at least 5 cells for d = 1, 8 cells for d = 2, 16 cells for d = 3, 20 cells for d = 4 and 26 cells for d = 5.

Fig. 3 .

 3 Fig. 3. Example of Delaunay (left) and Voronoi (right) meshes for the DST problem.

Fig. 4 .

 4 Fig. 4. Error curves for the DST problem for series of Delaunay meshes (empty symbols) and of Voronoi meshes (filled symbols) for the DMP detection process (top) and the u2 one (bottom).

20

 20

Fig. 5 .

 5 Fig. 5. Initial mesh and initial data for the SBR problem. The mesh is composed of 5190 triangles refined around the slotted disk. The resulting mesh is genuinely non-uniform.

5 Fig. 6 .

 56 Fig. 6. Profiles of the SBR solution for the initial/final exact solution, MOOD-P 5 with DMP detection process, MOOD-P 5 with u2 detection process, MOOD-P 5 without any limitation.

Fig. 7 .

 7 Fig. 7. Profiles of the SBR solution for the initial/final exact solution (top), for a limited MUSCL method (MLP) and MOOD-P 1 , MOOD-P 3 , MOOD-P 5 with u2 detection process.

Fig. 8 .

 8 Fig. 8. Convergence curves for the isentropic vortex. Top figures correspond to the reconstruction with primitive variables while bottom figures use reconstruction with conservative variables. The left column represents the L 1 -norm error and the right column the L ∞ -norm error. The PAD detection process has been used.

Fig. 9 .

 9 Fig.9. Comparison between WENO and MOOD methods on 100 × 10 quadrangles split into triangles -Top: results of the 4 th -order WENO method using Lax-Friedrichs flux from[START_REF] Hu | Weighted essentially non-oscillatory schemes on triangular meshes[END_REF] on conservative variables (left) and on characteristic variables (right) -Bottom: results of the 4 th -order MOOD-P 3 method with u2+PAD detection process using Lax-Friedrichs (left) and HLLC (right) fluxes.

Fig. 10 .

 10 Fig. 10. Comparison between the P 2 (top), P 3 (middle) and P 5 (bottom) polynomial reconstructions with the conservative variables using the same mesh. Physical Admissible Detection (PAD) and u2 Detector have been both used to prevent numerical oscillations

Fig. 11 .

 11 Fig. 11. Results of the MOOD method with P5. In the top figure, simulation has been carried out with the Physical Admissible Detection (PAD) Detection Process while we have both employed the PAD and u2 Detection in the middle figure. The left bottom and right bottom figure give a zoom of the solution with the PAD and u2+PAD Detection Process respectively.

Fig. 12 .

 12 Fig. 12. Noh problem at t final = 2.0 on a polygonal grid -Left: Density map and mesh -Right: Cell density as a function of cell radius vs exact solution -Top panels correspond to the PAD detection process -Bottom panels correspond to the u2+PAD detection process.

Fig. 13 .

 13 Fig. 13. Domain characteristics for the shock impacting a cylindrical cavity. Red arrows represent inflow and outflow boundary conditions.

Fig. 14 .

 14 Fig. 14. On top, we display the global view of the mesh where the different mesh zones are clearly visible. On bottom, zooms on the non-convex part of the mesh (on left) and on the junction between the polar part of the mesh and the quasi-uniform one (right). Non-conformity are clearly visible.

Fig. 15 .

 15 Fig. 15. Gradient density magnitude is shown at different times. Time 0 corresponds to the initial shock at position x=0.

Fig. 16 .

 16 Fig. 16. Zooms on different parts of the solution. On top, gradient density magnitude is shown at a late time when instabilities are well developed. On bottom, vortices at the entry of the cavity (left) and the instabilities (right) along the wall are displayed with density gradient magnitude in color and velocity vectors.

1 .

 1 Initialization at t n . The MOOD procedure starts by initializing the CellPD to d i = d max and by computing the coefficients of the polynomial reconstruction U i (x; d i ) on each cell. 2. Evaluation of EdgePD and values at Gauss points. We compute the EdgePD d ij on each edge and use polynomial function U Computation of candidate solution U h . Numerical fluxes are computed using the reconstructed solution at Gauss points and one time step is carried out to provide a candidate U h at time t n+1 = t n + ∆t n . 4. Check U i for A-eligibility. If d i = 0 we check the A-eligibility of each mean value U i with respect to the constraints set A. In the case U i is not A-eligible then CellPD d i is decremented. If all cells are A-eligible then the candidate solution is valid and we set U n+1 h

i (x; d ij ) and U j (x; d ij ) to compute approximations of U at Gaussian points on e ij .

3.

Table 1 L

 1 ).1 and L ∞ errors and convergence rate for the DST problem for the MOOD method with DMP and u2 detection process.

			DMP detec. process		u2 detec. process
	Deg./Type	Cell Nb	L 1 error		L ∞ error		L 1 error		L ∞ error
	P 2 /Delaunay	456	1.775E-01	-	2.629E-01	-	1.656E-01	-	1.549E-01	-
			2.303E-02 2.95 8.016E-02 1.71	2.351E-02 2.82 2.283E-02 2.76
			3.142E-03 2.87 2.522E-02 1.67	3.049E-03 2.95 2.995E-03 2.93
		29184	4.391E-04 2.84 8.082E-03 1.64	3.870E-04 2.98 3.784E-04 2.98
	P 2 /Voronoi	300	4.804E-01	-	5.278E-01	-	4.423E-01	-	4.339E-01	-
			7.483E-02 2.68 1.359E-01 1.96	7.482E-02 2.56 7.070E-02 2.62
			9.779E-03 2.94 3.432E-02 1.99	9.788E-03 2.93 9.348E-03 2.92
		19200	1.244E-03 2.97 1.039E-02 1.72	1.233E-03 2.99 1.176E-03 2.99
	Expected order			3		3		3		3
	P 3 /Delaunay	456	6.383E-02	-	1.801E-01	-	9.474E-03	-	1.007E-02	-
			8.369E-03 2.93 5.920E-02 1.61	5.751E-04 4.04 7.916E-04 3.67
			9.916E-04 3.08 2.057E-02 1.53	3.611E-05 3.99 4.664E-05 4.09
		29184	1.185E-04 3.06 7.146E-03 1.53	2.140E-06 4.08 2.774E-06 4.07
	P 3 /Voronoi	300	1.158E-01	-	2.826E-01	-	6.431E-02	-	5.961E-02	-
			2.263E-02 2.36 9.234E-02 1.61	4.017E-03 4.00 3.632E-03 4.04
			2.157E-03 3.39 2.787E-02 1.73	2.583E-04 3.96 2.539E-04 3.84
		19200	2.393E-04 3.17 9.295E-03 1.58	1.649E-05 3.97 1.718E-05 3.89
	Expected order			4		4		4		4
	P 5 /Delaunay	456	6.098E-02	-	1.691E-01	-	3.034E-04	-	3.715E-04	-
			9.660E-03 2.66 6.383E-02 1.41	6.796E-06 5.48 9.939E-06 5.22
			1.359E-03 2.83 2.399E-02 1.41	1.207E-07 5.82 1.831E-07 5.76
		29184	1.704E-04 3.00 8.574E-03 1.48	1.767E-09 6.09 2.836E-09 6.01
	P 5 /Voronoi	300	1.352E-01	-	2.610E-01	-	4.584E-03	-	4.955E-03	-
			2.213E-02 2.61 9.116E-02 1.52	7.327E-05 5.97 8.740E-05 5.83
			2.119E-03 3.38 2.914E-02 1.65	1.341E-06 5.77 1.573E-06 5.80
		19200	2.449E-04 3.11 1.005E-02 1.54	3.017E-08 5.47 3.703E-08 5.41
	Expected order			6		6		6		6

Table 2 L

 2 1 error for the SBR problem for different detection processes and polynomial degrees.

Table 3

 3 Minimal and maximal mean values for the SBR problem for different detection processes and polynomial degrees.

  the convergence curves

			Conservative variables	Primitive variables
	Deg. Cell Nb	L 1 error		L ∞ error		L 1 error		L ∞ error
	P 2	200	1.850E-02	-	2.680E-01	-	2.002E-02	-	2.884E-01	-
		800	6.519E-03 1.50 1.255E-01 1.09 7.621E-03 1.39 1.771E-01 0.70
		3200	1.444E-03 2.17 2.208E-02 2.51 1.536E-03 2.31 4.054E-02 2.13
		12800	2.504E-04 2.53 3.631E-03 2.60 2.554E-04 2.59 6.060E-03 2.74
		51200	3.347E-05 2.90 4.923E-04 2.88 4.540E-05 2.49 8.756E-04 2.79
	Expected order		3		3		3		3
	P 3	200	1.137E-02	-	1.880E-01	-	1.424E-02	-	2.384E-01	-
		800	2.504E-03 2.18 4.686E-02 2.00 3.530E-03 2.01 8.358E-02 1.51
		3200	3.524E-04 2.83 5.977E-03 2.97 5.666E-04 2.64 8.835E-03 3.24
		12800	1.947E-05 4.18 3.725E-04 4.00 1.377E-04 2.04 1.649E-03 2.42
		51200	1.069E-06 4.19 1.996E-05 4.22 3.460E-05 1.99 4.091E-04 2.01
	Expected order		4		4		4		4
	P 5	200	8.193E-03	-	1.200E-01	-	1.161E-02	-	1.915E-01	-
		800	1.762E-03 2.22 3.433E-02 1.81 2.492E-03 2.22 3.740E-02 2.36
		3200	6.767E-05 4.70 1.133E-03 4.92 5.482E-04 2.18 6.112E-03 2.61
		12800	1.011E-06 6.06 2.237E-05 5.66 1.382E-04 1.99 1.598E-03 1.94
		51200	2.583E-08 5.29 4.809E-07 5.54 3.462E-05 2.00 4.039E-04 1.98
	Expected order		6		6		6		6

Table 4 L

 4 1 and L ∞ errors and convergence rates for the isentropic vortex problem with MOOD and the Physical Admissible Detection Process. Comparison between conservative and primitive variables polynomial reconstructions for different polynomial degrees.

Table 5

 5 MeshP 2 -3 rd -order P 3 -4 th -order P 5 -6 th -order Memory storage and total number of iterations for the Double Mach problem according to the different configurations with the MOOD method.

	Detection

Table 6

 6 Total time and cost for one complete time step of a single cell for the Double Mach problem.


Adrien Kaladji

Antoine Lucas

Alain Cardon

Pascal Haigron

Computer-aided surgery: concepts and applications in vascular surgery

Keywords: Aneurysm, radiography, surgery, Aortography, Blood Vessel Prosthesis, Blood Vessel Prosthesis Implantation, instrumentation, methods, Endovascular Procedures, instrumentation, methods, Equipment Design, Humans, Predictive Value of Tests, Prosthesis Design, Surgery, Computer-Assisted, instrumentation, Treatment Outcome endovascular aneurysm repair (EVAR), computer-aided surgery, imaging

Computer-aided surgery makes use of a variety of technologies and information sources. The challenge over the last ten years has been to apply these methods to tissues that deform, as do vessels when relatively rigid flexible objects are introduced into them (Lunderquist rigid guidewire, aortic prosthesis, etc.) Three stages of computer-aided endovascular surgery are examined: sizing, planning and intra-operative assistance. Our work shows that an approach based on optimized use of the imaging data acquired during the various observation phases (pre and intra-operative), involving only lightweight computer equipment that is relatively transparent for the user, makes it possible to provide useful (i.e. necessary and sufficient) information at the appropriate moment, in order to aid decision-making and enhance the security of endovascular procedures.

INTRODUCTION

Computer-aided surgery (Fig. 1) is a discipline based on a variety of technologies and information sources 1,2 : Medical imaging, pre-[ ] or intra-operative (CT scan, MRI, ultrasonography, coelioscopic imaging) 3 ; micro-technologies, used both for action (micro-instruments) [ ] and observation (micro-sensors) 4 ; virtual reality 5 and augmented reality 6 , in particular interactive simulation and the possibility of

[ ] [ ] [ ]
enriching reality, as perceived by the user, with data or pre-processed models; medical robotics 7 used in master-slave mode and which is

[ ] now being applied to new fields, in particular that of robot-controlled catheterization 8 .

[ ]

Computer-assisted surgery tools were first developed for use in orthopaedics and neuro-surgery. These were applications using rigid tools in the presence of tissues that were also, in the first place, rigid or supposedly so (organs encapsulated in bony structures). The challenge over the last ten years has been to apply these methods to tissues that deform, as do vessels when relatively rigid flexible objects are introduced into them (Lunderquist rigid guidewire, aortic prosthesis, etc.).

Progress in medical imaging is often seen as synonymous with high capital investment in equipment derived from technological research. We present here three applications in regular clinical use that are essentially the result of research into signal and image processing and algorithm development, requiring relatively little investment in hardware (CT scan, interventional MRI, sophisticated robotics). Three stages of computer-aided endovascular surgery are examined: sizing, planning and intra-operative assistance.

Sizing

Pre-operative sizing, particularly in tricky locations (supra-aortic trunks), and whether for an aortic stent or a stenotic lesion, necessitates a perfect description of the lesion and measurements that will lead to the choice of the most appropriate device.

In practice the surgeon is faced with the use of viewers that are generally not well suited to this type of task. He calls the help of a radiologist or uses a lightweight imaging console that represents a significant financial investment. He may also use non-specific devices made available free of charge, all of which are equivalent to a slimmed-down, multi-purpose radiology console. Sometimes the vascular surgeon may delegate sizing to product specialists working for the manufacturer of the device.

The aim of our works was to develop a tool, running on a simple portable computer, enabling the surgeon to perform sizing himself by automating the thresholding, centre-line detection and sizing functions 9 . The software should not be time-consuming, enabling sizing to [ ] be carried out with seven clicks and without a previous experience of imaging workstations for the user. It is a decision-making tool that is action-oriented (placement of a thoracic or abdominal endoprosthesis or an endovalve, or stenting of arterial stenosis).

The software developed by THERENVA (Fig. 2) provides feasibility alarms (diameter of iliac arteries, sinuosity, angle of neck). This makes it possible to optimize ordering of available devices through a link to the manufacturer s catalogue. This tool was validated by an analysis of thirty-two patients and comparison with the measurements made by an endovascular radiologist using a General Electric Advanced Vessel Analysis system; in this study there was no significant difference between the sizing determined by the radiologist and that obtained by the surgeon on EndoSize in an average time of fifteen minutes (10 to 40 min. range) 10 .

[ ] is also fully automated (Fig. 3) and concerns three essential functions: Planning Pre-operative choice of operational viewing incidences, so as to limit injections of contrast medium and radiation doses in intra-operative fluoroscopy 11 . Depending on the existence of partnerships with the suppliers of operating-room imaging systems, this [ ] function can be associated with automatic steering of a motor-driven C-arm.

Preparation for endovascular navigation in lesions to which access is complex, particularly at carotid and supra-aortic level. This planning can be performed visually using virtual angioscopy 12 , i.e. virtual exploratory navigation through pre-operative imaging, and/or [ ] be based on explicit study of interactions between the endovascular material and tissue 13,14 . This implies that the mechanical behaviour

[ ]
of the material used has been modelled, that full CT images of the aortic arch and the supra-aortic vessels have been acquired

pre-operatively, and that the material has been computer-tested to predict its behaviour (length, torque, whip). Industrial development of this application is limited first by the need to obtain the mechanical parameters of the guidewires and catheters on the market (industrially sensitive data that manufacturers are legitimately reluctant to reveal), and also by the considerable workload involved in creating material data-bases for a somewhat limited number of users.

Training, for which the application described above turns out to be an excellent tool. Surgeons are of course familiar with catheterization simulators, but generally speaking the latter do not allow quick and easy testing of endovascular tools with the data of the actual patients to be treated. Catheterization simulation with integration of the pre-operative CT-scan of a planned patient would make it possible not only to optimize device selection but also to test the material and train oneself quickly to use the new tools that are becoming available on the market with ever-greater frequency.

Intra-operative aid

It is paradoxical to note that very precise information is available pre-operatively thanks to CT-scan, with visualization of complete lesions, aneurysm or stenosis, detection of the presence of thrombus or calcification, and measurement capability, while intra-operatively, only subtraction images are available, with a notable loss of parietal information with the current c-arms. This intra-operative 2D imaging, in the absence of rotational capability, exposes the surgeon to errors (most frequently parallax error), faulty interpretation of sinuous vessels, and leads to a search for the best angle of incidence involving repeated injections of contrast medium and additional doses of radiation.

Intra-operative imaging thus provides relatively little information compared to pre-or post-operative imaging. A certain number of complications endoleaks, malpositioning, covering of renal or internal iliac arteries result from sizing errors, failure to anticipate The device consists of an assistance console (figure 4) which, after an automated registration procedure, enables the surgeon not only to perform the intervention in the usual way using the scope screen but also to monitor the progression of the inserted device through the 3D-space derived from the pre-operative CT-scan. This registration and image implantation correspond perfectly to the definition of augmented reality. Everything takes place in real time, allowing surgeons to benefit from a genuine instrument panel on which they can see: progression under fluoroscopy, progression through the large-scale 3D reconstruction of the aneurysm or the supra-aortic vessels and CT-scan slices corresponding to the position of the device. Information concerning the vessel wall is then visible: thrombus, calcification.

This tool thus enables full localization to be performed using a single angle of incidence, and it is backed up by an alarm system based on the pre-operative sizing that alerts the operator in the event of deviation from the pre-established operating plan. Clinical evaluation is currently under way on ten aortic endoprostheses and four carotid angioplasties; it shows that device placement (endoprosthesis or stent) is optimized, with greater security; there is less use of fluoroscopy to look for different angles of view; and less contrast medium is injected.

The remaining technical problem concerns deformations. Introduction of an extrastiff guidewire or an aortic stentgraft naturally makes arteries, and in particular the iliac arteries, straighten and deform. Flexion/extension movements of the neck during carotid angioplasty can also modify the tortuosity of the common or internal carotid arteries. Pre-operative planning based on the pre-operative CT-scan may thus turn out to not match the intra-operative morphology 17 . To solve this problem, the aortoiliac deformations during EVAR have been

CONCLUSION

Computer-aided vascular surgery is a newly-emerging practice. Our work shows that an approach based on optimized use of the imaging data acquired during the various observation phases (pre-and intra-operative), involving only lightweight computer equipment that is relatively transparent for the user, makes it possible to provide useful (i.e. necessary and sufficient) information at the appropriate moment in order to aid decision-making. These software solutions can be integrated into a standard operating room, and all the more easily into a hybrid high-tech installation.

In the long term, the addition of data derived from 2D or 3D intra-operative ultrasound, an intra-corporeal magnetic localizer or interventional MRI, coupled with adequate pre-operative sizing and planning taking into account possible future deformations should enable us to reach our intended goal, i.e. endovascular surgery with limited use of X-rays and with an already-available angio-navigation station at the heart of this merging of data.

Acknowledgements:

This article is based on research that was partially funded the French national research agency (ANR) through the Tecsan program (project ANGIOVISION n ANR-09-TECS-003)

  faulty intra-operative positioning. The assistance tool 15,16 developed by the INSERM team and THERENVA was tested experimentally by Rennes CHU and was [ ] awarded first prize in the AGBM 2008 Innovative Medical Technologies competition.

  element simulation. This model is based on patient specific geometry from the preoperative CT-scan and is still assessed on a large cohort to represent various anatomies. According to this model imaging fusion, including the deformations of the aortoiliac structure, between preoperative CT and intraoperative fluoroscopy with a standard c-arm could be feasible.

  Fig. 1 Computer Aided Surgery (CAS)

  Fig. 4EndoNaut augmented reality device

  S , Cinquin P , Troccaz J . Editor: Roux C , Coatrieux JL . Computer Integrated Surgery and Therapy: State of the Art . Contemporary Perspectives in Three-Dimensional Biomedical Imaging . chapter 10 : IOS Press ; Amsterdam 1997 ; 239 -310 2 . Galloway RL . The process and development of image-guided procedures . Annual Review of Biomedical Engineering . 2001 ; 3 : (1 ) 83 -108 3 . Peters TM . Image Guidance for Surgical Procedures . Physics in Medicine and Biology . 2006 ; 51 : (14 ) 505 -540 4 . Haigron P , Luo L , Coatrieux JL . Issues in image-guided therapy . IEEE Eng Med Biol Mag . 2009 ; 28 : (4 ) 96 -8 5 . Fuchs P , Moreau G . Le Trait de la r alit virtuelle Pelerin JL , Kulik C , Goksu C , Coatrieux JL , Rochette M . Fluid-structure interaction applied to the simulation of Abdominal Aortic Aneurysms . Proceedings International IEEE Engineering In Medicine and Biology Conference New York August 2006 15 . Haigron P , Lucas A . Method for assisting and guiding the navigation of a tool in anatomical structures . Brevet FR 02/04898 . Universit de Rennes Goksu C , Haigron P , Acosta O , Lucas A . Endovascular navigation based on real/virtual environments cooperation for computer assisted TEAM procedures . Proceedings of SPIE Medical Imaging 2004:Visualization, Image-Guided Procedures, and Display 2004 San Diego, USA SPIE ; Bellingham, WA 2004 ; 257 -6617 . Carrell TW , Modarai B , Brown JR , Penney GP . Feasibility and limitations of an automated 2D-3D rigid image registration system for complex endovascular aortic procedures . J Endovasc Ther . 2010 ; Aug 17 : (4 ) 527 -33

	°References:
	1 . Lavall e
	é	é é	. Les Presses de l Ecole des Mines de Paris, deuxi me dition ' è é	. 2003 ;
	6 . Eggers G , Salb T , Hoppe H , Kahrs L , Ghanai S , Sudra G , Raczkowsky J , Dillmann R , Worn H , Hassfeld S , Marmulla R . Intraoperative augmented reality : the
	surgeons view . Stud Health Technol Inform . 2005 ; 111 : 123 -5
	7 . Taylor R , Stoianovici D . Medical robotics in computer-integrated surgery . IEEE Trans Robotics and Automation . 19 : (5 ) 2003 ; 765 -781
	8 . Servatius H , Steven D , Rostock T , Hoffman B , Drewitz I , M llerleile ü	K , Aydin M , Meinertz T , Willems S . Robotic Navigation Minimizes the X-ray Exposure in
	Pulmonary Vein Isolation: A Randomized Comparison with the Conventional Approach . Heart Rhythm . 2009 ; 6 : (5 ) S341 -
	9 . Acosta O , G ksu
	Ann Vasc Surg . 2010 ; Oct 24 : (7 ) 912 -20
	11 . Goksu G , Haigron P , Zhang H , Soulas T , Le Certen G , Lucas A . 3D intraoperative localization for endovascular navigation guidance . Surgetica 2002, Computer-aided '
	medical interventions: tools and applications . Sauramps Medical ; 2002 ; 323 -9
	12 . Haigron P , Bellemare ME , Acosta O , Goksu C , Kulik C , Rioual K , Lucas A . Depth-map-based scene analysis for active navigation in virtual angioscopy . IEEE
	Transactions on Medical Imaging . 2004 ; 23 : (11 ) 1380 -90
	13 . Guilloux V , Haigron P , G ksu

é ö G , Lucas A , Kulik C , Rolland Y , Haigron P . Geometrical description of anatomical structures based on virtual exploratory navigation . Surgetica . 2005 ; Sauramps Medical ; 2005 : 437 -45 10 . Kaladji A , Lucas A , Kervio G , Haigron P , Cardon A . Sizing for endovascular aneurysm repair: clinical evaluation of a new automated three-dimensional software . ö C , Kulik C , Lucas A . Simulation of guide-wire navigation in complex vascular structures . Proceedings of SPIE Medical Imaging 2006:Visualization, Image-Guided Procedures, and Display 2006 San Diego, USA SPIE ; Bellingham, WA 2006 ; 6141 : (07 ) 1 -11 14 . é ; 1 : 2002 ; 16 .