Abderrhramane Boularas
email: boularas@laplace.univ-tlse.fr

Stéphane Clain
email: clain@math.uminho.pt

Fulbert Baudoin
email: fulbert.baudoin@laplace.univ-tlse.fr

A sixth-order finite

Keywords: Finite Volume, high-order, convection-diffusion, polynomial reconstruction, Poisson equation, curved boundary

come   L'archive ouverte pluridisciplinaire

Introduction

Very high-order finite volume method (higher than the second-order) for elliptic and parabolic operators is a recent trend and has received considerable attention during the last decade. The coupling of the Euler system with a viscous term [START_REF] Ollivier-Gooch | High-Order ENO schemes for unstructured Meshes based on least-squares reconstruction[END_REF][START_REF] Barth | Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction[END_REF][START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Ivan | High-order solution-adaptative central essentially non-oscillatory (CENO) method for viscous flows[END_REF][START_REF] Toro | ADER finite volume schemes for nonlinear reaction-diffusion equations[END_REF][START_REF] Michalak | Gooch Unstructured High-Order Accurate Finite[END_REF] or the shallow-water system with turbulence [START_REF] Cea | Vázquez-Cendón Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts[END_REF] are, among others, strong motivations to design efficient and accurate schemes for elliptic operators in the finite volume context. There exists a large literature on secondorder approximations and convergence analysis [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection diffusion problem[END_REF][START_REF] Coudière | Convergence rate of a finite volume scheme for the linear convectiondiffusion equation on locally refined meshes[END_REF][START_REF] Eymard | The finite volume method, Handbook for Numerical Analysis[END_REF][START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF][START_REF] Hernández | High-order finite volume schemes for the advection-diffusion equation[END_REF][START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF][START_REF] Manzini | A finite volume method for advection-diffusion problems in convectiondominated regimes[END_REF][START_REF] Coudière | The discrete duality finite volume method for convection-diffusion problems[END_REF] but few studies have been done for very high-order [START_REF] Hidalgo | ADER Schemes for Nonlinear Systems of Stiff Advection-Diffusion-Reaction Equations[END_REF][START_REF] Zambra | Finite Volume Schemes of High Order of Accuracy in Time and Space for Flows in Unsaturated Porous Media[END_REF]. In a recent paper [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF], a new high-order finite volume method for convectiondiffusion problem has been developed for two-dimensional geometries providing up to sixthorder approximations. The technique is based on specific polynomial reconstructions to evaluate the fluxes across the cell interfaces. Nevertheless, the original technique does no take into account curved boundaries (non-polygonal domains) and the method accuracy is lowered to a second-order one due to the approximation of arcs by segments. The goal of the present study is, on the one hand, to achieve an extension for the three-dimensional case and, on the other hand, to develop a new class of polynomial representation for the boundary to preserve an effective sixth-order approximation even with curved boundaries.

When dealing with second-order finite volume schemes for elliptic operator, the curved boundary is substituted by a polygon or a polyhedra and the Dirichlet condition is evaluated at the vertices or the midpoint of the boundary edges. Very few papers in the finite volume community tackles the question of curved boundary since the problem is relevant for at least third-order schemes. In [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF], the authors propose a polynomial reconstruction which exactly matches the Dirichlet condition at the Gauss points on the curved boundary. The procedure then involves two new linear constraints added to the matrix system we use to compute the polynomial coefficients. Such an algorithm is extended to the Neumann condition on curved boundary. Another approach consists to apply a local mapping px, yq Ñ pξ, ηq to transform the cell into a new computational domain which matches the boundary [START_REF] Wang | Extension of the spectral volume method to high-order boundary representation[END_REF]. Integrals therefore involves the Jacobian matrix of the transformation and polynomial reconstruction procedure is performed in the local basis leading to a rather complex algorithm. The method we propose in the present study is related to the works of [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF] but the Dirichlet condition is enforced in a different way by using the mean value on the edge as a free parameter to fit with the boundary condition at the Gauss points. On the other hand, very high-order method leads to an important computational effort due to the reconstruction process and a large part of the time consuming derives from the polynomial coefficients evaluation. We proposed a new procedure to dramatically reduce the computational cost by a pre-evaluation of local vector for each Gauss point one can identify as a partial assembly procedure. It results that the polynomial evaluation at the Gauss points is reduced to a simple inner product between the local vector and the data saving a lot of memory and time.

In the present document, we do not deal with the convective part on purpose since the main difficulty concerns the diffusive contribution so we only focus on the Poisson problem. The second section is devoted to the generic very high-order finite volume scheme while the third section focus on the polynomial reconstruction and detail all the improvements we propose.

Numerical experiments are presented in the four section to assess the scheme accuracy and robustness as well as the efficiency of the iterative solver coupling with a new preconditioning strategy.

Generic high-order Finite Volume Scheme

We consider an open bounded domain Ω of R 2 or R 3 with a piecewise regular curved boundary. Γ D and Γ N define a partition of the boundary BΩ where we shall prescribe the Dirichlet and the Neumann conditions respectively. We intend to compute accurate approximations of function V solution of the Poisson equation

´∇.pε∇V q " g, in Ω, ε∇V.n " 0, on Γ N , V " V D on Γ D , (1)
with ε a positive function which may present some discontinuities, V D is a given function defined on Γ D and g is the source term. We assume that boundary Γ N is composed of lines (2D case) or planes (3D case) while Γ D is curved. In some applications, we shall split the domain into two sub-domains Ω 1 and Ω 2 shared by an interface Γ and we denote by ε ℓ " ε |Ω ℓ , ℓ " 1, 2 the restrictions on each sub-domain. Functions ε 1 and ε 2 are regular but ε presents a discontinuity at the interface Γ. Such an assumption is required when dealing with two different materials (for instance the air and a dielectric material).

We introduce the notations to derive the finite volume scheme (see Fig. 1). For the twodimensional case, the domain is divided into non-overlapping convex polygonal cells c i , i " 1, . . . , I and we denote by e ij " c i X c j the common interface shared by two adjacent cells while n ij stands for the unit normal vector from c i towards c j . To handle the edges on the boundary, we introduce the notation e iD which corresponds to an edge of cell c i which belongs to Γ D while n iD is the unit outward normal vector. In the same way we define e iN and n iN for the edge associated to Γ N . At last ν i is the index subset of the adjacent cells including symbols D, N when dealing with a boundary edge. To perform accurate numerical integration, we denote by q ij,r , r " 1, . . . , R 2 the Gauss points on edge e ij while ξ r , r " 1, . . . , R 2 represent the respective weights. We adopt a similar notation setting q iD,r (resp. q iN,r ) for the Gauss points on edge e iD (resp, e iN ).

The three-dimensional situation is more complicated due to the existence of possible noncoplanar interfaces. For the sake of simplicity, we shall only consider tetrahedral cells c i , i " 1, . . . , I with triangular faces f ij " c i X c j with n ij the unit normal vector from c i pointing to c j . As in the two-dimensional case, we use symbol D (resp. N ) to represent a boundary face f iD (resp. f iN ) and its associated outward normal vector n iD (resp. n iN ). We also introduce the set ν i to represent the index of the neighbour cells which share a common face, including D, N for the boundaries. Since the face are triangular, we denote by q ij,r , r " 1, . . . , R 3 the Gauss points with their associated weights ξ r .

Remark 2.1 We use the same notation q ij,r and ξ r for the two-and three-dimensional case for the sake of simplicity since it will not give rise to any confusion.

e ij n ij c i c j q ij,r n iD e iD f ij n ij c i c j q ij,r
Fig. 1. Notation for the two-dimensional (left) and three-dimensional mesh

Integration by part on cell c i Ă R 2 (resp. c i Ă R 3 ) and substituting the exact integration on edge (resp. faces) with their respective quadrature formula provide the relations

G i " ÿ jPν i |e ij | |c i | R 2 ÿ r"1 ξ r F ij,r ´ḡ i , (resp.) G i " ÿ jPν i |f ij | |c i | R 3 ÿ r"1 ξ r F ij,r ´ḡ i ,
where |e ij |, |f ij |, |c i | represent the respective measures of the geometrical entities while ḡi "

1 |c i | ż c i
gpxq dx stands for the mean value of function g on cell c i and F ij,r " εpq ij,r q∇V pq ij,r q.n ij .

We substitute the exact (but unknown) flux by a numerical approximation F ij,r at each Gauss points we shall detail in the next section and introduce the residual form

G i " ÿ jPν i |e ij | |c i | R 2 ÿ r"1 ξ r F ij,r ´gi , (resp.) G i " ÿ jPν i |f ij | |c i | R 3 ÿ r"1 ξ r F ij,r ´gi ,
where g i is a sixth-order approximation of the source term mean value using the Gauss points on cells c i .

Polynomial reconstructions

Polynomial reconstructions in the finite volume context has been introduced by [START_REF] Barth | Recent developments in high order k-exact reconstruction on unstructured meshes[END_REF] but we use the technique initially proposed for hyperbolic system [START_REF] Clain | A high-order polynomial finite volume method for hyperbolic system of conservation laws with Multi-dimensional Optimal Order Detection (MOOD)[END_REF][START_REF] Diot | Improved detection criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The MOOD method in the three-dimensional case: Very-High-Order Finite Volume Method for Hyperbolic Systems[END_REF] and adapted for the convectiondiffusion situation in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF]. Let us denote by V i an approximation of Vi "

1 |c i | ż c i V pxq dx and
V " pV i q P R I be the vector which gathers all the components (reps. V with components Vi ). The section aims to define the numerical fluxes across edges or faces in function of vector V . Basically the idea is to build local polynomial representations based on the mean value approximations V i and the boundary conditions. We have to distinguish several situations to adapt the polynomial reconstruction whether the edge (or the face) is on the boundary, on the interface Γ or is an inner edge. We rephrase the method presented in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF] for the sake of consistency and to introduce the notations but we highlight that several important new improvements will also be detailed: curved boundaries are now considered, a faster algorithm to compute values and derivatives at the Gauss points is designed and a new preconditioning matrix is proposed.

Inner face situation

We first consider the situation of cells and faces f ij (or e ij ) which do not belong to the boundary. The case for e iD and f iD as well as e iN and f iN will be treated in the next subsections. We present the reconstruction for faces (3D case) but a very similar reconstruction is available for edges as presented in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF].

Conservative reconstruction for Cell

We first start with the conservative reconstruction associated to a cell. Such a reconstruction is usually dedicated to the convection contribution but, in the present study, we need it for cells in contact with the interface Γ. Let c i be a cell, d the degree of the polynomial reconstruction we intend to build, and Spc i , dq the associated stencil constituted of cells belonging to the same sub-domain of cell c i . We assume that the approximation V i on c i is given and the polynomial function of degree d writes

p V i px; dq " V i `ÿ 1ď|α|ďd R d,α i ! px ´mi q α ´M α i ) , (2)
with α " pα 1 , α 2 , α 3 q, |α| " α 1 `α2 `α3 , x " px 1 , x 2 , x 3 q, m i the cell centroid, and we adopt the convention x α " x α 1 1 x α 2 2 x α 3 3 . We set M α i "

1 |c i | ż c i px ´bi q α dx such that the conservation property 1 |c i | ż c i p V i px; dq dx " V i (3)
holds. To fix the coefficients R d,α i of polynomial function [START_REF] Barth | Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction[END_REF] we assume that values V ℓ on cells c ℓ , ℓ P Spc i , dq are known, and we introduce the functional

E i pR d i ; dq " ÿ ℓPSpc i ,dq " 1 |c ℓ | ż c ℓ p V i px; dq dx ´Vℓ ı 2 , (4)
where R d i is the vector gathering all the coefficients R d,α i . Existence and uniqueness of the functional minimizer p R d i is guaranteed and provides the best approximation. We refer to [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF] for a detailed description of the method and the pre-conditioning strategy to obtain vector p R d i .

Conservative reconstruction for a face on interface Γ

We assume that f ij Ă Γ such that f ij " c i X c j with c i Ă Ω 1 and c j Ă Ω 2 . To compute an accurate diffusion flux across interface Γ with a continuous transition of function V , we consider a polynomial function under its conservative form

q V j px; dq " V j `ÿ 1ď|α|ďd R d,α j ! px ´mj q α ´M α j
) .

where we assume that V j is a relevant information on cell c j . We slightly modify the former functional to minimize, namely

E j pR d j ; dq " ÿ ℓPSpc j ,dq " 1 |c ℓ | ż c ℓ q V j px; dq dx ´Vℓ ı 2 `ωij " 1 |f ij | ż f ij q V j px; dq ds ´Vij ı 2 (5)
introducing a positive weight ω ij while V ij is a relevant information on face f ij . The functional takes the continuity condition into account across Γ, controlled by parameter ω ij to be fixed in the numerical applications. The minimizing vector q R d j is solution of an over-determined system, in the least-square sense, pre-conditioned as in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF].

Non-conservative reconstruction for an inner face

An inner face of Ω 1 or Ω 2 is a face of the form f ij " c i X c j such that cells c i , c j both belong to the same sub-domain. For a given degree d, we associate a stencil Spf ij , dq constituted of cells of the same sub-domain than f ij . We consider the polynomial reconstruction

r V ij px; dq " ÿ 0ď|α|ďd R d,α ij px ´mij q α , (6)
where m ij is the centroid of face f ij . Assuming that we know V ℓ on cells c ℓ , ℓ P Spf ij , dq, we introduce the functional

E ij pR d ij ; dq " ÿ ℓPSpf ij ,dq ω ij,ℓ " 1 |c ℓ | ż c ℓ r V ij px; dq dx ´Vℓ ı 2 , (7)
where ω ij,ℓ are positive weights and R d ij is the vector gathering all the coefficients R d,α ij . As mentioned in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF], weights are introduced to provide the positivity preserving property (see the numerical test 4.1) and are defined in the following way: if ℓ " i, j, ω ij,ℓ " ω ě 1 else ω ij,ℓ " 1. Vector r R d ij stands for the unique vector minimizing the functional and provides the best approximation.

Polynomial reconstruction for the Dirichlet condition on curved boundary

When dealing with curved boundaries, the substitution of domain Ω with the polygonal (or polyhedral) domain Ω h " ď i"1,...,I c i dramatically reduces the method order due to a rough approximation of the boundary with line segments. A specific treatment of the polynomial reconstruction associated to edge e iD or face f iD is required. For the Neumann case, we have assumed that Γ N is polygonal (or polyhedral) in the present study for the sake of simplicity hence no particular treatment is required.

Second-order approach

As a first approach, we recall the simple second-order method where one uses the mean value of the Dirichlet condition on the boundary edge. We present the two-dimensional situation since the extension for the three-dimensional case is straightforward. Let e iD be an edge situated on the boundary Γ D and set the mean value ViD "

1 |e iD | ż e iD
φ D psq ds. In [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF] the following conservative polynomial function of degree d has been considered

p V iD px; dq " V iD `ÿ 1ď|α|ďd R d,α iD ! px ´miD q α ´M α iD ) , (8)
taking V iD " ViD , m iD the midpoint of edge e iD and M α iD "

1 |e iD | ż e iD px ´miD q α ds such that the conservative property 1 |e iD | ż e iD p V iD px; dq ds " V iD holds.
To fix the coefficients, we introduce the functional

E iD pR d iD ; dq " ÿ ℓPSpf iD ,dq ω iD,ℓ " 1 |c ℓ | ż c ℓ p V iD px; dq dx ´φℓ ı 2 , (9)
where ω iD,ℓ are positive weights and vector p R d iD stands for the unique vector minimizing the functional and providing the best approximation. The three-dimensional situation is identical (just substitute e iD with f iD ).

If the boundary is a line (a plane), we get a d`1th-order method but such an approach achieves, at most, a second-order approximation when we substitute the mean value on a curved arc by the mean value on edge e ij . The key-point is to evaluate the polynomial reconstruction with a better V iD choice in formulae [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection diffusion problem[END_REF], different to the candidate ViD proposed above, in order to provide better approximations of p V iD px; dq in the sense we shall specify in the sequel.

3.2.2

The two-dimensional case

Geometrical ingredients

To provide a better approximation of the boundary condition, a local parametrization is introduced and a new quadrature formula is used to perform accurate numerical integrations on the arc. Let e be a generic boundary edge. We denote by v 1 , v 2 the vertices and by e " v 1 v 2 the segment with length

|v 1 v 2 | while Ŋ v 1 v 2 represents the boundary arc between v 1 and v 2 with length | Ŋ v 1 v 2 | (see Fig. 2 left panel).
We introduce the edge parametrization qptq " p1 ´tqv 1 `tv 2 , t P r0, 1s which satisfies |q 1 ptq| " |v 1 v 2 | while pptq is a parametrization of the arc such that

pp0q " v 1 , pp1q " v 2 , |p 1 ptq| " | Ŋ v 1 v 2 | is constant.
Let us denote by q 1 , . . ., q R 2 the Gauss points on edge e associated to parameters t 1 , . . ., t R 2 , then p 1 " ppt 1 q, . . ., p R 2 " ppt R 2 q are the corresponding Gauss points on the boundary arc.

v 1 v 2 q r p r c v 1 v 2 v 3 c q r p r
Fig. 2. Geometrical ingredients for the sixth-order approximation of the Dirichlet condition with curved boundary: the two-dimensional (left) and three-dimensional case (right). Points q r are the Gauss points on the edge (face) while points p r are the Gauss points on the curved boundary.

Indeed, using the quadrature rule for the numerical integration over the arc, one has

ż Ŋ v 1 v 2 V D ppqdp " ż 1 0 V D ppptqq|p 1 ptq|dt « R 2 ÿ r"1 ξ r V D pp r q|p 1 pt r q| " | Ŋ v 1 v 2 | R 2 ÿ r"1 ξ r V D pp r q.
Notice that the following property then holds for r "

1, ¨¨¨, R 2 | Ŋ v 1 p r | |v 1 q r | " | Ŋ v 2 p r | |v 2 q r | " | Ŋ v 1 v 2 | |v 1 v 2 |
which, in practice, enable to determine the location of points p r on the arc.

The polynomial reconstruction on boundary

Let us now consider u e iD the boundary arc associated to the edge e iD and assume that mean value approximations V i , i " 1, . . . , I are given. We consider the following linear operator deriving from equation ( 8)

V iD Ñ p V iD px; d, V iD q P P d
where V iD represents an approximation of the mean value on e iD . We explicitly mention the variable V iD in the polynomial function arguments since the dependence is of crucial importance in the procedure we shall implement. The main difficulty we have to face is that the Dirichlet condition is defined on u e iD and not on e iD . To overcome the problem, we introduce the functional

HpV iD q " R 2 ÿ r"1 ´p V iD pp ij,r ; d, V iD q ´VD pp ij,r q ¯2
which corresponds to the error between the polynomial approximation and the real Dirichlet condition at the Gauss points on the boundary arc. Since we are dealing with a quadratic functional, existence and uniqueness of the minimum V ‹ iD is guaranteed and p V iD px; dq " p V iD px; d, V ‹ iD q will be the polynomial reconstruction we use to compute the diffusive flux on edge e iD .

To determine V ‹ iD in practice, we propose the following simple algorithm. We consider the sequence V k iD initialized with V 0 iD " ViD and given by (1) with V k iD in hand, compute the associated polynomial function p V iD px; d, V k iD q, (2) evaluate the errors δ k r " V D pp ij,r q ´ViD pp ij,r ; d, V k iD q at the Gauss points of the arc,

(3) update the mean value on e iD with V k`1 iD

" V k iD `R2 ÿ r"1 ξ r δ k r , ( 4
) stop if |V k`1 iD ´V k iD | ă ǫ D ViD
where the tolerance ǫ D has been prescribed and set V ‹ iD " V k`1 iD , else goto step [START_REF] Ollivier-Gooch | High-Order ENO schemes for unstructured Meshes based on least-squares reconstruction[END_REF].

Numerical experiments shows that we quickly converge with two or three steps using ǫ D " 10 ´12 . Furthermore, when | ViD | ă ǫ D , we use the absolute error criterion |V k`1 iD ´V k iD | ă ǫ D in place of the relative error criterion.

Remark 3.1 Notice that the method only requires the arc length |u e iD |, the Gauss points q iD,r on the edge and the associated Gauss points p iD,r on the boundary arc. No geometrical transformation is performed which provides a very simple method, easy to implement. The main difficulty is the determination of p iD,r for a general arc of curve.

The three-dimensional case

Let f iD be a triangular face situated on the boundary Γ D and set the mean value ViD "

1 |f iD | ż f iD V D psq ds.
As in the two-dimensional case, the polynomial reconstruction p V iD with V iD " ViD only provides at most a second-order approximation with curved boundary. We propose an extention of the method to achieve an effective sixth-order accuracy.

Geometrical ingredients Let f be a boundary triangle face characterized by vertices

v 1 , v 2 , v 3 . We set f " v 1 v 2 v 3 with area |v 1 v 2 v 3 | while u f " Ŕ v 1 v 2 v 3 represents the boundary surface between v 1 , v 2 , v 3 with area | Ŕ v 1 v 2 v 3 | (see Fig. 2 right panel).
Let △ " tps, tq P r0, 1s 2 such that s `t ď 1u. We introduce the parametrization qps, tq " p1 ´t ´sqv 1 `sv 2 `tv 3 , ps, tq P △ of face f which provides

|B s q ^Bt q| " |v 1 v 2 ^v1 v 3 | " 2|v 1 v 2 v 3 |.
On the other hand, we assume that pps, tq, ps, tq P △ is a parametrization of the surface u f such that

pp0, 0q " v 1 , pp1, 0q " v 2 , pp0, 1q " v 3 , |B s p ^Bt p| " 2| Ŕ v 1 v 2 v 3 | is constant.
Let denote by q 1 , . . ., q R 3 the Gauss points on f associated to parameters ps 1 , t 1 q, . . ., ps R 3 , t R 3 q, then p 1 " pps 1 , t 1 q, . . ., p R 3 " pps R 3 , t R 3 q are the corresponding Gauss points on the boundary surface and the quadrature rule for the numerical integration over the surface writes

ż Ŕ v 1 v 2 v 3 V D ppqdp " ż △ V D ppps, tqq|B s p ^Bt p|dsdt psee that|∆| " 1 2 q « 1 2 R 2 ÿ r"1 ξ r V D pp r q|B s pps r , t r q ^Bt pps r , t r q| " | Ŕ v 1 v 2 v 3 | R 2 ÿ r"1 ξ r V D pp r q.
Notice that the following property holds for r " 1, ¨¨¨, R 3

| Ŕ v 1 v 2 p r | |v 1 v 2 q r | " | Ŕ v 2 v 3 p r | |v 2 v 3 q r | " | Ŕ v 3 v 1 p r | |v 3 v 1 q r | " | Ŕ v 1 v 2 v 3 | |v 1 v 2 v 3 | (10)
which enable, in practice, to determine points p r . Indeed, the main difficulty is the determination of the Gauss points p r associated to q r for any surface Ŕ v 1 v 2 v 3 . To this end, assume that there exists a function such that for each point

p P Ŕ v 1 v 2 v 3 provides the three areas | Ő v 1 v 2 p|, | Ő v 2 v 3 p|, | Ő v 3 v 1 p|.
Based on that function, we determine point p r associated to q r using relation [START_REF] Eymard | The finite volume method, Handbook for Numerical Analysis[END_REF] where we introduce the functional

p Ñ ˆ| Ő v 1 v 2 p| |v 1 v 2 q r | ´χ˙2 `ˆ| Ő v 2 v 3 p| |v 2 v 3 q r | ´χ˙2 `ˆ| Ő v 3 v 1 p| |v 3 v 1 q r | ´χ˙2 with χ " | Ŕ v 1 v 2 v 3 | |v 1 v 2 v 3 |
. The Gauss point p r corresponds to the minimum of the quadratic functional.

The polynomial reconstruction on boundary

Let us denote by u f iD the boundary surface associated to the triangular face f iD and assume that mean value approximations V i , i " 1, . . . , I are given. We consider one more time the following linear operator

V iD Ñ p V iD px; d, V iD q P P d
where V iD represents an approximation of the mean value on f iD while the Dirichlet condition are defined on u f iD . As in the two-dimensional case, we introduce the functional

HpV iD q " R 3 ÿ r"1 ´p V iD pp ij,r ; d, V iD q ´VD pp ij,r q ¯2
which corresponds to the error between the polynomial approximation and the real Dirichlet condition at the Gauss points of the boundary surface. Existence and uniqueness of the minimum V ‹ iD is obtained and p V iD px; dq " p V iD px; d, V ‹ iD q will be the polynomial reconstruction. The algorithm to provide V ‹ iD is the same as the one proposed in the two-dimensional case.

Numerical fluxes

With the the polynomial reconstructions in hand, we explicit the fluxes in function of the face. One has to distinguish four cases.

(1) Assume that f ij is an inner face of Ω 1 or Ω 2 , then we define the flux at the quadrature point by F ij,r " ´εpq ij,r q∇ r V ij pq ij,r ; dq.n ij .

(2) Assume that f iD belongs to Γ D , then the flux at the quadrature point writes F iD,r " ´εpq iD,r q∇ p V iD pq iD,r ; dq.n iD .

Notice that we compute the flux at the Gauss points of the edge q iD,r and not at the Gauss points p iD,r of the boundary. (3) Assume that f iN is on Γ N , then the flux at the quadrature point writes F iN,r " 0.

(4) Assume that f ij Ă Γ. We first compute p V i px; dq for the cell in c i Ă Ω 1 which shares the face f ij . Then, we evaluate the quantity

V ij " R 3 ÿ r"1 ξ r p V i pq ij,r
; dq which provides the mean value of p V i over f ij . At last, we compute q V j px; dq on the other side of the interface for cell c j Ă Ω 2 . To determine the flux across the face, we then use the formula F ij,r " ´ε2 pq ij,r q∇ q V j pq ij,r ; dq.n ji .

For any vector V P R I , we define the residual operator on cell c i by

G i pV q " ÿ jPνpiq |f ij | |c i | R 3 ÿ r"1 ζ r F ij,r pV q ´gi (11)
and introduce the affine operator GpV q " pG 1 pV q, . . . , G I pV qq T from R I into R I . The solution is vector V ‹ such that GpV ‹ q " 0 while the consistency error is given by GpV q with V " pV i q the vector of the mean values of the exact solution. We get a matrix-free problem we solve using the GMRES method (see [START_REF] Saad | Iterative methods for spares linear systems[END_REF][START_REF] Sadd | GMRES: a general minimal residual algorithm for solving nonsymmetric linear systems[END_REF] for a detail review of the method).

Improvements

Several modifications have been introduced to increase the method efficiency and reduce the computational cost in regards to the former version developed in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF].

Faster polynomial evaluation

Polynomial reconstructions are always evaluated at the Gauss points so it is a waste of computational resource to evaluate several times the same monomial functions such as px ´mi q α at the same points. To cut the time consumption, we take advantage of the linearity of the reconstruction with respect to vector V . New notations are required to present the fast evaluation since we need local indexations for the stencil and the polynomial coefficients.

Polynomial r

V ij We start with polynomial function r V ij associated to stencil Spf ij , dq. Let #Spf ij , dq be the number of elements of the stencil. We define a local indexation setting λ P t1, . . . , #Spf ij , dqu Ñ ℓpλq P Spf ij , dq.

On the other hand, we introduce the one-parameter indexation for the multi-index α. Let denote #d " dpd `1qpd `2q 6 and M pdq " tα P N 3 ; |α| ď du. We define the one-to-one mapping by κ P t1, . . . , #pd `1qu Ñ αpκq " pα 1 , α 2 , α 3 q P M pdq corresponding to the following sequence κ α ˇˇ1 p0,0,0q , 2 p1,0,0q , 3 p0,1,0q , 4 p0,0,1q , 5 p2,0,0q , 6 p1,1,0q , 7 p0,2,0q , 8 p1,0,1q , 9 p0,1,1q , 10 p0,0,2q , 11 p3,0,0q , 12 p2,1,0q , 13 p1,2,0q , 14 p0,3,0q , 15 p2,0,1q , . . . such that αpκq, κ P t1, . . . , #pd `1qu ranges all the coefficients R k,α of any polynomial function of degree d. The reciprocal operator from M pdq onto t1, . . . , #pd `1qu is given by α Ñ κpαq " #|α| `α3 pα 1 `α2 q `α3 pα 3 `3q 2 `α2 `1.

Thanks to the mapping, we introduce the local one-parameter indexation for coefficients setting r R ij,κ " r R α ij with κ " κpαq and we denote by r R ij " p r R ij,κ q κ"1,...,#pd`1q the vector which gathers all the components. In the same way, for a given stencil Spf ij , dq, we denote by W ij,λ " V ℓ with ℓ " ℓpλq and define W ij " pW ij,λ q λPt1,...,#Spf ij ,dqu the vector which gathers all the components. At last, the minimization of functional [START_REF] Cea | Vázquez-Cendón Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts[END_REF] writes in the least square sense (see [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF]) r

A ij r R ij " W ij with r A ij " p r A ij,λκ qλ"1,...,#Spf ij ,dq κ"1,...,#pd`1q
where the coefficients are determined in function of the stencil and the geometrical data of the mesh. The problem is equivalent to the matrix-vector product

r R ij " r A : ij W ij , r A : ij " p r A : ij,κλ q κ"1,...,#pd`1q
λ"1,...,#Spf ij ,dq with r A : ij the Moore-Penrose pseudo-inverse matrix computed during a pre-processing stage. Moreover the corresponding tables λ Ñ ℓpλq and κ Ñ αpκq are also stored in memory at the beginning of the computation.

For each Gauss point q ij,k , we define the vector

Q ij,r " ¨pq ij,r ´mij q αp1q pq ij,r ´mij q αp2q . . . pq ij,r ´mij q αp#pd`1qq ‹ ‹ ‹ ‹ ‹ ‹ ' P R #pd`1q .
and one has

r V ij pq ij,k , dq " #pd`1q ÿ κ"1 r R ij,κ pq ij,r ´mij q αpκq " Q T ij,r r R ij " pQ T ij,r r A : ij qW ij .
We also introduce the column vector

r χ ij,r " pr χ ij,rλ q λ"1,...,#Spf ij ,dq " Q T ij,r r A : ij P R #Spf ij ,dq ,
and the polynomial value at the Gauss point is then evaluated with the product

r V ij pq ij,r ; dq " #Spf ij ,dq ÿ λ"1 r χ ij,rλ V ℓpλq " r χ T ij,r W ij . (12)
We now deal with the gradient of the polynomial reconstruction. We introduce the #pd `1qlength vector of elements of R 3 (or a 3 ˆ#pd `1q real matrix)

JQ ij,r " ¨¨0 0 0 ‹ ‹ ‹ ' , ¨1 0 0 ‹ ‹ ‹ ' , ¨0 1 0 ‹ ‹ ‹ ' , ¨0 0 1 ‹ ‹ ‹ ' , . . .
¨α1 pq ij,r ´mij q α´p1,0,0q α 2 pq ij,r ´mij q α´p0,1,0q α 3 pq ij,r ´mij q α´p0,0,1q

‹ ‹ ‹ ' , . . . ‹ ‹ ‹ '
with the rule α ´β " pα 1 ´β1 , α 2 ´β2 , α 3 ´β3 q. In other words, for any κ " 1, ¨¨¨, #pd `1q, we have the R 3 -value entries JQ ij,rκ " ¨α1 pκqpq ij,r ´mij q αpκq´p1,0,0q α 2 pκqpq ij,r ´mij q αpκq´p0,1,0q

α 3 pκqpq ij,r ´mij q αpκq´p0,0,1q ‹ ‹ ‹ ' .
Applying the same procedure, we deduce that

∇ r V ij pq ij,k , dq " #pd`1q ÿ κ"1 JQ ij,rκ p r A : ij W ij q κ .
For each Gauss point q iD,k P f iD , we define vector

Q iD,r "
¨pq iD,r ´miD q αp1q pq iD,r ´miD q αp2q

. . .

pq iD,r ´miD q αp#pd`1qq ‹ ‹ ‹ ‹ ‹ ‹ ' P R #pd`1q .
and one has

p V iD pq iD,k , dq " #pd`1q ÿ κ"1 p R iD,κ pq iD,r ´miD q αpκq " Q T iD,r p R iD " pQ T iD,r p A : iD qW iD .
Setting p χ iD,r " pp χ iD,r,λ q λ"0,...,#Spf iD ,dq " Q T iD,r p A : iD P R #Spf iD ,dq`1 , the polynomial value at each Gauss point is then evaluated with the product

p V iD pq iD,r ; dq " V iD `#Spf iD ,dq ÿ λ"1 p χ iD,rλ V ℓpλq .
For the gradient evaluation, we introduce the vector of elements of R 3

JQ iD,r " ¨¨0 0 0 ‹ ‹ ‹ ' , ¨1 0 0 ‹ ‹ ‹ ' , ¨0 1 0 ‹ ‹ ‹ ' , ¨0 0 1 ‹ ‹ ‹ ' , . . .
¨α1 pq iD,r ´miD q α´p1,0,0q α 2 pq iD,r ´miD q α´p0,1,0q α 3 pq iD,r ´miD q α´p0,0,1q

‹ ‹ ‹ ' , . . . ‹ ‹ ‹ '
and we write

∇ p V iD pq iD,k , dq " #pd`1q ÿ κ"1 JQ iD,rκ p p A : iD W iD q κ .
Setting for λ " 0, ¨¨¨, #Spf iD , dq

x Jχ iD,rλ " #pd`1q ÿ κ"1 JQ iD,rκ p A : iD,κλ P R 3 ,
the polynomial gradient at each Gauss point then corresponds to the product

∇ p V iD pq iD,r ; dq " x Jχ iD,r0 V iD `#Spf iD ,dq ÿ λ"1 x Jχ iD,rλ V ℓpλq " #Spf iD ,dq ÿ λ"1
x Jχ iD,rλ V ℓpλq .

A new preconditioning matrix

Due to the linearity of the problem and the reconstruction process, operator GpV q is an affine mapping V Ñ GpV q " AV ´b with b P R I and A a I ˆI real numbers matrix. In [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF], a new the system P : pi, iqP pi, iq `ÿ jPν i P : pi, jqP pj, iq " 1, P : pi, jqP pj, jq `P : pi, iqP pi, jq " 0, @j P ν i .

After some algebraic manipulations, we obtain an explicit expression for P : , namely P : pi, jq " ´P pi, jq P : pi, iq P pj, jq , j P ν i , with P : pi, iq " 1 P pi, iq ´ÿ jPνpiq P pi, jqP pj, iq P pj, jq .

Numerical tests

The present section is dedicated to the quantitative and qualitative assessments of the scheme robustness and accuracy. We first tackle the stability question, then we propose convergence tests for regular and discontinuous diffusion coefficients. Numerical simulations with curved boundary are carried out to show the critical issues that a specific treatment is mandatory for the polynomial reconstruction to preserve the optimal order. At last, we present some numerical experiences deriving from a concrete physical application to highlight the performances of the improvements proposed in section 3.4

To compare the numerical approximation with the exact solution, we introduce two kinds of errors, namely:

E 1 " I ÿ i"1 |c i ||φ i ´φ i |, E 8 " max i"1,¨¨¨,I |φ i ´φ i |,
while the convergence order between two successive meshes of N 1 and N 2 cells is given by

O 1 " d lnpE 1 1 {E 1 2 q lnpN 2 {N 1 q , O 8 " d lnpE 8 1 {E 8 2 q lnpN 2 {N 1 q ,
with d the space dimension. In this study, we employ triangular meshes for the two-dimensional geometries and tetrahedral meshes for the three-dimensional simulations. Figures are obtained using the gmsh software [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF].

Maximum principle

We first assess the scheme ability to preserve the positivity or satisfy the maximum principle. Indeed, in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF], numerical tests show that the choice of the weights in the non-conservative reconstruction process is of crucial importance. We have proceeded in the same way as [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF],

experimenting several weight values and check the non-negativity of the inverse matrix of the global linear system. To this end, let consider the unit-length cube. Function V px, y, zq " sinpπxq sinpπyq sinpπzq is the solution of problem ´∆V " g with g " 3π 2 sinpπxq sinpπyq sinpπzq and homogeneous boundary condition V D " 0 hence the exact solution ranges in the interval r0, 1s. The polynomial reconstruction process on face f ij depends on weights ω ij,ℓ we introduce in the functional [START_REF] Cea | Vázquez-Cendón Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts[END_REF]. We fix the values in the following way: if ℓ " i, j we prescribe ω ij,ℓ " ω whereas ω ij,ℓ " 1 for the other cells. All the numerical tests have been performed with the P 1 reconstruction but the conclusion holds with a P 3 or a P 5 reconstruction.

Table 1 Minimum and maximum of the approximations with respect to the weight with a 2321 cells mesh (left) and a 6327 cells mesh (right). We have use the P 1 reconstruction to compute the diffusive flux. Table 1 provides the minimum and the maximum of the approximations in function of ω for two different meshes. We plot in Fig. 3 three simulations performed with three different values of ω. We clearly observe that the maximum principle is not fulfilled for ω " 1.7 and ω " 2 but, in the second case, we obtain a lower number of non-physical cells i.e. cells associated to a negative value or a value larger to 1. For ω large enough, (here ω " 4), the inverse matrix has all the coefficients non-negative and the numerical solution ranges in r0, 1s. In the following tests, we shall use the value ω " 4 to fulfil the maximum principle and provide the stability.

Convergence test

Convergence tests are carried out to check the order of accuracy of the scheme in the threedimensional context. In the unit-length cube domain, we compute approximations of the solution for problem ´∆V " g with g " 3π 2 sinpπxq sinpπyq sinpπzq and homogeneous boundary condition V D " 0 that we compare with the exact solution V " sinpπxq sinpπyq sinpπzq. We have experimented three reconstructions and reported in Table 2 the error in L 1 -and L 8 -norm with the associated convergence orders. We obtain an effective second-order, fourth-order and sixth-order of convergence using the P 1 , P 3 and P 5 reconstructions respectively. To highlight the efficiency of very high-order methods, a simple extrapolation shows that the error of 4.84e ´08 obtained with the sixth-order method with 32408 tetrahedra will require around 63 billions of cells to achieve the same error with a second-order method and a huge computational effort.

Simulation with discontinuous dielectric coefficients

We tackle the important situation where the domain is constituted of two sub-domains with different dielectric coefficients leading to a discontinuity of the electrical field at the interface. Of course the test can be generalized to several domains with discontinuous coefficients. We consider the unit-length cube such that ε 1 " 1 in Ω 1 " tx, y, z ă 0.5u and ε 2 " 10 in Ω 2 " tx, y, z ą 0.5u i.e. the discontinuity interface Γ is located in the plane z " 0.5. The potential function V is composed of two functions V 1 " V |Ω 1 and V 2 " V |Ω 2 and the flux conservation condition across the interface Γ writes

ε 1 ∇V 1 .n Γ " ε 2 ∇V 2 .n Γ , (14)
while we assume the continuity of the potential function V 1 " V 2 on Γ.

One can check that the functions with a " 1

V 1 px, y, zq " A ε 1 ˆ1 π ˙2 sinpπzq `a ε 1 z `VD1 , V 2 px, y, zq " A ε 1 ˆ1 π ˙2 sinpπzq `a ε 2 pz ´1q `VD2 x z Ω 1 Ω 2 Γ Γ D2 Γ D1 Γ N
1 2ε 1 `1 2ε 2 ˜A ˆ1 π ˙2 ˆ1 ε 2 ´1 ε 1 ˙`V D2 ´VD1 ¸,
are the solution of the electrical problem with the source term g " A sinpπzq with homogeneous Neumann condition on the lateral sides Γ N (see Fig. 4) while we prescribe the Dirichlet condition V px, y, 0q " V D1 on Γ D1 and V px, y, 1q " V D2 on Γ D2 . All the numerical simulations are carried out using A " 1, V D1 " 1 and V D2 " 10 Table 3 shows that we obtain an effective second-, fourth-and sixth-order of convergence when using the P 1 , P 3 and P 5 reconstructions respectively. The low regularity of the solution in the vicinity of Γ is well-treated by the specific reconstruction given in section 3.1.2 and does not keep the scheme from achieving very high-order approximations.

A non-homogeneous test case

We now assess the convergence order when dealing with a non-constant smooth dielectric coefficient and non-homogeneous boundary conditions. The goal is to check that the scheme effectively provides the optimal convergence order. To this end, we assume that the solution is V px, y, zq " lnp1 `x `yqzp1 ´zq and the dielectric coefficient is given by εpx, y, zq " 1 `x `y.

Hence, one has to solve ´∇ ¨pε∇V q " g with the right-hand side source term gpx, y, zq " 2 lnp1 `x `yqp1 `x `yq and the Dirichlet condition is given by the exact solution on the whole boundary. We display a cut of the numerical solution in Fig. 5 while we report in Table 4 the L 1 and L 8 errors and convergence rates. Effective second-, fourth-and sixth-order of accuracy are obtained for the L 1 norm while the accuracy with the L 8 norm is slightly worse but no instabilities are reported. To highlight the efficiency of the very high-order scheme, we observe that the sixthorder scheme cuts the error with an order of six magnitudes (in L 1 -norm) in relation with the second-order one (for example with the 32048 cells mesh). An important goal of the present study is to achieve very high accuracy even with curved boundary. Indeed, a formal sixth-order scheme may be spoiled by a rough approximation of the boundary condition leading to, at most, a second-order method. We start with the twodimensional case and assess the efficiency of the methods detailed in section 3.2.2 where the straightforward approximation on edge (Dirichlet condition applied on edges) is substituted by an algorithm to recover the optimal convergence order. We consider an annulus of external radius R 1 " 10 ´6 and internal radius r 2 " 10 ´7. We prescribe Dirichlet condition on the boundaries, namely V " V 1 , V " V 2 on the external and internal circles respectively. The solution of problem ´∆V " 0 is given by V px, yq " a lnp ?

x 2 `y2 q`b with a " V 1 ´V2 lnpR 1 q ´lnpR 2 q and b "

V 2 lnpR 1 q ´V1 lnpR 2 q lnpR 1 q ´lnpR 2 q (see Fig. 6). Table 5 reports the errors and convergence orders using the simple approximation of the Dirichlet condition as proposed in section 3.2.1 (the Dirichlet condition is directly involved in the mean value over the face) while we report in Table 6 the errors and convergence orders when using the correction proposed in section 3.2.2.

Numerical evidences show that prescribing the Dirichlet on the edge in a straightforward man- ner provides at most a second-order scheme where the limitation is clearly due to a poor evaluation of the Dirichlet condition define on a curved boundary. The technique we propose to overcome that limitation provides an excellent result since we recover the optimal sixth-order in the case of the P 5 reconstruction (and the fourth-order with the P 3 reconstruction). Such an example underlines that a correct approximation of the boundary condition in the case of curved boundary is crucial.

Three-dimensional problem with curved boundary

We now turn to the three-dimensional case following the strategy proposed for the two-dimensional annulus. The domain is situated between two spheres of radius R 1 " 1 and R 2 " 0.1 and we prescribe the Dirichlet condition V " V 1 and V " V 2 respectively. The solution of problem

´∆V " 0 is given by V px, y, zq " a r `b with a " R 1 R 2 pV 1 ´V2 q R 2 ´R1 and b " R 2 V 2 ´R1 V 1 R 2 ´R1 as
shown in Fig. 7.

In Table 7, we report the errors and convergence orders using the straightforward approximation (the Dirichlet condition is directly involved in the mean value over the face) while we report in Table 8 the errors and convergence orders using the correction proposed in section 3.2.3. Alike the two-dimensional case, the poor approximation of the boundary ruins the efforts of the polynomial reconstruction to achieve high accuracy whereas the corrected version recover the optimal orders. To end the section we propose a concrete example and assess the efficiency of the new polynomial reconstruction detailed in section 3.4.1. The micro-tip is the main piece in Atomic Force with V the electrical potential, ρ the charge density, ε the dielectric coefficient and E " ´∇V the electric field. We prescribe V " 0 on Γ 0 , V " V b on Γ 1 , ∇V ¨n " 0 on Γ N , ε a " 8.85 10 ´12 , ε d " 10ε a .

We here highlight the two main difficulties of the exercise: the dielectric coefficient is discontinuous at the interface Γ which shares the two domains and the tip extremity is spherical hence the boundary condition is not supported by the polyhedral mesh boundary. To overcome the difficulties, we have applied the techniques tested in section 4.3 and section 4.6. Numerical simulations have been carried out with a 88645 tetrahedra mesh where the tip extremity has been locally refined. Fig. 9 displays the electrical potential (left panel) and the electrical field (right panel) for a bias V b " 25 V .

To assess the performances of the improvements proposed in section 3.4, we have computed the electrical potential following two different methods. Method M1 uses the fast reconstruction where polynomial functions are pre-calculated at the Gauss points (see section 3.4.1) and used the preconditioning matrix proposed in section 3.4.2. Method M2 is the original (or former) method where we recompute all the polynomial coefficients while the preconditioning matrix is based on an ersatz of Patankar matrix (see [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF]). Table 9 provides the computational time (in second) and the number of iterations for both the methods using the P 1 , P 3 and P 5 reconstructions. Numerical simulations are carried out with three levels of meshes: very coarse, coarse and fine to evaluate the computational performances. Notice that method M2 with the P 5 reconstructions has fail (the GMRES procedure does not converge) or the computational time is too large. Clearly, the new version requires three or four times less iterations to converge (the residual norm in the GMRES method is lower than 10 ´10 ) and the computational time is cut by ten or more. We also observe that the computational effort of the new method is less sensitive to the polynomial degree in comparison with the original method. Notice that we do not have taken advantage of the high capacity of parallelization of the scheme (each evaluation of the flux is independent) which will provide very good speed-up.

Fig. 3 .

 3 Fig. 3. Maximum principle property for different weights ω: the adjacent cells have a weight equal to 1.7 (left), the adjacent cells have a weight equal to 2 (middle) while the right panel displays the approximation for ω " 4 (very close to the exact solution). White cells correspond to negative values while black cells correspond to values larger than one.

Fig. 4 .

 4 Fig. 4. The domain is constituted of two sub-domains with constant dielectric coefficients with a discontinuity at the interface Γ. We prescribe Dirichlet conditions on Γ D1 and Γ D2 and homogeneous Neumann condition on the lateral faces Γ N .

Fig. 5 .

 5 Fig. 5. A view of the solution for the Poisson equation with non-constant dielectric coefficient and a non-homogeneous boundary condition.

Fig. 6 .

 6 Fig. 6. Solution of the Poisson equation for the annulus problem.

Fig. 7 .

 7 Fig. 7. Numerical approximation of the electrical potential between two spheres.

Fig. 8 .

 8 Fig. 8. The micro-tip.

Fig. 8

 8 Fig. 8 right panel sketches a longitudinal cut of the three-dimensional computational domain Ω constituted of air (domain Ω a q and a dielectric (domain Ω d ). The boundary is partitioned into Γ 0 where we prescribe the null voltage condition, Γ 1 where we prescribe the bias voltage while we have an homogeneous Neumann condition on Γ N since no current crosses the boundary. The model then writes ´∇ ¨pε∇V q " ρ

Fig. 9 .

 9 Fig. 9. The electrical potential (left) and electrical field (right) around the tip.

Table 2 L

 2 1 and L 8 errors and convergence rates for a smooth solution of the Poisson problem.

	Nb of		P 1				P 3				P 5
	Cells	err 1		err 8		err 1		err 8		err 1		err 8
	2321	4.31e-03	-	1.93e-02	-	2.03e-04	-	1.72e-03	-	1.63e-05	-	7.61e-05	-
	6371	2.28e-03	1.9	1.18e-02	1.5	4.95e-05	4.2	4.64e-04	3.9	1.52e-06	7.1	1.52e-05	4.8
	16020	1.24e-03	2.0	6.09e-03	2.2	1.51e-05	3.9	1.88e-04	2.9	1.92e-07	6.7	2.90e-06	5.4
	32408	7.54e-04	2.1	4.33e-03	1.5	5.41e-6	4.4	7.01e-05	4.2	4.84e-08	5.9	8.34e-07	5.3

Table 3

 3 Poisson problem with discontinuous coefficients. L 1 and L 8 errors and convergence rates for the P 1 , P 3 and P 5 reconstructions.

	Nb of		P 1				P 3				P 5
	Cells	err 1		err 8		err 1		err 8		err 1		err 8
	2495	4.55e-02	-	2.70e-01	-	1.23e-04	-	7.58e-04	-	1.30e-05	-	9.50e-06	-
	5404	2.81e-02	1.87	1.79e-01	1.6	5.15e-05	3.38	2.78e-04	3.9	3.36e-06	5.25	1.57e-06	7.0
	11546	1.59e-02	2.25	1.01e-01	2.25	1.74e-05	4.29	1.24e-04	3.19	6.25e-07	6.65	4.26e-07	5.65
	24848	9.2e-03	2.14	9.58e-02	2.15	5.54e-06	4.48	3.84e-05	4.6	1.52e-07	5.55	1.08e-07	4.9

Table 4 L

 4 1 and L 8 errors and convergence rates using non-constant dielectric coefficient and non-homogeneous Dirichlet condition.

	Nb of		P 1				P 3				P 5
	Cells	err 1		err 8		err 1		err 8		err 1		err 8
	2321	1.03e-03	-	4.76e-03	-	1.72e-06	-	2.01e-05	-	2.75e-08	-	6.50e-06	-
	6327	5.01e-04	2.16	2.85e-03	1.53	4.41e-07	4.07	5.07e-06	4.12	3.2 e-09	6.43	1.39e-06	4.61
	16020	2.7e-04	2.0	1.39e-03	2.32	1.32e-07	3.9	2.05e-06	2.92	5.7 e-10	5.57	2.49e-07	5.55
	32048	1.7e-04	2.0	9.5e-04	1.65	4.82e-08	4.36	8.71e-07	3.7	1.32e-10	6.32	7.7e-07	5.08

Table 5 L

 5 1 and L 8 errors and convergence rates for the two-dimensional annulus using the straightforward approximation of the Dirichlet condition.

	Nb of		P 1				P 3					P 5
	Cells	err 1		err 8		err 1		err 8		err 1		err 8
	804	9.26E-02	-	1,95E-01	-	3,80E-02	-	1,50E-01	-	5.10E-02	-	1,38E-01	-
	2178	3.09E-02	2.0	7.50E-02	1.9	9.46E-03	2.8	1.54E-02	4.1	9.60E-03	3.0	1.46E-02	4.10
	8226	7.17E-03	2.2	2.14E-02	1.9	2.09E-03	2.3	4.43E-03	1.9	1.98E-03	2.4	4.19E-03	1.88
	24502	2.92E-03	1.7	9.98E-03	1.1	8.05E-04	1.4	2.50E-03	0.9	7.79E-04	1.3	2.46E-03	0.75

Table 6 L

 6 1 and L 8 errors and convergence rates for the two-dimensional annulus using the correction method for the Dirichlet condition.

	Nb of		P 1				P 3				P 5
	Cells	err 1		err 8		err 1		err 8		err 1		err 8
	804	6.40E-02	-	1,66E-01	-	1.98E-02	-	1.39E-01	-	2.40E-02	-	1.28E-01	-
	2178	2.20E-02	2.2	6.65E-02	1.9	1.00E-03	5.1	7.48E-03	5.9	8.84E-04	6.6	4.34E-03	6.8
	8226	5.19E-03	2.2	1.81E-02	2.0	4.87E-05	4.6	6.59E-04	3.7	4.99E-06	7.8	6.33E-05	6.0
	24502	1.20E-03	2.1	4.35E-03	2.1	3.60E-06	3.7	7.49E-04	3.1	9.09E-08	5.9	2.26E-06	5.5

Table 7 L 1

 71

	Nb of		P 1				P 3				P 5
	Cells	err 1		err 8		err 1		err 8		err 1		err 8
	2924	5.19E-01	--	9.11E-01	--	6.41E-02	--	3.17E-01	--	2.24E-01	--	7.83E-01	--
	7917	3.02E-01	1.63	5.87E-01	1.32	2.52E-02	2.81	1.59E-01	2.07	8.77E-02	2.82	2.96E-01	2.93
	17901	1.96E-01	1.59	4.30E-01	1.15	1.27E-02	2.52	7.72E-02	2.66	4.39E-02	2.54	9.15E-02	4.32
	31976	1.50E-01	1.38	3.06E-01	1.75	7.47E-03	2.75	4.33E-02	2.99	2.33E-02	3.28	5.90E-02	2.27

and L 8 errors and convergence rates for the three-dimensional case using the straightforward method for the Dirichlet condition

Table 8 L

 8 1 and L 8 errors and convergence rates for the three-dimensional case using the correction method for the Dirichlet condition

	Nb of		P 1				P 3				P 5
	Cells	err 1		err 8		err 1		err 8		err 1		err 8
	2924	4.17E-01	-	8.34E-01	-	3.94E-02	-	2.77E-01	-	1.55E-01	-	5.60E-01	-
	7917	2.06E-01	2,1	5.20E-01	1,4	1.07E-02	3.9	9.47E-02	3.2	7.33E-03	9.2	2.89E-02	8.9
	17901	1.21E-01	1,9	3.27E-01	1,7	3.50E-03	4.1	2.90E-02	4.4	1.17E-03	6.8	4.54E-03	6.8
	31976	8.11E-02	2,1	2.40E-01	1,6	1.72E-03	3.7	1.40E-02	3.8	4.11E-04	5.4	1.73E-03	5.0

Table 9

 9 Computational time (in seconds) and number of iterations between the M1 (improved) and the M2 (original) method.

	Nb of		P 1				P 3				P 5
	Cells	M1		M2		M1		M2		M1		M2
		time	iter	time	iter	time	iter	time	iter	time	iter	time	iter
	3663	3	51	98	215	15	105	181	288	38	200	-	-
	11339	13	75	361	237	27	147	711	340	112	237	-	-
	22590	55	93	697	311	65	183	1263	512	204	321	-	-

Microscopy (AFM) set-up[START_REF] Binnig | Atomic Force Microscope[END_REF]. Fig.8left panel represents a general view of the micro-tip while the middle panel shows a zoom of the tip extremity which is the critical part of the piece. We aim at computing the electrostatic field since electrostatic interactions at the nano-scale are among those of greatest strength and longer interaction distance, compared to the other relevant forces such as van der Waals and capillarity[START_REF] Colchero | Resolution enhancement and improved data interpretation in electrostatic force microscopy[END_REF]. We use a realistic model based on the Poisson equation to reproduce the electrostatic interaction in AFM working conditions[START_REF] Cappella | Force-distance curves by atomic force microscopy[END_REF].Γ N Γ 0 Γ 1 Ω a Ω d Γ

Acknowledgements

The first and third authors are supported by Region Midi-Pyrénées (France). The second author research was financed by FEDER Funds through Programa Operacional Factores de Competitividade -COMPETE and by Portuguese Funds through Fundação para a Ciência e a Tecnologia, within the project PTDC/MAT/121185/2010, the project FCT-ANR/MAT-NAN/0122/2012 and the strategic programme PEst-OE/MAT/UI0013/2014.

Setting matrix Ă

Jχ ij,r with entries

the polynomial gradient at each Gauss point then corresponds to the product ∇ r V ij pq ij,r ; dq "

The main advantage is that coefficients r χ ij,rλ P R and Ă Jχ ij,rλ P R 3 are pre-computed just one time at the initial stage. Hence, evaluations of r V ij pq ij,r ; dq or ∇ r V ij pq ij,r ; dq turn to be very efficient and the storage is strongly reduced since we just need several small vectors for each Gauss point to use relations [START_REF] Hernández | High-order finite volume schemes for the advection-diffusion equation[END_REF] or [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF]. We dramatically reduce the computational effort since neither the polynomial coefficients, nor the monomial expressions have to be recomputed. Moreover the Moore-Penrose matrices used in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF] are now useless. We save a lot of memory since the coefficients storage is really lower that the storage of the Moore-Penrose matrix r A : ij and the polynomial coefficients r R ij for each face f ij .

Polynomial p

V iD We now turn to the case of a conservative reconstruction for the faces associated to the Dirichlet condition. The main difference with the former case is the indexation since we the value on face f iD is prescribed. Let #Spf iD , dq be the number of stencil elements. We define a local indexation λ P t1, . . . , #Spf iD , dqu Ñ ℓpλq P Spf iD , dq and introduce vector W iD P R #Spf iD ,dq`1 adding the extra-component W 0 " V iD and W iD,λ " V ℓpλq for λ P t1, . . . , #Spf iD , dqu. Using similar notations, vector

is evaluated with p R iD " p A : iD W iD where p A : iD " p p A : iD,κλ q κ"1,...,#pd`1q λ"0,...,#Spf iD ,dq .

Notice that the conservative property yields p R iD,1 " V iD .

preconditioning matrix method has been proposed taking advantage that the problem derived from the finite volume method. With coefficients x Jχ ij,rλ in hand, we propose an other efficient preconditioning matrix, that is a matrix P : P R IˆI such as one solves the problem P : GpV q " 0 in place of GpV q " 0. Clearly the two problems are equivalent (P : is assumed non-singular) but we expect a better conditioning problem by multiplying with P : . To this end, we slightly modified some previous notations and introduce new ingredients. For a given cell c i and j P νpiq, the face f ij is associated to a stencil Spf ij , dq. Let #Spf ij , dq be the number of elements of the stencil. We redefine the local indexation setting λ P t1, . . . , #Spf ij , dqu Ñ ℓpi, j; λq P Spf ij , dq to highlight the dependence of function ℓpi, j; .q to face f ij . Since the function is one-to-one, there exists a unique reciprocal function ℓ P Spf ij , dq Ñ λpi, j; ℓq P t1, . . . , #Spf ij , dqu.

Index i always belongs to Spf ij , dq by construction, then λpi, j; iq is the local index associated to i. For j P νpiq, if f ij is an inner cell, i.e. cell c j exists, then by construction j P Spf ij , dq and λpi, j; jq is the local index associated to j. If f ij is on the boundary, we skip the index and denote by νpiq Ă νpiq the index set such that f ij only are inner faces. Matrix P P R IˆI is then defined as follow: for the diagonal coefficient we set

and the extra-diagonal coefficients for j P ν i writes

Matrix P is not necessary symmetric but have a symmetric structure i.e. the non-null entries are symmetric. In some extent, matrix P is the main part of matrix A where we have performed the assembly procedure only for specific entries, namely the diagonal and the extra-diagonal entries corresponding to the adjacent cells.

To provide the preconditioning matrix, we consider a matrix P : with the same null entries than matrix P (the other entries still unknown) and denote by C " P : P the product of the two matrices. Since a face f ij , j P ν i , only shares two cells, the following equalities hold Cpi, iq " P : pi, iqP pi, iq `ÿ jPν i P : pi, jqP pj, iq, Cpi, jq " P : pi, jqP pj, jq `P : pi, iqP pi, jq, @j P ν i .

To fix the P : matrix coefficients, we state that C " I only for the non-null entries and we get

Conclusions

We have proposed several extensions of the sixth-order finite volume method initially developed for two-dimensional geometries in [START_REF] Clain | A sixth-order finite volume method for multidomain convection-diffusion problem with discontinuous coefficients[END_REF]: an accurate method to prescribe the Dirichlet condition on bended boundary for the two-and the three-dimensional case which provide the optimal order; a new strategy to reduce the computational effort. Numerical tests have been carried to assess the technique efficiency in particular the question of curved boundary has been tackled. Extension of the method for the Stokes problem is the next step to demonstrate the versatility and the competitiveness of very high-order finite volume methods.