François Dubois
email: francois.dubois@u-psud.fr

Isabelle Greff
email: isabelle.greff@univ-pau.fr

Charles Pierre
email: charles.pierre@univ-pau.fr

Greff Pierre Dubois

Raviart-Thomas finite elements of Petrov-Galerkin type

Keywords: inf-sup condition, finite volumes, mixed formulation. Subject classification: 65N08, 65N12, 65N30

come

Introduction

Finite volume methods are very popular for the approximation of conservation laws. The unknowns are mean values of conserved quantities in a given family of cells, also named "control volumes". These mean values are linked together by numerical fluxes. The fluxes are defined and computed on interfaces between two control volumes. They are defined with the help of cell values on each side of the interface. For hyperbolic problems, the computation of fluxes is obtained by linear or nonlinear interpolation (see e.g. Godunov et al. [START_REF] Godounov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF]).

This paper addresses the question of flux computation for second order elliptic problems. To fix the ideas, we restrict ourselves to the Laplace operator. The computation of flux is held by differentiation: the interface flux must be an approximation of the normal derivative of the unknown function at the interface between two control volumes. The computation of diffusive fluxes using finite difference formulas on the mesh interfaces has been addressed by much research for more than 50 years, as detailed below. Observe that for problems involving both advection and diffusion, the method of Spalding and Patankar [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF] defines a combination of interpolation for the advective part and derivation for the diffusive part.

The well known two point flux approximation (see Faille, Gallouët and Herbin [START_REF] Faille | Des mathématiciens découvrent les volumes finis[END_REF][START_REF] Herbin | An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh[END_REF]) is based on a finite difference formula applied to two scalar unknowns on each side of the interface. These unknowns are ordered in the normal direction of the interface considering a Voronoi dual mesh of the original mesh, [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques[END_REF]. When the mesh does not satisfy the Voronoi condition, the normal direction of the interface does not coincide with the direction of the centres of the cells. The tangential component of the gradient needs to be introduced. We refer to the "diamond scheme" proposed by Noh in [START_REF] Noh | CEL: A time-dependent, two-space-dimensional, coupled Euler-Lagrange code[END_REF] in 1964 for triangular meshes and analysed by Coudière, Vila and Villedieu [START_REF] Coudière | Convergence rate of a finite volume scheme for a two-dimensional convection-diffusion problem[END_REF]. The computation of diffusive interface gradients for hexahedral meshes was studied by Kershaw [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF], Pert [START_REF] Pert | Physical constraints in numerical calculations of diffusion[END_REF] and Faille [START_REF] Faille | A control volume method to solve an elliptic equation on a two-dimensional irregular mesh[END_REF]. An extension of the finite volume method with duality between cells and vertices has also been proposed by Hermeline [START_REF] Hermeline | Une méthode de volumes finis pour les équations elliptiques du second ordre[END_REF] and Domelevo and Omnes [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF].

The finite volume method has been originally proposed as a numerical method in engineering [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF][START_REF] Rivas | Be03, programme de calcul tridimensionnel de la transmission de chaleur et de l'ablation[END_REF]. Eymard et al. (see e.g. [START_REF] Eymard | Finite volume methods. Handb. Numer. Anal[END_REF]) proposed a mathematical framework for the analysis of finite volume methods based on a discrete functional approach. Even if the method is non consistent in the sense of finite differences, they proved convergence. Nevertheless, a natural question is the reconstruction of a discrete gradient from the interface fluxes. This question has been first considered for interfaces with normal direction different to the direction of the neighbour nodes by [START_REF] Noh | CEL: A time-dependent, two-space-dimensional, coupled Euler-Lagrange code[END_REF][START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF][START_REF] Pert | Physical constraints in numerical calculations of diffusion[END_REF][START_REF] Faille | A control volume method to solve an elliptic equation on a two-dimensional irregular mesh[END_REF]. From a mathematical point of view, a natural condition is the existence of the divergence of the discrete gradient: how to impose the condition that the discrete gradient belongs to the space H(div) ? If this mathematical condition is satisfied, it is natural to consider mixed formulations. After the pioneering work of Fraeijs de Veubeke [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF], mixed finite elements for two-dimensional space were introduced by Raviart and Thomas [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF] in 1977. They will be denoted as "RT" finite elements in this contribution.

The discrete gradient built from the RT mixed finite element is non local. Precisely, this discrete gradient for the mixed finite element method of a scalar shape function u is defined as the unique p := ∇ h u ∈ RT so that (p, q) 0 = -(u, div q) 0 for all q ∈ RT. With this definition, the flux component of p for a given mesh interface cannot be computed locally using only the values of u in the interface neighbourhood. This is not suitable for the discretisation of a differentiation operator that is essentially local. In their contribution [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF], Baranger, Maitre and Oudin proposed a mass lumping of the RT mass matrix to overcome this difficulty. They introduced an appropriate quadrature rule to approximate the exact mass matrix. With this approach, the interface flux is reduced to a true two-point formula. Following the idea in [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF], for general diffusion problems, further works have investigated the relationships between local flux expressions and mixed finite element methods. Arbogast, Wheeler and Yotov in [START_REF] Arbogast | Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences[END_REF] present a variant of the classical mixed finite element method (named expanded mixed finite element). They shown that, in the case of the lowest order Raviart-Thomas elements on rectangular meshes, the approximation of the expanded mixed finite element method using a specific quadrature rule leads to a cell-centered scheme on the scalar unknown. That scheme involves local flux expressions based on finite difference rules. The results in [START_REF] Arbogast | Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences[END_REF] were extended by Wheeler and Yotov in [START_REF] Wheeler | A multipoint flux mixed finite element method[END_REF] for the classical mixed finite element method. The multipoint flux approximation methods propose to evaluate local fluxes with finite difference formula, see e.g. Aavatsmark [START_REF] Aavatsmark | An introduction to multipoint flux approximations for quadrilateral grids[END_REF]. That method has been later shown in Aavatsmark et al. [START_REF] Aavatsmark | Convergence of a symmetric mpfa method on quadrilateral grids[END_REF] to be equivalent on quadrangular meshes with the mixed finite element method with low order elements implemented with a specific quadrature rule. Local flux computation using the Raviart-Thomas basis functions has also been developed by Younès et al. in [START_REF] Younès | A new formulation of the mixed finite element method for solving elliptic and parabolic PDE with triangular elements[END_REF]. That question has been further investigated by Vohralík in [START_REF] Vohralík | Equivalence between lowest-order mixed finite element and multi-point finite volume methods on simplicial meshes[END_REF]. He shows that the mixed finite element discrete gradient p = ∇ h u can be computed locally with the help both of u and of the source term f h := -div(∇ h u) (that depends on ∇ h u). More precisely, with a slight modification of the discrete source term f h in the finite volume method, it has been proven in [START_REF] Chavent | On the finite volume reformulation of the mixed finite element method for elliptic and parabolic PDE on triangles[END_REF][START_REF] Vohralík | Mixed finite element methods: implementation with one unknown per element, local flux expressions, positivity, polygonal meshes, and relations to other methods[END_REF] that the two discrete gradients defined either with the mixed finite element or with the finite volume method are identical. Moreover, in case of a vanishing source term f = 0, the two discrete gradients are identical without any modification of the discrete source term f h (see [START_REF] Younès | A new formulation of the mixed finite element method for solving elliptic and parabolic PDE with triangular elements[END_REF], [START_REF] Chavent | On the finite volume reformulation of the mixed finite element method for elliptic and parabolic PDE on triangles[END_REF], [START_REF] Younès | From mixed finite elements to finite volumes for elliptic PDEs in two and three dimensions[END_REF]).

Our purpose is to build a discrete gradient with a local computation on the mesh interfaces, that is conformal in H(div). Our paradigm is to define this discrete gradient only using the scalar field and without considering the source term. On the contrary of the previously discussed works [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF][START_REF] Arbogast | Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences[END_REF][START_REF] Aavatsmark | Convergence of a symmetric mpfa method on quadrilateral grids[END_REF][START_REF] Wheeler | A multipoint flux mixed finite element method[END_REF], the expression of that discrete gradient will not be obtained through an approximation of a discrete mixed problem using quadrature rules. It will be obtained from the variational setting itself. The main idea is to choose a test function space that is L 2 -orthogonal with the shape functions, i.e. in duality with the Raviart-Thomas space. With a Petrov-Galerkin approach the spaces of the shape and test functions are different. It is now possible to insert duality between the shape and test functions and then to recover a local definition of the discrete gradient, as we proposed previously in the one-dimensional case [START_REF] Dubois | Finite volumes and mixed Petrov-Galerkin finite elements: the unidimensional problem[END_REF]. The stability analysis of the mixed finite element method emphasises the "inf-sup" condition [START_REF] Ladyzhenskaya | The Mathematical Theory of Viscous Incompressible Flow[END_REF][START_REF] Babuška | Error-bounds for finite element method[END_REF][START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers[END_REF]. In his fundamental contribution, Babuška [START_REF] Babuška | Error-bounds for finite element method[END_REF] gives general inf-sup conditions for mixed Petrov-Galerkin (introduced in [START_REF] Petrov | Application of Galerkin's method to the problem of stability of flow of a viscous fluid[END_REF]) formulation. The inf-sup condition guides the construction of the dual space.

In this contribution we extend the Petrov-Galerkin formulation to two-dimensional space dimension with Raviart-Thomas shape functions. In section 1, we introduce notations and general backgrounds. The discrete gradient is presented in section 2. Dual Raviart-Thomas test functions for the Petrov-Galerkin formulation of Poisson equation are proposed in section 3. In section 4, we retrieve the four point finite volume scheme of Herbin [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh[END_REF] for a specific choice of the dual test functions. Section 5 is devoted to the stability and convergence analysis in Sobolev spaces with standard finite element methods.

Background and notations

In the sequel, Ω ⊂ R 2 is an open bounded convex with a polygonal boundary. The spaces L 2 (Ω), H 1 0 (Ω) and H(div, Ω) are considered, see e.g. [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]. The L 2 -scalar products on L 2 (Ω) and on [L 2 (Ω)]

2 are similarly denoted (•, •) 0 .

Meshes

A conformal triangle mesh T of Ω is considered, in the sense of Ciarlet in [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. The angle, vertex, edge and triangle sets of T are respectively denoted T -1 , T 0 , T 1 and T 2 . The area of K ∈ T 2 and the length of a ∈ T 1 are denoted |K| and |a|.

Let K ∈ T 2 . Its three edges, vertexes and angles are respectively denoted a K,i , W K,i and θ K,i , (for 1 ≤ i ≤ 3) in such a way that W K,i and θ K,i are opposite to a K,i (see figure 1). The unit normal to a K,i pointing outwards K is denoted n K,i . The local scalar products on K are introduced as, for

 K θ K,i n K,i W K,i a K,i Figure 1 -Mesh notations for a triangle K ∈ T 2 θ a,K n a K L a a n a K N a S a W a E a ∂Ω a ∈ T 1 i , ∂ c a = (K, L) W a θ a,L a ∈ T 1 b , ∂ c a = (K)
f i ∈ L 2 (Ω) or p i ∈ [L 2 (Ω)] 2 : (f 1 , f 2 ) 0,K = K f 1 f 2 dx or (p 1 , p 2 ) 0,K = K p 1 • p 2 dx .
Let a ∈ T 1 . One of its two unit normal is chosen and denoted n a . This sets an orientation for a. Let S a , N a be the two vertexes of a, ordered so that (n a , S a N a ) has a direct orientation. The sets T 1 i and T 1 b of the internal and boundary edges respectively are defined as,

T 1 b = a ∈ T 1 , a ⊂ ∂Ω , T 1 i = T 1 \T 1 b . Let a ∈ T 1 i .
Its coboundary ∂ c a is made of the unique ordered pair K, L ∈ T 2 so that a ⊂ ∂K ∩ ∂L and so that n a points from K towards L. In such a case the following notation will be used:

a ∈ T 1 i , ∂ c a = (K, L)
 and we will denote W a (resp. E a ) the vertex of K (resp. L) opposite to a (see figure 2). Let a ∈ T 1 b : n a is assumed to point towards the outside of Ω. Its coboundary is made of a single K ∈ T 2 so that a ⊂ ∂K, which situation is denoted as follows:

a ∈ T 1 b , ∂ c a = (K)
and we will denote W a the vertex of K opposite to a. If a ∈ T 1 is an edge of K ∈ T 2 , the angle of K opposite to a is denoted θ a,K .

Finite element spaces

Relatively to a mesh T are defined the spaces P 0 and RT . The space of piecewise constant functions on the mesh is denoted by P 0 subspace of L 2 (Ω). The classical basis of P 0 is made of the indicators 1l K for K ∈ T 2 . To u ∈ P 0 is associated the vector (u K ) K∈T 2 so that u = K∈T 2 u K 1l K . The space of Raviart-Thomas of order 0 introduced in [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF] is denoted by RT and is a subspace of H(div, Ω). It is recalled that p ∈ RT if and only if p ∈ H(div, Ω) and for all K ∈ T 2 , p(x

) = α K + β K x, for x ∈ K, where α K ∈ R 2 and β K ∈ R are two constants. An element p ∈ RT is uniquely determined by its fluxes p a := a p • n a ds for a ∈ T 1 . The classical basis {ϕ a , a ∈ T 1 } of RT is so that b ϕ a • n b ds = δ ab
for all b ∈ T 1 and with δ ab the Kronecker symbol. Then to p ∈ RT is associated its flux vector (p a ) a∈T 1 so that, p = a∈T 1 p a ϕ a . The local Raviart-Thomas basis functions are defined, for K ∈ T 2 and i = 1, 2, 3, by:

ϕ K,i (x) = 1 4|K| ∇|x -W K,i | 2 on K and ϕ K,i = 0 otherwise. ( 1
)
With that definition:

ϕ a = ϕ K,i -ϕ L,j if a ∈ T 1 i , ∂ c a = (K, L) and a = a K,i = a L,j ϕ a = ϕ K,i if a ∈ T 1 b , ∂ c a = (K) and a = a K,i . (2)
The support of the RT basis functions is supp(ϕ

a ) = K ∪ L if a ∈ T 1 i , ∂ c a = (K, L) or supp(ϕ a ) = K in case a ∈ T 1 b , ∂ c a = (K)
. This provides a second way to decompose p ∈ RT as,

p = K∈T 2 3 i=1 p K,i ϕ K,i , where p K,i = εp a if a = a K,i with ε = n a • n K,i = ±1.
For simplicity we will denote ϕ K,a = ϕ K,i for a ∈ T 1 such that a ⊂ ∂K and a = a K,i . The divergence operator div : RT → P 0 is given by,

div p = K∈T 2 (div p) K 1l K , (div p) K = 1 |K| 3 1=1 p K,i . (3)



The two unbounded operators,

∇ : L 2 (Ω) ⊃ H 1 0 (Ω) → [L 2 (Ω)] 2 and div : [L 2 (Ω)] 2 ⊃ H(div, Ω) → L 2 (Ω) together satisfy the Green formula: for u ∈ H 1 0 (Ω) and p ∈ H(div, Ω): (∇u, p) 0 = -(u, div p) 0 . Identifying L 2 (Ω) and [L 2 (Ω)]
2 with their topological dual spaces using the L 2 -scalar product yields the following property,

∇ = -div ,
that is a weak definition of the gradient on H 1 0 (Ω). Consider a mesh of the domain and the associated spaces P 0 and RT as defined in section 1. We want to define a discrete gradient: ∇ T : P 0 → RT , based on a similar weak formulation. Starting from the divergence operator div : RT → P 0 , one can define div : (P 0 ) → (RT ) , between the algebraic dual spaces of P 0 and RT . The classical basis for P 0 is orthogonal for the L 2 -scalar product. Thus, P 0 is identified with its algebraic dual (P 0 ) . On the contrary, the Raviart-Thomas basis {ϕ a , a ∈ T 1 } of RT is not orthogonal. For this reason, a general identification process of (RT ) to a space RT = span (ϕ a , a ∈ T 1 ) is studied. We want it to satisfy,

ϕ a ∈ H(div, Ω), (ϕ a , ϕ a ) 0 = 0, (4)
so that RT ⊂ H(div, Ω), together with the orthogonality property,

(ϕ a , ϕ b ) 0 = 0 for a, b ∈ T 1 , a = b. (5)
The discrete gradient is defined with the diagram,

RT div ---→ P 0 Π     id RT ← --- div P 0 , ∇ T = -Π -1 • div : P 0 → RT, (6)
where Π : RT → RT is the projection defined by Πϕ a = ϕ a for any a ∈ T 1 .

Various choices for RT are possible. The first choice is to set RT = RT , and therefore to build {ϕ a , a ∈ T 1 } with a Gram-Schmidt orthogonalisation process on the Raviart-Thomas basis. Such a choice has an important drawback. The dual base function ϕ a does not conserve a support located around the edge a. The discrete gradient matrix will be a full matrix related with the Raviart-Thomas mass matrix inverse. This is not relevant with the definition of the original gradient operator that is local in space. This choice corresponds to the classical mixed finite element discrete gradient that is known to be associated with a full matrix [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]. In order to overcome this problem, Baranger, Maitre and Oudin [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF] have proposed to lump the mass matrix of the mixed finite element method. They obtain a discrete local gradient. Other methods have been proposed by  Thomas-Trujillo [START_REF] Thomas | Mixed finite volume methods[END_REF][START_REF] Thomas | Finite volume methods for elliptic problems: convergence on unstructured meshes[END_REF], by Noh [START_REF] Noh | CEL: A time-dependent, two-space-dimensional, coupled Euler-Lagrange code[END_REF], and analysed by Coudière, Vila and Villedieu [START_REF] Coudière | Convergence rate of a finite volume scheme for a two-dimensional convection-diffusion problem[END_REF]. Another approach is to add unknowns at the vertices, as developed by Hermeline [START_REF] Hermeline | Une méthode de volumes finis pour les équations elliptiques du second ordre[END_REF] and Domelevo-Omnes [START_REF] Domelevo | A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids[END_REF].

A second choice, initially proposed by Dubois and co-workers [START_REF] Dubois | Finite volumes and mixed Petrov-Galerkin finite elements: the unidimensional problem[END_REF][START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF][START_REF] Borel | Boundary conditions for Petrov-Galerkin finite volumes[END_REF][START_REF] Dubois | Dual Raviart-Thomas mixed finite elements[END_REF], is investigated in this paper. The goal is to search for a dual basis satisfying equation ( 4) and in addition to the orthogonality property (5), the localisation constraint,

∀ a ∈ T 1 , supp(ϕ a ) ⊂ supp(ϕ a ), (7)
in order to impose locality to the discrete gradient. We observe that due to the H(div)conformity, we have continuity of the normal component on the boundary of the coboundary of the edge a:

ϕ a • n b = 0 if a = b ∈ T 1 . (8)
With such a constraint [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF] the discrete gradient of u ∈ P 0 will be defined on each edge a ∈ T 1 only from the two values of u on each side of a (as detailed in proposition 1). In this context it is no longer asked to have ϕ a ∈ RT so that RT = RT : thus, this is a Petrov-Galerkin discrete formalism.

3 Raviart-Thomas dual basis Definition 1. (ϕ a ) a∈T 1 is said to be a Raviart-Thomas dual basis if it satisfies (4), the orthogonality condition (5), the localisation condition [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF] and the following flux normalisation condition:

∀ a, b ∈ T 1 , b ϕ a • n b ds = δ ab , (9)
as for the Raviart-Thomas basis functions ϕ a , see section 1.

In such a case, RT = span(ϕ a , a ∈ T 1 ) is the associated Raviart-Thomas dual space, Π : ϕ a ∈ RT → ϕ a ∈ RT the projection onto RT and ∇ T = -Π -1 div : P 0 → RT the associated discrete gradient, as described in diagram (6).

Computation of the discrete gradient

Proposition 1. Let (ϕ a ) a∈T 1 be a Raviart-Thomas dual basis. The discrete gradient is given for u ∈ P 0 , by the relation

∇ T u = a∈T 1 p a ϕ a with, if a ∈ T 1 i , ∂ c a = (K, L), p a = u L -u K (ϕ a , ϕ a ) 0 if a ∈ T 1 b , ∂ c a = (K), p a = -u K (ϕ a , ϕ a ) 0 . ( 10
)
The formulation of the discrete gradient only depends on the coefficients (ϕ a , ϕ a ) 0 . The discretisation of the Poisson equation (see the next subsection) also only depends on these coefficients.



The result of the localisation condition ( 7) is, as expected, a local discrete gradient: its value on an edge a ∈ T 1 only depends on the values of the scalar function u on each sides of a. The discrete gradient on the external edges expresses a homogeneous Dirichlet boundary condition. At the continuous level, the gradient defined on the domain H 1 0 (Ω) is the adjoint of the divergence operator on the domain H(div, Ω). That property is implicitly recovered at the discrete level. This is consistent since the discrete gradient is the adjoint of the divergence on the domain RT .

Proof. Condition (9) leads to b ϕ a • n b ds = b ϕ a • n b ds, for any a, b ∈ T 1 . Then the divergence theorem implies that ∀ p ∈ RT, ∀ K ∈ T 2 , K div p dx = K div(Πp) dx ,
and so proves

∀ (u, p) ∈ P 0 × RT, (div p, u) 0 = (div(Πp), u) 0 . (11)
Let us prove that,

∀ u ∈ P 0 , ∀ q ∈ RT , (∇ T u, q) 0 = -(u, div q) 0 . (12)
From property (5) one can check that,

∀ q 1 , q 2 ∈ RT , (Π -1 q 1 , q 2 ) 0 = (q 1 , Π -1 q 2 ) 0 .
Now consider u ∈ P 0 and q ∈ RT . We have with [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF], (u, div q) 0 = (u, div(Π -1 q)) 0 = (div u, Π -1 q) 0 = (Π -1 div u , q) 0 , which gives [START_REF] Coudière | Convergence rate of a finite volume scheme for a two-dimensional convection-diffusion problem[END_REF] by definition of the discrete gradient.

We can now prove [START_REF] Chavent | On the finite volume reformulation of the mixed finite element method for elliptic and parabolic PDE on triangles[END_REF]. Let u ∈ P 0 and p = ∇ T u ∈ RT that we decompose as

∇ T u = a∈T 1
p a ϕ a . For any a ∈ T 1 , with (5), ∇ T u, ϕ a 0 = p a ϕ a , ϕ a 0 , and meanwhile with equation ( 12) and ( 11) successively,

∇ T u, ϕ a 0 = -u, div ϕ a 0 = -u, div ϕ a 0 .
Finally, div ϕ a is explicitly given by,

if a ∈ T 1 i , ∂ c a = (K, L) : div ϕ a = 1 |K| 1l K - 1 |L| 1l L , ( 13
) if a ∈ T 1 b , ∂ c a = (K) : div ϕ a = 1 |K| 1l K .
This yields relations [START_REF] Chavent | On the finite volume reformulation of the mixed finite element method for elliptic and parabolic PDE on triangles[END_REF]. 

Petrov-Galerkin discretisation of the Poisson problem

Consider the following Poisson problem on Ω,

-∆u = f ∈ L 2 (Ω), u = 0 on ∂Ω. ( 14
)
Consider a mesh T and a Raviart-Thomas dual basis (ϕ a ) a∈T 1 as in definition 1 leading to the space RT . Let us denote V = P 0 × RT and V = P 0 × RT . The mixed Petrov-Galerkin discretisation of equation ( 14) is: find (u, p) ∈ V so that,

∀ (v, q) ∈ V , (p, q) 0 + (u, div q) 0 = 0 and -(div p, v) 0 = (f, v) 0 . (15)
The mixed Petrov-Galerkin discrete problem [START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF] reformulates as: find (u, p) ∈ V so that,

∀ (v, q) ∈ V , Z (u, p), (v, q) = -(f, v) 0 .
where the bilinear form Z is defined for (u, p) ∈ V and (v, q) ∈ V by,

Z (u, p), (v, q) ≡ (u, div q) 0 + (p, q) 0 + (div p, v) 0 . ( 16
)
Proposition 2 (Solution of the mixed discrete problem). The pair (u, p) ∈ V is a solution of problem [START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF] if and only if

∇ T u = p, -div(∇ T u) = f T , (17)
where f T ∈ P 0 is the projection of f , defined by,

f T = K∈T 2 f K 1l K , f K = 1 |K| K f dx.
If (ϕ a , ϕ a ) > 0 for all a ∈ T 1 , then problem (15) has a unique solution.

Proposition 2 shows an equivalence between the mixed Petrov-Galerkin discrete problem [START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF] and the discrete problem [START_REF] Eymard | Finite volume methods. Handb. Numer. Anal[END_REF]. Problem [START_REF] Eymard | Finite volume methods. Handb. Numer. Anal[END_REF] actually is a finite volume problem. Precisely, with [START_REF] Chavent | On the finite volume reformulation of the mixed finite element method for elliptic and parabolic PDE on triangles[END_REF], it becomes: find u ∈ P 0 so that, for all K ∈ T 2 :

a∈T 1 i , ∂ c a=(K,L) u L -u K (ϕ a , ϕ a ) 0 + a∈T 1 b , ∂ c a=(K) -u K (ϕ a , ϕ a ) 0 = |K|f K .
It is interesting to notice that this problem only involves the coefficients (ϕ a , ϕ a ) 0 that are going to be computed later.

Proof. Let u ∈ P 0 , denote p = ∇ T u ∈ RT and assume that div p = f T . Then using relation [START_REF] Coudière | Convergence rate of a finite volume scheme for a two-dimensional convection-diffusion problem[END_REF], equation [START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF] clearly holds. Conversely, consider (u, p) ∈ V a solution of problem [START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF]. Relation [START_REF] Coudière | Convergence rate of a finite volume scheme for a two-dimensional convection-diffusion problem[END_REF] implies that p = ∇ T u, as a result, -div(∇ T u) = f T . We assume that (ϕ a , ϕ a ) > 0 for all a ∈ T 1 and prove existence and uniqueness.

It suffices to prove that u = 0 is the unique solution when f T = 0. In such a case, div(∇ T u) = 0, and using successively [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] and ( 12):

0 = -(div(∇ T u), u) 0 = -(div(Π ∇ T u), u) 0 = (Π ∇ T u, ∇ T u) 0 = a∈T 1 p 2 a (ϕ a , ϕ a ).
As a result p a = 0 for all a ∈ T 1 and p = ∇ T u = 0. From ( 15) it follows that for all q ∈ RT we have (u, div q) 0 = 0. Thus with [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] we also have (u, div q) 0 = 0 for all q ∈ RT . Since div(RT ) = P 0 it follows that u = 0.

Retrieving the four point finite volume scheme

In this section we present sufficient conditions for the construction of Raviart-Thomas dual basis. These conditions will allow to compute the coefficients (ϕ a , ϕ a ) 0 . We start by introducing the normal flux g on the edges, and the divergence of the dual basis δ K on K ∈ T 2 . Let g : (0, 1) → R be a continuous function so that, 1 0 g ds = 1, 1 0 g(s)s 2 ds = 0, g(0) = 0 and g(s) = g(1 -s).

(

) 18
On a mesh T are defined g K,i : a K,i → R for K ∈ T 2 and i = 1, 2, 3 as,

g K,i (x) = g(s)/|a K,i | for x = sS K,i + (1 -s)N K,i . (19)
For K ∈ T 2 is denoted δ K ∈ L 2 (K) a function that satisfies,

K δ K dx = 1 and K δ K (x)|x -W K,i | 2 dx = 0 for i = 1, 2, 3. ( 20
)
To a family (ϕ K,i ) of functions on Ω for K ∈ T 2 and for i = 1, 2, 3 is associated the family (ϕ a ) a∈T 1 so that,

ϕ a = ϕ K,i -ϕ L,j if a ∈ T 1 i , ∂ c a = (K, L) and a = a K,i = a L,j ϕ a = ϕ K,i if a ∈ T 1 b , ∂ c a = (K) and a = a K,i . (21)
This is the same correspondence as in (2) between the Raviart-Thomas local basis functions (ϕ K,i ) and the Raviart-Thomas basis functions (ϕ a ) a∈T 1 . Similarly, we will denote ϕ K,a = ϕ K,i for a ∈ T 1 such that a ⊂ K and a = a K,i .

Theorem 1. Consider a family (ϕ K,i ) K∈T 2 , i=1, 2, 3 of local basis functions on Ω that satisfy supp ϕ K,i ⊂ K [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh[END_REF] and independently on i,

div ϕ K,i = δ K , on K. ( 23
)
On ∂K, the normal component is given by

ϕ K,i • n K = g K,i on a K,i 0 otherwise , (24)
where g K,i and δ K satisfy equations ( 18), ( 19) and [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF]. Let (ϕ a ) a∈T 1 be constructed from the local basis functions (ϕ K,i ) K,i with equation [START_REF] Godounov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF]. Then (ϕ a ) a∈T 1 is a Raviart-Thomas dual basis as defined in definition 1. Moreover, the coefficients (ϕ a , ϕ a ) 0 only depend on the mesh T geometry,

a ∈ T 1 i , ∂ c a = (K, L) ⇒ (ϕ a , ϕ a ) 0 = (cotan θ a,K + cotan θ a,L ) /2, a ∈ T 1 b , ∂ c a = (K) ⇒ (ϕ a , ϕ a ) 0 = cotan θ a,K /2 . ( 25
)
Notations are recalled on figure 3. We will also denote g a,K = g K,i for a ∈ T 1 such that a ⊂ K and a = a K,i . Corollary 1. Assume that the mesh satisfies the Delaunay condition: for all internal edge a ∈ T 1 we have the angle condition θ a,K + θ a,L < π (denoting ∂ c a = (K, L)). Also assume that for any boundary edge a, θ a,K < π/2 (denoting ∂ c a = (K)). Then with ( 25), (ϕ a , ϕ a ) 0 > 0 and proposition 2 ensures the existence and uniqueness of the solution to the discrete problem. Moreover, the mixed Petrov-Galerkin discrete problem [START_REF] Eymard | Finite volume methods. Handb. Numer. Anal[END_REF] for the Laplace equation ( 14) coincides with the four point finite volume scheme defined and analysed in Herbin [START_REF] Herbin | An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh[END_REF].

Therefore, the Raviart-Thomas dual basis does not need to be constructed. Whatever are δ K and g that satisfy equations ( 18), ( 19) and ( 20), the coefficients (ϕ a , ϕ a ) 0 will be unchanged. They only depend on the mesh geometry and are given by equation [START_REF] Ladyzhenskaya | The Mathematical Theory of Viscous Incompressible Flow[END_REF]. Practically, this means that neither the (ϕ a ) a∈T 1 nor δ K and g need to be computed. Such a dual basis will be explicitly computed in section 5.1. The numerical scheme will always coincide with the four point volume scheme. Finally, this theorem provides a new point of view for the understanding and analysis of finite volume methods. Theorem 1 gives sufficient conditions in order to build Raviart-Thomas dual basis. In the sequel we will focus on such Raviart-Thomas dual basis, though more general ones may exist: this will not be discussed in this paper.

Proof of corollary 1. We have the general formula cotan θ 1 +cotan θ 2 = sin(θ 1 +θ 2 )/(sin θ 1 sin θ 2 ) that ensures that (ϕ a , ϕ a ) 0 > 0 under the assumptions in the corollary. For the equivalence between the two schemes, it suffices to prove that cotan θ a,K /2 = d a,K /|a| where d a,K denotes the distance between the edge a and the circumcircle centre C of K. Denote S and N the two vertexes of a. Then the angle SCN = 2θ a,K . The distance d a,K is equal to CH with H the orthogonal projection of C on a. The triangle SCN being isosceles, H is also the midle of [SN ]. In the right angled triangle SCH we have SCH = SCN /2 = θ a,K and cotan SCH = CH/SH = d a,K /(|a|/2) which gives the result.

Proof of theorem 1. Consider as in theorem 1 a family (ϕ K,i ) K∈T 2 , i=1, 2, 3 that satisfy, ( 22), ( 23) and [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF] for δ K and g K,i such that the assumptions [START_REF] Faille | A control volume method to solve an elliptic equation on a two-dimensional irregular mesh[END_REF], [START_REF] Faille | Des mathématiciens découvrent les volumes finis[END_REF] and [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF] are true. Let (ϕ a ) a∈T 1 be constructed from the local basis functions (ϕ K,i ) K,i with equation [START_REF] Godounov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF].

Let us first prove that (ϕ a ) a∈T 1 is a Raviart-Thomas dual basis as in definition 1. Consider an internal edge a ∈ T 1 , a = (K|L). With [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF], we have supp ϕ a = K ∪ L and relation [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF] holds. With [START_REF] Godounov | Résolution numérique des problèmes multidimensionnels de la dynamique des gaz[END_REF],

ϕ a |K = ϕ K,a ∈ H(div, K), ϕ a |L = -ϕ L,a ∈ H(div, L).
The normal flux ϕ a • n a is continuous across a = K ∩ L since g K,a = g L,a and with [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF]. Moreover, ϕ a • n = 0 on the boundary of K ∪ L due to [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF]. Therefore ϕ a belongs to H(div, Ω). With formula [START_REF] Ladyzhenskaya | The Mathematical Theory of Viscous Incompressible Flow[END_REF] and the angle condition made in theorem 1, (ϕ a , ϕ a ) 0 = 0 and so ( 4 23), [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF] and the fact that W K,b is opposite to b and so is a vertex of a. This implies the orthogonality condition [START_REF] Arbogast | Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences[END_REF] with the assumptions in [START_REF] Faille | A control volume method to solve an elliptic equation on a two-dimensional irregular mesh[END_REF] and [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF].

±4|K|(ϕ a , ϕ b ) 0 = - K div ϕ a |x -W K,b | 2 dx + ∂K ϕ a • n|s -W K,b | 2 ds = - K δ K |x -W K,b | 2 dx + 1 0 g(s) s 2 ds, using (
It remains to prove [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers[END_REF]. In the case where a, b ∈ T 1 are two distinct edges, b ϕ a •n b ds = 0. Assume that a ∈ T 1 is an edge of K ∈ T 2 . We have n a = εn K,a with ε = ±1. With relation [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF] and the divergence formula,

a ϕ a • n a ds = a (εϕ K,a ) • (εn K,a ) ds = ∂K ϕ K,a • n ds = K div ϕ K,a dx.
This ensures that a ϕ a • n a ds = 1 with relation [START_REF] Hermeline | Une méthode de volumes finis pour les équations elliptiques du second ordre[END_REF] and the first assumption in [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF]. We successively proved ( 4), ( 5), ( 7) and ( 9) and then (ϕ a ) a∈T 1 is a Raviart-Thomas dual basis.

Let us now prove [START_REF] Ladyzhenskaya | The Mathematical Theory of Viscous Incompressible Flow[END_REF]. Let a ∈ T 1 an internal edge with the notations in figure 3. The Raviart-Thomas basis function ϕ a has its support in K ∪ L, so that

(ϕ a , ϕ a ) 0 = K ϕ a • ϕ a dx + L ϕ a • ϕ a dx.
With the local decompositions ( 2) and ( 21) we have,

(ϕ a , ϕ a ) 0 = K ϕ K,a • ϕ K,a dx + L ϕ L,a • ϕ L,a dx .
By relation [START_REF] Aavatsmark | An introduction to multipoint flux approximations for quadrilateral grids[END_REF], W being the opposite vertex to the edge a in the triangle K,

4|K| K ϕ K,a • ϕ K,a dx = K ϕ K,a ∇|x -W | 2 dx = - K div ϕ K,a |x -W | 2 dx + ∂K ϕ K,a • n K |x -W | 2 dσ.
By hypothesis [START_REF] Hermeline | Une méthode de volumes finis pour les équations elliptiques du second ordre[END_REF] and [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF], and using [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF],

4|K| K ϕ K,a • ϕ K,a dx = K δ K |x -W | 2 dx + a g K,a |x -W | 2 dσ = a g K,a |x -W | 2 dσ.
Let H be the orthogonal projection of the point W on the edge a. We have |x -W | 2 = W H 2 + |x -H| 2 and with ( 18) and ( 19), a g K,a dσ = |a| 1 0 g(s)/|a| ds = 1 and so, 4|K|

K ϕ K,a • ϕ K,a dx = W H 2 + a g K,a |x -H| 2 dσ.
Let s and s respectively be the curvilinear coordinates of x and H on a with origin S, then

4|K| K ϕ K,a • ϕ K,a dx = W H 2 + |a| 2 1 0 (s -s) 2 g(s)ds.
The assumptions in [START_REF] Faille | A control volume method to solve an elliptic equation on a two-dimensional irregular mesh[END_REF] on g imply that 2

1 0 g(s)sds = 1. By expanding (s -s) 2 = s 2 -2ss + s 2 we get, 1 0 (s -s) 2 g(s)ds = s 2 -s . It follows that, 4|K| K ϕ K,a • ϕ K,a dx = W H 2 + |a|s |a|(s -1) = W H 2 + -→ SH • --→ N H = --→ W S • --→ W N .



Raviart-Thomas finite elements of Petrov-Galerkin type Some trigonometry results in K leads to sin θ K,a = 2|K| W S•W N . As a result,

4|K| K ϕ K,a • ϕ K,a dx = 2|K| cotan θ K,a ,
this gives [START_REF] Ladyzhenskaya | The Mathematical Theory of Viscous Incompressible Flow[END_REF].

Stability and convergence

In this section we develop a specific choice of dual basis functions. We provide for that choice technical estimates and prove a theorem of stability and convergence. With theorem 1, this leads to an error estimate for the four point finite volume scheme. We begin with the main result in theorem 2. Theorem 3 provides a methodology in order to get the inf-sup stability conditions. The inf-sup conditions need technical results that are proved in subsections 5.1 to 5.2. We will need the following angle condition. Angle assumption. Let θ and θ chosen such that

0 < θ < θ < π/2 (26)
We consider meshes T such that all the angles of the mesh are bounded from below and above by θ and θ respectively:

∀ θ ∈ T -1 , θ ≤ θ ≤ θ . (27)
With that angle condition, the coefficients (ϕ a , ϕ a ) in ( 25) are strictly positive. With proposition 2 this ensures the existence and uniqueness for the solution (u T , p T ) of the mixed Petrov-Galerkin discrete problem [START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF].

Theorem 2 (Error estimates). We assume that Ω ⊂ R 2 is a bounded polygonal convex domain and that f ∈ H 1 (Ω). Under the angle hypotheses ( 26) and ( 27), there exists a constant C independent on T satisfying [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF] and independent on f so that the solution (u T , p T ) of the mixed Petrov-Galerkin discrete problem (15) satisfies,

u T 0 + p T H(div,Ω) ≤ C f 0 .
Let u be the exact solution to problem ( 14) and p = ∇u the gradient, the following error estimates holds,

u -u T 0 + p -p T H(div,Ω) ≤ Ch T f 1 , (28)
with h T the maximal size of the edges of the mesh.

Proof. We prove that the unique solution of the mixed Petrov-Galerkin (15) continuously depends on the data f . The bilinear form Z defined in ( 16) is continuous, with a continuity constant M independent on the mesh T ,

|Z(ξ, η)| ≤ M ξ L 2 ×H div η L 2 ×H div , ∀ ξ ∈ V, η ∈ V .
The following uniform inf-sup stability condition: there exists a constant β > 0 independent on T such that,

∀ ξ ∈ V, so that ξ L 2 ×H div = 1, ∃ η ∈ V , η L 2 ×H div ≤ 1 and Z(ξ, η) ≥ β, (29)
is proven in theorem 3 under some conditions. Moreover, the two spaces V and V have the same dimension. Then the Babuška theorem in [START_REF] Babuška | Error-bounds for finite element method[END_REF], also valid for Petrov-Galerkin mixed formulation, applies. The unique solution ξ T = (u T , p T ) of the discrete scheme (15) satisfies the error estimates, and

ξ -ξ T L 2 ×H div ≤ 1 + M β inf ζ∈V ξ -ζ L 2 ×H div ,
with ξ = (u, p), u the exact solution to the Poisson problem ( 14) and p = ∇u. In our case, this formulation is equivalent to

u -u T 0 + p -p T H(div,Ω) ≤ C inf v∈P 0 u -v 0 + inf q∈RT p -q H(div,Ω) (30)
for a constant C = 1 + M β dependent of T only through the lowest and the highest angles θ and θ . With the interpolation operators Π 0 : L 2 (Ω) → P 0 and Π RT :

H 1 (Ω) 2 → RT 0 u -u T 0 + p -p T H(div,Ω) ≤ C u -Π 0 u 0 + p -Π RT p H(div,Ω) .
On the other hand we have the following interpolation errors:

u-Π 0 u 0 ≤ C 1 h T u 1 , p-Π RT p 0 ≤ C 2 h T p 1 , div p-Π RT p 0 ≤ C 1 h T div p 1 .
On the left, we have the Poincaré-Wirtinger inequality where the constant C 1 = 1/π is independent on the mesh, due to [START_REF] Payne | An optimal Poincaré inequality for convex domains[END_REF]. The third inequality is the same as the first one since Π 0 div p = div Π RT p. For the second inequality, the constant C 2 has been proven in [START_REF] Acosta | The maximum angle condition for mixed and nonconforming elements: Application to the stokes equations[END_REF] to be dominated by 1/ sin θ with θ the maximal angle of the mesh. Then ,

u -u T 0 + p -p T H(div,Ω) ≤ C h T u 1 + p 1 + div p 0 ,
with a constant C only depending on the maximal angle θ . Since -∆u = f in Ω, with f ∈ H 1 (Ω) and Ω convex, then u ∈ H 2 (Ω) and u 2 ≤ c f 0 . Moreover p = ∇u and div p = -f leads to

u -u T 0 + p -p T H(div,Ω) ≤ C h T 2 f 0 + f 1 .
Finally, we get

u -u T 0 + p -p T H(div,Ω) ≤ C h T f 1 ,
that is exactly [START_REF] Payne | An optimal Poincaré inequality for convex domains[END_REF].

Theorem 3 (Abstract stability conditions). Assume that the projection Π : RT → RT , such that Πϕ a = ϕ a in diagram (6) satisfies, for any p ∈ RT :

(p, Πp) 0 ≥ A p 2 0 , (H1)

Πp 0 ≤ B p 0 , (H2)
(div p, div Πp) 0 ≥ C div p 2 0 , (H3)
div Πp 0 ≤ D div p 0 (H4)

where A, B, C, D > 0 are constants independent on T . Then the uniform discrete inf-sup condition [START_REF] Pert | Physical constraints in numerical calculations of diffusion[END_REF] holds: there exists a constant β > 0 independent on T such that,

∀ ξ ∈ V, so that ξ L 2 ×H div = 1, ∃ η ∈ V , η L 2 ×H div ≤ 1 and Z(ξ, η) ≥ β.
This result has been proposed by Dubois in [START_REF] Dubois | Petrov-Galerkin finite volumes[END_REF]. For the completeness of this contribution, the proof (presented in the preprint [START_REF] Dubois | Dual Raviart-Thomas mixed finite elements[END_REF]) is detailed in Annex A.

In order to prove the conditions (H1), (H2), (H3) and (H4), one needs some technical lemmas on some estimations of the dual basis functions so that theorem 3 holds. It is the goal of the next subsections.

A specific Raviart-Thomas dual basis

Choice of the divergence For K a given triangle of T 2 , we propose a choice for the divergence δ K of the dual basis functions ϕ K,i , 1 ≤ i ≤ 3 in [START_REF] Hermeline | Une méthode de volumes finis pour les équations elliptiques du second ordre[END_REF]. We know from [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF] that this function has to be L 2 (K)-orthogonal to the three following functions: |x -W K,i | 2 for i =1, 2, 3 and that its integral over K is equal to 1. We propose to choose δ K as the solution of the leastsquare problem: minimise K δ 2 K dx with the constraints in [START_REF] Fraeijs De Veubeke | Displacement and equilibrium models in the finite element method[END_REF]. It is well-known that the solution belongs to the four dimensional space

E K = span (1l K , |x -W K,i | 2 , 1 ≤ i ≤ 3)
and is obtained by the inversion of an appropriate Gram matrix. Lemma 1. For the above construction of δ K , we have the following estimation:

|K| K δ 2 K dx ≤ ν, with ν = 8 • 3 5 • 23 5 1 tan 4 θ .
The proof of this result is technical and has been obtained with the help of a formal calculus software. It is detailed in Annex C.

Choice of the flux on the boundary of the triangle

A continuous function g : (0, 1) → R satisfying the conditions (18) can be chosen as the following polynomial:

g(s) = 30s (s -1) (21s 2 -21s + 4). ( 31
)
Construction of the Raviart-Thomas dual basis For a triangle K and an edge a of K, we construct now a possible choice of the dual function ϕ K,a satisfying ( 22), ( 23) and [START_REF] Kershaw | Differencing of the diffusion equation in lagrangian hydrodynamic codes[END_REF]. Let F K,a be an affine function that maps the reference triangle K into the triangle K such that the edge â ≡ [0, 1]×{0} is transformed into the given edge a ⊂ ∂K. Then the mapping K x -→ x = F K,a ( x) ∈ K is one to one. We define x = F K,a (x) for any x ∈ K and the right hand side δ K (x) = 2 |K| δ K (x). With g defined in [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF], let us define g ∈ H 1/2 (∂ K) according to g := g on â = [0, 1] × {0} 0 elsewhere on ∂ K .

Since K δ K dx = 1 = ∂ K g dγ, the inhomogeneous Neumann problem ∆ζ K = δ K in K , ∂ζ K ∂n = g on ∂ K (32)
is well posed. The dual function ϕ K,a is defined according to

ϕ K,a (x) = 1 det(dF K,a ) dF K,a ∇ζ K . (33)
These so-defined functions satisfy the hypotheses ( 22), ( 23) and ( 24) of theorem 1. Let us now estimate their L 2 -norm.

L 2 -norm of the Raviart-Thomas dual basis

An upper bound on the L 2 norm of the Raviart-Thomas dual basis will be needed in order to prove the stability conditions in theorem 3. This bound is given in lemma 3. It only involves the mesh minimal angle θ . Lemma 2. For K ∈ T 2 and a ∈ T 1 , a ⊂ ∂K, we have ϕ K, a 0 K ≤ µ where µ is essentially a function of the smallest angle θ of the triangulation.

Proof. Since the reference triangle K is convex and g ∈ H 1/2 (∂ K), the solution ζ K of the Neumann problem [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF] satisfies the regularity property (see for example [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. 1[END_REF]) ζ K ∈ H 2 ( K), continuously to the data:

ζ K 2, K ≤ C K δ K 0, K + g 1/2, ∂ K .
Moreover thanks to lemma 1,

δ K 2 0, K = K δ K 2 d x = K (2|K|δ K ) 2 1 det(dF K,a ) dx = 2 |K| K δ 2 K dx ≤ 2 ν and then ∇ζ K 0, K ≤ C K √ 2ν + g 1/2, ∂ K .
Since the dual function ϕ K,a is defined by [START_REF] Rivas | Be03, programme de calcul tridimensionnel de la transmission de chaleur et de l'ablation[END_REF] and dF K,a 2 ≤ 8|K| sin θ from direct geometrical computations on the triangle K, we obtain

ϕ K,a 2 0,K ≤ 1 2 |K| 2 8 |K| sin θ ∇ζ K 2 0, K (2 |K|) . Then ϕ K,a 2
0,K ≤ (µ ) 2 , with (µ ) 2 = 4 sin θ C 2 K √ 2ν + g 1/2, ∂ K 2 .
Lemma 3. For K ∈ T 2 and q ∈ RT :

Πq 2 0,K ≤ 3(µ ) 2 3 i=1 q 2 K,i .
Proof. We have for a triangle K, Πq = 3 i=1 q K,i ϕ K,i , and so

Πq 2 0,K ≤ 3 i=1 |q K,i | ϕ K,i 0,K 2 ≤ 3 i=1 |q K,i | 2 3 i=1 ϕ K,i 2
0,K
Then lemma 2 applies and: Πq 2 0,K ≤ 3(µ ) 2 3 i=1 q 2 K,i .

Local Raviart-Thomas mass matrix

The proof of the stability conditions in theorem 3 involves lower and upper bounds of the eigenvalues of the local Raviart-Thomas mass matrix. We will need the following result proved in Annex B.

The hypotheses of theorem 3 are satisfied

Let us finally prove that the conditions (H1), (H2), (H3) and (H4) of theorem 3 hold. The proof relies on lemma 4, lemma 3 and lemma 1 involving the mesh independent constants λ , λ , µ and ν. In the following, p denotes an element of RT and K a fixed mesh triangle. It is recalled that on K, p =

3 i=1 p K,i ϕ K,i .
Condition (H1). Using the orthogonality property (5), and relation [START_REF] Ladyzhenskaya | The Mathematical Theory of Viscous Incompressible Flow[END_REF] successively, leads to

(Πp, p) 0,K = 3 i=1 p 2 K,i (ϕ K,i , ϕ K,i ) 0,K = 1 2 3 i=1 p 2 K,i cotan θ K,i ≥ 1 2 cotan θ 3 1=1 p 2 K,i .
Lemma 4 gives a lower bound, (Πp, p) 0,K ≥ cotan θ 2λ p 2 0,K . Summation over all K ∈ T 2 gives (H1) with,

A = cotan θ 2λ = 2 5 cotan θ tan θ .
Condition (H2). Using successively lemma 3 and lemma 4 we get,

Πp 2 0,K ≤ 3(µ ) 2 3 i=1 p 2 K,i ≤ 3(µ ) 2 λ p 2 0,K .
With the values of λ given in lemma 3 this implies (H2) with,

B = 3(µ ) 2 λ = 12 tan θ µ .
Condition (H3). Relation ( 11) induces (div Πp, div p) 0,K = div p 2 0,K since div p is a constant on K, and as a result inequality (H3) indeed is an equality with

C = 1.
Condition (H4). With equation (3) we get div p 2 0,K =

Conclusion

In this contribution we present a way to define a local discrete gradient of a piecewise constant function on a triangular mesh. This discrete gradient is obtained from a Petrov-Galerkin formulation and belongs to the Raviart-Thomas function space of low order. We have defined suitable dual test functions of the Raviart-Thomas basis functions. For the Poisson problem, we can interpret the Petrov-Galerkin formulation as a finite volume method. Specific constraints for the dual test functions enforce stability. Then the convergence can be established with the usual methods of mixed finite elements. It would be interesting to try to extend this work in several directions: the three-dimensional case, the case of general diffusion problems and also the case of higher degree finite element methods.

Proof. Let u ∈ P 0 be a discrete scalar function supposed to be constant in each triangle K of the mesh T . Let ψ ∈ H 1 0 (Ω) be the variational solution of the Poisson problem ∆ ψ = u in Ω , ψ = 0 on ∂Ω .

(

) 40
Since Ω is convex, the solution ψ of the problem (40) belongs to the space H 2 (Ω) and there exists some constant G > 0 that only depends on Ω such that

ψ 2 ≤ G u 0 .
Then the field ∇ψ belongs to the space H 1 (Ω) × H 1 (Ω). It is in consequence possible to interpolate this field in a continuous way (see e.g. Roberts and Thomas [START_REF] Roberts | Mixed and Hybrid Methods[END_REF]) in the space H(div, Ω) with the help of the fluxes on the edges:

p a = a ∂ψ ∂n a dγ , p = a ∈T 1 p a ϕ a ∈ RT .
Then there exists a constant L > 0 such that

p H div ≤ L u 0 . (41)
The two fields div p and u are constant in each element K of the mesh T . Moreover, we have: Then we have exactly, div p = u in Ω because this relation is a consequence of the above property for the mean values. Let now Π p be the interpolate of p in the "dual space" RT and q = 1 C Π p,

q = 1 C Π p = 1 C a∈T 1
p a ϕ a with Π p = a∈T 1 p a ϕ a .

We have as a consequence of (36) and div p = u that, ( u , div q ) 0 = 1 C ( div p , div Π p ) ≥ div p 2 0 = u 2 0 that establishes [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques[END_REF]. Moreover, we have due to equations ( 35), ( 37) and ( 41):

q 0 = 1 C Π p 0 ≤ B C p 0 ≤ BL C u 0 , div q 0 = 1 C div Π p 0 ≤ D C div p 0 = D C u 0 .
Then the two above inequalities establish the estimate [START_REF] Vohralík | Mixed finite element methods: implementation with one unknown per element, local flux expressions, positivity, polygonal meshes, and relations to other methods[END_REF] with With the properties (49) and (50), tr(G K ) = 15 4|K| ρ 2 K . This leads to the value of λ thanks to lemma 5.

F = 1 C √ B 2 L 2 + D
Derivation of λ . In order to compute λ , we want to find a lower bound for the smallest eigenvalues of the Gram matrix G K . The characteristic polynomial is given by

P (λ) = -det(λI -G K ) = -[λ 3 -tr(G K )λ 2 + Rλ -det G K ]
where R := 3 i=1 R i with R i := ϕ i 2 0 ϕ i+1 2 0 -(ϕ i , ϕ i+1 ) 2 0,K with the usual notation if i = 3, ϕ i+1 = ϕ 1 . Since P (λ) is of degree 3 with positive roots, the smallest root λ is such that λ ≥ det(G K ) R . As G K is a Gram matrix, the determinant of G K is the square of the volume of polytope generated by the basis function:

det(G K ) = vol(ϕ 1 , ϕ 2 , ϕ 3 ) 2 .
We expand each basis function on the orthogonal basis made of the three vector fields: -→ i , -→ j , x -G. Then the volume can be computed via a 3 by 3 elementary determinant. This leads to det(G K ) = ρ 2 K 16|K| .

The explicit computation of R i with help of (50) leads to

R i =

Figure 2 -

 2 Figure 2 -Mesh notations for an internal edge (left) and for a boundary edge (right)

Figure 3 -

 3 Figure 3 -Co-boundary of the edge a ∈ T 1 .

  ) holds. Consider two distinct edges a, b ∈ T 1 . If a and b are not two edges of a same triangle K ∈ T 2 , then ϕ a and ϕ b have distinct supports so that (ϕ a , ϕ b ) 0 = 0. If a and b are two edges of K ∈ T 2 , then (ϕ a , ϕ b ) 0 = K ϕ a • ϕ b dx. With the definition (1) of the local RT basis functions and using the Green formula,

Figure 4 -

 4 Figure 4 -Affine mapping F K,a between the reference triangle K and the given triangle K.

Lemma 4 .

 4 For p ∈ RT and K ∈ T 2 : for two constants λ and λ only depending on θ in (26),

3 1=1 p K,i 2 /

 32 |K| and with condition[START_REF] Hermeline | Une méthode de volumes finis pour les équations elliptiques du second ordre[END_REF], div Πp = δ K (x) 3 1=1 p K,i . Therefore we get,

2

 2 and the proposition is proven. on K. One can easily prove that, is the trace of G K .

Annex A Proof of theorem 3

In this section, we consider meshes T that satisfy the angle conditions [START_REF] Patankar | Numerical Heat Transfer and Fluid Flow[END_REF] parametrised by the pair 0 < θ < θ < π 2 . We suppose that the interpolation operator Π defined in section 1 by Π : RT -→ RT with Πϕ a = ϕ a satisfies the following properties: there exist four positive constants A, B, C and D only depending on θ and θ such that for all q ∈ RT (q, Πq) ≥ A q 2 0 ,

(div q, div Πq) 0 ≥ C div q 2 0 ,

div Πq 0 ≤ D div q 0 .

Let us first prove the following proposition relative to the lifting of scalar fields.

Proposition 3 (Divergence lifting of scalar fields). Under the previous hypotheses [START_REF] Roberts | Mixed and Hybrid Methods[END_REF], [START_REF] Thomas | Finite volume methods for elliptic problems: convergence on unstructured meshes[END_REF], [START_REF] Thomas | Mixed finite volume methods[END_REF] and [START_REF] Vohralík | Equivalence between lowest-order mixed finite element and multi-point finite volume methods on simplicial meshes[END_REF], there exists some strictly positive constant F that only depends of the minimal and maximal angles θ and θ such that for any mesh T and for any scalar field u constant in each element K of T , (u ∈ P 0 ), there exists some vector field q ∈ RT , such that

Proof of theorem 3

We suppose that the dual Raviart-Thomas basis satisfies the Hypothesis (34) to [START_REF] Vohralík | Equivalence between lowest-order mixed finite element and multi-point finite volume methods on simplicial meshes[END_REF]. We introduce the constant F > 0 such that [START_REF] Vohralík | Mixed finite element methods: implementation with one unknown per element, local flux expressions, positivity, polygonal meshes, and relations to other methods[END_REF] and ( 39) are realised for some vector field q ∈ RT for any u ∈ P 0 :

• We set a = 1 2

with the constants F , A, B and D introduced in ( 42), ( 34), ( 35) and ( 37) respectively. We shall prove that for

the inf-sup condition

Then we have after an elementary algebra: a F + a 2 = 1. In consequence,

thanks to the relations (43) and (45):

• Consider now ξ ≡ (u, p) satisfying the hypothesis of unity norm in the product space:

(48) Then at last one of these terms is not too small and due to the three terms that arise in relation (48), the proof is divided into three parts.

Then, div v 0 = 1 and

and the relation ( 44) is satisfied in this particular case.

(ii) If the conditions div p 0 ≤ β and p 0 ≥ α are satisfied, we set

We check that η L 2 ×H div ≤ 1:

Then

is exactly the inequality (47). Then the relation ( 44) is satisfied in this second case.

(iii) If the last conditions div p 0 ≤ β and p 0 ≤ α are satisfied, we first remark that the first component u has a norm bounded below: from (46),

with a discrete vector field q satisfying the inequalities [START_REF] Younès | A new formulation of the mixed finite element method for solving elliptic and parabolic PDE with triangular elements[END_REF]. Then,

because, due to (46) we have the following inequalities:

Then the relation (44) is satisfied in this third case and the proof is completed.

Annex B Proof of lemma 4

We first recall the statement of lemma 4.



Raviart-Thomas finite elements of Petrov-Galerkin type Lemma 4. For p ∈ RT and K ∈ T 2 :

for two constants λ and λ only depending on θ in [START_REF] Noh | CEL: A time-dependent, two-space-dimensional, coupled Euler-Lagrange code[END_REF],

The following technical result will be necessary for the proof of lemma 4.

Lemma 5. The gyration radius of a triangle K is defined as,

Proof. Let A i and a i , 1=1, 2, 3, be respectively the three vertices and edges of the triangle K. One can check that: 36

K , that gives the lower bound.

On the other hand, using the definition of the tangent,

, that gives the upper bound.

Proof of lemma 4

For a triangle K, the local RT mass matrix is G K := [(ϕ K,i , ϕ K,j ) 0,K ] 1≤i,j≤3 . Explicit computation obtained by Baranger-Maitre-Oudin in [START_REF] Baranger | Connection between finite volume and mixed finite element methods[END_REF] gives some properties on the gyration radius:

where θ i are the angles of the triangle K and lead to information on the Raviart-Thomas basis as follows:

where k is the third index of the triangle K (k = i, j, 1 ≤ i, j , k ≤ 3).

Derivation of λ . The triangle K ∈ T 2 is fixed and p ∈ RT rewrites p =

Annex C Proof of lemma 1

We express the function δ K as a linear combination of the functions 1l K and |x-W K,i | 2 , for 1 ≤ i ≤ 3. Thanks to the conditions (20), we solve formally a 4 by 4 linear system (with the help of a formal calculus software) in order to explicit the components. We can then compute the integral I given by,

The result is a symmetric function of the length |a i | of the three edges of the triangle K. It is a ratio of two homogeneous polynomials of degree 12. More precisely I reads,

where N and D respectively are homogeneous polynomials of degree 12 and 4. The exact expressions of D and N are,

with the following definitions,

and where σ p is the sum of of the three edges length |a j | to the power p:

The lemma 1 states an upper bound of I. To prove it, we look for an upper bound of N and a lower bound of D.

The denominator D in (51) is the difference of two positive expressions. We remark that,

We have on the one hand,

) and on the other hand a 2 i a 2 j ≤ 1 2 a 4 i + a 4 j . Then by summation Σ 2,2,0 ≤ σ 4 .

(54)

In the expression of D in (51), we split the term relative to σ 4 into two parts:

Then thanks to (53), , so ending the proof of lemma 1.