V Ruas

F A Radu

Hermite analogs of the lowest order Raviart-Thomas mixed method for convection-diffusion equations

 proved to be a valid alternative to attain higher order approximation of the primal variable, while keeping intact the matrix structure and the quality of the discrete flux variable of the original RT 0 method. Non trivial extensions of this method are studied here, that can be viewed as Hermite analogs of the two Douglas & Roberts' versions of the RT 0 method, to solve convection-diffusion equations. A detailed convergence study is carried out for one of the Hermite methods, and numerical results illustrate the performance of both of them, as compared to each other and to the corresponding mixed methods.

Introduction

A rather great amount of numerical solution techniques for the convection-diffusion equations are available today. Nevertheless the fact that these equations lie on the basis of the mathematical modeling of countless physical phenomena, keeps encouraging specialists in the search for efficient methodology to solve this class of problems. This is particularly true of convection dominated processes, which often reveal demerits of widespread computational techniques, even when the problem to solve is linear.

This work is primarily aimed at carrying out a complete mathematical study of the Hermite finite element method proposed in [START_REF] Ruas | A Hermite finite element method for convection-diffusion equations[END_REF], to solve the convection-diffusion equations. More specifically such a method is an extension to the convection-diffusion equations of the Hermite finite element method introduced in [START_REF] Ruas | Hermite finite elements for second order boundary value problems with sharp gradient discontinuities[END_REF] for pure diffusion problems. We recall that the latter method can be regarded as a variant of the well-established Raviart-Thomas mixed method of the lowest order, also known as the RT 0 method, to solve the diffusion equations in two-and three-dimensional space. The method of [START_REF] Ruas | A Hermite finite element method for convection-diffusion equations[END_REF] in turn is a Hermite analog applied to the case of the convection-diffusion equations, of one of the RT 0 method's extensions proposed by Douglas & Roberts in this journal [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF] to solve more general second order boundary-value problems. A Hermite analog of the other version of the RT 0 method due to the latter authors, also studied in [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF], is considered in this work, though from a purely numerical point of view.

Historically Hermite finite elements have mostly been used to solve fourth order partial differential equations, because minimum continuity of solution derivatives across inter-element boundaries is required in this case. However the construction of such elements can be rather laborious, as shown in [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. It is noticeable in this respect that the recent technique of the virtual element led to feasible constructions of C k functions for k ≥ 1 on meshes consisting of polygonal elements of arbitrary shape [START_REF] Beirão Da Veiga | Basic principles of virtual element methods[END_REF], though by means of polynomials of rather high degree.

On the other hand Hermite interpolation has been showing to be a good alternative to solve several kinds of field problems modeled by second order boundary value problems in many respects. An outstanding demonstration of such an assertion is provided by the isogeometric analysis (IGA) introduced in the last decade (see e.g. [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF]). In this case advantage is taken from data satisfying high continuity requirements supplied by CAD, for a subsequent finite element analysis. However IGA in connection with triangular or tetrahedral meshes is incipient, in spite of the undeniable geometric flexibility of this type of partitions. This is a good reason to study Hermite finite elements methods defined upon triangles or tetrahedra to solve second order partial differential equations, which are low order and easy to implement at a time. That is what we do in this work, by focusing more particularly the representation of fluxes for the simulation of phenomena or processes of the convection-diffusion type.

In practice quantities directly depending on partial derivatives of the variable in terms of which an equation is expressed, i.e., the primal variable, are often more important than this unknown itself. Among them one might quote the flux in a porous medium flow or in heat flow. As far as methods allowing to enforce the continuity of normal derivatives or normal fluxes across the boundaries of triangular or tetrahedral cells are concerned, both mixed finite elements and finite volumes have been playing a prominent role since long. In particular the RT 0 method is a popular numerical tool to solve diffusion-like problems providing flux continuity across inter-element boundaries. As recalled above, Douglas & Roberts [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF] extended the method to linear second order boundary-value problems including the convection-diffusion equations. The main drawback of these methods however is the poor representation of the primal variable by piecewise constant functions. That is why many authors attempted to enhance the quality of approximation of the primal variable through post-processing or hybridization, among other techniques.

In contrast to such approaches a Hermite analog of the RT 0 method can be used in the direct solution of the diffusion equations. It proved to be a well adapted alternative to attain higher convergence rates without increasing the computational effort, as shown in [START_REF] Ruas | Hermite finite elements for diffusion phenomena[END_REF]. This rather well-off experience encouraged the authors to further study two finite element methods of the Hermite type based on a quadratic interpolation, to solve the convection-diffusion equations. Both can be viewed as Hermite analogs of the Douglas & Roberts' extensions of the RT 0 method. However the former have either identical or better convergence properties than the latter, according to the norm under consideration. Once again this is achieved at practically the same implementation cost.

An outline of this paper is as follows. In Section 2 we introduce some notations and specify the model problem, which all the studies conducted here apply to. The method proposed in [START_REF] Ruas | A Hermite finite element method for convection-diffusion equations[END_REF] is described in more detail in Section 3, where we recall that it is uniformly stable with respect to a suitable working norm. In Section 4 we apply these results, which immediately lead to first order error estimates in the same norm. We also prove, through duality techniques, that the method's convergence order in the L 2 -norm is two, in contrast to the first order ones that hold for the mixed extension of the RT 0 method to convection-diffusion-reaction equations in non divergence form proposed in [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF]. In Section 5 we consider a variant of the method studied in Sections 3 and 4, which can be viewed as the Hermite analog of the Douglas & Roberts mixed method [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF], applying to the C-D equations in divergence form. Both methods are compared, either to each other or to the corresponding Douglas & Roberts methods, by checking their convergence properties and accuracy in different senses, in the light of numerical experiments reported in Section 6. In Section 7 we conclude with some comments on the whole work.

Notations and model problem

Let Ω be a bounded Lipschitz domain of N , N = 2, 3, with boundary Γ. Referring to [START_REF] Adams | Sobolev Spaces[END_REF], in the sequel we employ the following notations: S being a proper subset of Ω, we denote the standard norm of Sobolev spaces H m (S) (resp. W m,p (S) for 1 ≤ p ≤ ∞, p = 2), for any non negative integer m by • m,S (resp. • m,p,S ) including L 2 (S) ≡ H 0 (S) (resp. L p (S) ≡ W 0,p (S)). For m > 0 the standard semi-norm of H m (Ω) (resp. W m,p (Ω) for 1 ≤ p ≤ ∞, p = 2), i.e. the standard norm of

H m 0 (Ω) (resp. W m,p 0 (Ω) for 1 ≤ p ≤ ∞, p = 2), is denoted by | • | m (resp. | • | m,p
). Further for any non negative integer m and 1 ≤ p ≤ ∞, p = 2, the standard norm of W m,p (Ω) will be denoted by • m,p ; moreover ∀f, g ∈ L 2 (S), (f, g) S := S f g dS and f S := [(f, f ) S ] 1/2 and we set (f, g)

:= Ω f g dx ∀f, g ∈ L 2 (Ω) and f := [(f, f )] 1/2 .
Let f be a given function in L 2 (Ω), K be a tensor assumed to be constant, symmetric and positivedefinite and w ∈ [C 0 ( Ω)] N denote a velocity field. In this work we study as a model the following equation:

Find

u ∈ H 1 0 (Ω) such that -∇ • K∇u + w • ∇u = f in Ω. (1)
Equation ( 1) is assumed to have a unique solution, which is guaranteed in some important particular cases. For instance, the case where the divergence of w is bounded in Ω and C D := ∇ • w 0,∞ is sufficiently small. Indeed by the Divergence Theorem we easily obtain:

-(∇ • K∇u, u) + (w • ∇u, u) = (K∇u, ∇u) -(∇ • wu, u)/2 ≥ λ min ∇u 2 -C D u 2 /2,
where λ min is the smallest eigenvalue of K. It follows that whenever

λ min > C D C 2 P /2, where C P is the constant of the Friedrichs-Poincaré inequality v ≤ C P ∇v ∀v ∈ H 1 0 (Ω), (1)
has a unique solution by the Lax-Milgram Theorem. Notice that in many applications w is solenoidal, and in this case existence and uniqueness of a solution to (1) always hold true. Another situation where this problem has a unique solution arises when a convection-diffusion process is taking place at a low Péclet number Pé (see e.g. [START_REF] Ikeda | Maximum Principle in Finite Element Models for Convection-diffusion Phenomena[END_REF]). Just to make ideas clear we can take Pé = w 0,∞ C P /(2λ min ). Then in case Pé < 1/2, again by the Lax-Milgram Theorem, (1) will have a unique solution.

A Hermite solution method

Henceforth we assume that Ω is a polygon if N = 2 or a polyhedron if N = 3, and that we are given a finite element partition T h of Ω, consisting of triangles or tetrahedra according to the value of N , and belonging to a regular family of partitions (cf. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). h denotes the maximum diameter of the elements of T h .

In the following paragraphs we define two finite element spaces U h and V h associated with T h . Let w h be the constant field in each element of T ∈ T h whose value in T is w(x T ), where x T is the position vector of the centroid of T , and w 1 h be the standard continuous piecewise linear interpolate of w at the vertices of T h . We further introduce the operators Π T : L 2 (T ) -→ L 2 (T ) given by Π T [v] := T vdx/meas(T ) for T ∈ T h , and

Π h : L 2 (Ω) -→ L 2 (Ω) by Π h [v] |T = Π T [v |T ] ∀T ∈ T h .
Now throughout this work we will work with the following, Local algebraic structure of the Hermite finite element spaces:

Every function v ∈ V h (resp. ∈ U h ) is such that in each element T ∈ T h it is expressed by v |T = x t {K -1 [ax/2 + b]} + d, (2)
where x represents the space variable, b is a constant vector of N and a and d are two real coefficients. In every N -simplex T we associate with a quadratic function v of the form (2):

Sets D T and E T of local degrees of freedom for the Hermite finite element spaces: F being an edge if N = 2 or a face if N = 3 belonging to the boundary ∂T of an N -simplex T , and n F being the unit normal vector on F oriented in a given manner for each F ⊂ ∂T , we set:

   D T := {∪ F U F } ∪ U T where U F (v) = F K∇v • n F ds/meas(F ); U T (v) = T v dx/meas(T ). (3)
   E T := {∪ F V F } ∪ V T where V F (v) = F (K∇v + w h Π T [v]) • n F ds/meas(F ); V T (v) = T v dx/meas(T ). (4)
The canonical basis functions associated with these sets of degrees of freedom are as follows. First we note that ∀v ∈ This should be enough to determine the N + 2 basis functions associated with a given N -simplex T , corresponding to the sets U T and V T of degrees of freedom, for spaces U h and V h respectively, since the RT 0 method is well-known (cf. [START_REF] Thomas | Sur l'analyse numérique des méthodes d'éléments finis hybrides et mixtes[END_REF]). However for the sake of clarity we exhibit them below.

V h or ∈ U h , ∇v |T for T ∈ T h is expressed by K -1 [a T x + b T ] for certain a T ∈
T being an element of T h let x T i be the position vector of the i-th vertex S T i of T , F T i be the face of T opposite to S T i and h T i be the length of the corresponding height of T , for i = 1, . . . , N + 1. We have:

Local basis functions ϕ T i for space U h : The local basis function ϕ T i associated with the degree of freedom U F T i and the basis function ϕ T N +2

associated with the degree of freedom U T are given by:

                     ϕ T i = x t {K -1 [a T i x/2 + b T i ]} + d T i for i = 1, . . . , N + 2, where a T i = [h T i ] -1 , b T i = -x T i a T i , d T i = -T x t {K -1 [a T i x/2 + b T i ]} dx/meas(T ),    for i = 1, . . . , N + 1; a T N +2 = 0, b T N +2 = (0, . . . , 0) t , d T N +2 = 1. (5)
Local basis functions ψ T i for space V h : Akin to the case of U h , the ψ T i 's are functions of the form (2) for i = 1, . . . , N + 2. Since by definition the mean values of all the first N + 1 ϕ T i 's vanish, the local basis functions ψ T i of V h associated with the degree of freedom V F T i for i = 1, . . . , N + 1, together with its local basis function ψ T N +2 associated with the degree of freedom V T are given by:

ψ T i = ϕ T i for i = 1, . . . , N + 1 ψ T N +2 = [x T -x] t K -1 w h + 1. (6)
Next we define, Hermite finite element spaces U h and V h : As already stated in our previous definitions for element's local algebraic structure, for every interface F (an inner edge for N = 2 and an inner face for N = 3) of two elements in T h , n F is oriented in a given manner for both of them. Then every function in v ∈ V h (resp. ∈ U h ) is such that its restriction to every T ∈ T h is a N + 2 coefficient quadratic function of the form (2), whose degrees of freedom of the type V F (resp. U F ) coincide on both sides of every interface F of a pair of elements in T h .

We proceed by setting the discrete variational problem [START_REF] Burman | Bubble stabilized discontinuous Galerkin method for parabolic and elliptic problems[END_REF] below, aimed at approximating (1), whose bi-linear form a h and linear form L h are given by ( 8):

Find

u h ∈ U h such that for all v ∈ V h a h (u h , v) = L h (v), (7)
holds, where ∀u ∈ U h and ∀v ∈ V h ,

       a h (u, v) := T ∈T h [(∇ • K∇u -w 1 h • ∇u, Π T [v]) T +(∇u, K∇v + w h Π T [v]) T + (u, ∇ • K∇v) T ]; L h (v) := -(f, Π h [v]). (8)
It is noteworthy that whenever w 1 h = w h the two terms involving w in the expression of a h cancel out. Hence apparently problem's dependence on the convective velocity is not taken into account by formulation [START_REF] Burman | Bubble stabilized discontinuous Galerkin method for parabolic and elliptic problems[END_REF] if w happens to be constant. However one should bear in mind that even in this case such a dependence remains implicit therein, through the construction of space V h . Now let us consider the space

V := {v|v ∈ H 1 (Ω); ∇ • K∇v ∈ L 2 (Ω)}.
Clearly a h can be extended to (U h + V ) × (V h + V ). Then we further introduce the functional • h : U h + V h + V -→ given by:

v h :=   (Π h v, Π h v) + T ∈T h {(∇v, ∇v) T + (∇ • K∇v, ∇ • K∇v) T }   1/2 . ( 9
)
The expression • h obviously defines a norm over V , U h and V h . In this manner, it is not difficult to establish the continuity of a h over (U h + V ) × (V h + V ) with a mesh independent constant M (cf. the proof of Proposition 3.1 hereafter):

a h (u, v) ≤ M u h v h . (10)
On the other hand there is no way for a h to be coercive. Hence we resort to an inf-sup condition for a h over U h × V h [START_REF] Babuška | The finite element method with Lagrange multipliers[END_REF], which directly implies that (7) has a unique solution. More specifically the following stability result was proved in [START_REF] Ruas | A Hermite finite element method for convection-diffusion equations[END_REF].

Proposition 3.1 ([23]) If h is sufficiently small and w ∈ [W 1,∞ (Ω)] N , there exists a constant α > 0 independent of h such that ∀u ∈ U h sup v∈V h \{0} a h (u, v) v h ≥ α u h . (11)
In the next section we derive estimates for u -u h h using a modified Strang Lemma for non coercive problems given in [START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF]. In this aim we have to consider the following auxiliary problem:

Find u * h ∈ U h such that for all v ∈ V h a * h (u * h , v) = L h (v), (12)
holds, where ∀u ∈ U h + V and ∀v ∈

V h + V , a * h (u, v) := T ∈T h [(∇ • K∇u -w • ∇u, Π T [v]) T +(∇u, K∇v + w h Π T [v]) T + (u, ∇ • K∇v) T ]. (13)
Similarly to the case of problem (7) (cf. [START_REF] Ruas | A Hermite finite element method for convection-diffusion equations[END_REF]) we can prove, Theorem 3.2 Problem ( 12) has a unique solution and moreover there exists a constant

C * independent of h such that u * h h ≤ C * f . (14)
Before proving Theorem 3.2 we establish a stability result for problem [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF].

Proposition 3.3 If h is sufficiently small and w ∈ [W 1,∞ (Ω)] N , there exists a constant α * > 0 inde- pendent of h such that ∀u ∈ U h sup v∈V h \{0} a * h (u, v) v h ≥ α * u h . ( 15
) Proof: Given u ∈ U h define v := v 1 + v 2 + v 3 , where v i ∈ V h for i = 1, 2, 3 are defined as follows: v 1 = θ 1 w 1 , θ 1
being a non negative constant to be specified, and w 1 being defined by

Π h [w 1 ] = Π h [u], together with (K∇w 1 + w h Π h [w 1 ]) • n T = (K∇u) • n T for every T ∈ T h , where n T is the outer normal on ∂T . v 2 equals θ 2 ∇ • K∇u in every T ∈ T h
, where θ 2 is a non negative constant to be specified. v 3 is constructed by applying Theorem 4 of [START_REF] Raviart | Mixed Finite Element Methods for Second Order Elliptic Problems[END_REF]. According to it there exists a field p ∈ Q h := {q | ∃u ∈ U h such that q |T = K∇u |T ∀T ∈ T h }, satisfying for a constant C independent of h:

∇ • p = Π h [u] in Ω; p ≤ C Π h [u] . (16)
Then recalling that the normal traces over the faces of the elements in T h of fields belonging to Q h are constant [START_REF] Raviart | Mixed Finite Element Methods for Second Order Elliptic Problems[END_REF] v 3 is defined in such a way that ∀T ∈ T h , (K∇v

3 + w h Π T [v 3 ]) • n T = p • n T ∀T ∈ T h and Π h [v 3 ] = -θ 1 Π h [u]. It is clear that ∇ • K∇w 1 = ∇ • K∇u.
Moreover by construction we have,

∂T K∇w 1 • n T w 1 dS -(∇ • Kw 1 , w 1 ) T = (K∇w 1 , ∇w 1 ) T = (K∇u -w h Π T [u], ∇w 1 ) T ∀T ∈ T h .
(17) Then λ and Λ being the smallest and the largest eigenvalue of K, after straightforward manipulations it follows that

∇w 1 T ≤ λ -1 (Λ ∇u T + w 0,∞ Π T [u] T ). (18)
This implies that for

C1 = λ -1 (Λ 2 + w 2 0,∞ ) 1/2 , T ∈T h ∇w 1 2 T ≤ C2 1 u 2 h , which immediately yields, v 1 h ≤ C 1 u h , with C 1 = θ 1 (1 + C2 1 ) 1/2 . ( 19
)
As for v 2 we trivially have,

v 2 h ≤ C 2 u h , with C 2 = θ 2 . ( 20
)
On the other hand by construction

v 3 fulfills ∇ • [K∇v 3 ] |T = ∇ • p |T = Π T [u]
, ∀T ∈ T h , and hence,

λ ∇v 3 2 T ≤ (K∇v 3 , ∇v 3 ) T = ∂T (K∇v 3 )•n T v 3 dS -(v 3 , ∇•p) T = (θ 1 w h Π T [u]+p, ∇v 3 ) T (21)
It easily follows that

v 3 h ≤ C 3 Π h [u] ≤ C 3 u h , where C 3 = [( C + θ 1 w 0,∞ ) 2 λ -2 + θ 2 1 + 1] 1/2 . ( 22
)
Now taking into account ( 16) and ( 17), after straightforward calculations we obtain,

a * h (u, v 1 ) = θ 1 T ∈T h {2(∇ • K∇u, Π T [u]) T + (K∇u, ∇u) T -(w • ∇u, Π T [u]) T }; (23)
a * h (u, v 2 ) = θ 2 T ∈T h { ∇ • K∇u 2 T +([w h -w] • ∇u, ∇ • K∇u) T }; (24)
a * h (u, v 3 ) = T ∈T h { Π T [u] 2 T +θ 1 [(w • ∇u, Π T [u]) T -(∇ • K∇u, Π T [u]) T ] + (p, ∇u) T }. ( 25
)
Then, recalling the definition of w h , there exists a mesh independent constant C W (cf. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]) such that ww h 0,∞,T ≤ C W h ∇w 0,∞,T ∀T ∈ T h . Using this fact, together with ( 23)-( 24)-( 25), simple manipulations lead to:

                     a * h (u, v) ≥ Π h [u] 2 + T ∈T h θ 1 λ ∇u 2 T + θ 2 2 ∇ • K∇u 2 T -p T ∇u T - T ∈T h θ 1 ∇ • K∇u T Π T [u] T ) + θ 2 C 2 W h 2 |w| 2 1,∞ ∇u 2 T /2 ≥ Π h [u] 2 /4 + T ∈T h {(θ 1 λ -C2 -θ 2 C 2 W h 2 |w| 2 1,∞ /2) ∇u 2 T +(θ 2 /2 -θ 2 1 /2) ∇ • K∇u 2 T } (26) Now if we assume that h 2 ≤ β(C W |w| 1,∞ ) -2 with β ≤ 4λ 2 /[D + (D 2 + 8λ 2 ) 1/2 ] for D = 1 + 4 C2 , we may choose θ 1 > 0 satisfying θ 1 λ -C2 -βθ 2 /2 ≥ 1/4 with θ 2 = 1/2 + θ 2 1 .
It follows from ( 26), [START_REF] Raviart | Mixed Finite Element Methods for Second Order Elliptic Problems[END_REF], [START_REF] Ruas | Automatic generation of triangular finite element meshes[END_REF] and [START_REF] Ruas | Hermite finite elements for diffusion phenomena[END_REF], that,

a * h (u, v) ≥ u 2 h /4; v h ≤ C u h , with C = [3(C 2 1 + C 2 2 + C 2 3 )] 1/2 . ( 27
)
This immediately yields [START_REF] Ikeda | Maximum Principle in Finite Element Models for Convection-diffusion Phenomena[END_REF] with α * = 1/(4C).

Proof of Theorem 3.2: Since V h is a finite dimensional space, according to [START_REF] Babuška | The finite element method with Lagrange multipliers[END_REF] the existence and uniqueness of a solution to [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF] follows from [START_REF] Ikeda | Maximum Principle in Finite Element Models for Convection-diffusion Phenomena[END_REF]. Moreover, combining [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF] and [START_REF] Ikeda | Maximum Principle in Finite Element Models for Convection-diffusion Phenomena[END_REF] we easily obtain,

α * u * h h ≤ sup v∈V h \{0} L h (v) Π h [v] . ( 28
) Since L h (v) = Ω f Π h [v]
dx from (28) we finally derive [START_REF] Grisvard | Elliptic Problems in Non Smooth Domains[END_REF] with

C * = [α * ] -1 .

Convergence results

Henceforth we denote by ∇ h the operator from

V + U h + V h onto L 2 (Ω) defined by [∇ h w] |T = ∇[w |T ] ∀T ∈ T h , ∀w ∈ V + U h + V h .
Notice that for any function u ∈ V + U h , ∇ • K∇u is well-defined in L 2 (Ω) (cf. [START_REF] Thomas | Sur l'analyse numérique des méthodes d'éléments finis hybrides et mixtes[END_REF]) and hence there is no need to use the operator ∇ h in this case.

In order to study the convergence of u h to u in appropriate norms we first note that from the properties of V h and equation ( 1) we easily infer that u satisfies

a * h (u, v) = L h (v) ∀v ∈ V h . (29)
From the continuity of a * h and the uniform stability result proved in Proposition 3.3, we may apply the generalized First and Second Strang's inequality for the weakly coercive case, namely, inequality (32) of [START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF]. In the case under study this writes,

u -u h h ≤ M * α * inf w∈U h u -w h + 1 α sup v∈V h \{0} [a * h -a h ](u * h , v) v h (30)
where M * is a constant such that

a * h (u, v) ≤ M * u h v h ∀u ∈ V + U h , ∀v ∈ V + V h . ( 31
)
Proposition 4.1 There exists a constant M * independent of h such that (31) holds.

Proof: Since (u, ∇ • K∇v) T = (Π T [u], ∇ • K∇v) T ∀T ∈ T h and ∀v ∈ V h , we trivially have,    a * h (u, v) ≤ ( ∇ • K∇ h u + w 0,∞ ∇ h u ) Π h [v] + ∇ h u (Λ ∇ h v + w 0,∞ Π h [v] ) + T ∈T h Π T [u] T ∇ • K∇v T . (32)
(32) immediately yields (31) with

M * = 2 max[1, 2 w 0,∞ , Λ].
Next we prove the validity of the following a priori error estimate for the method under study:

Theorem 4.2 Assume that w ∈ [W 1,∞ (Ω)] N and h is sufficiently small. Then if u ∈ H 2 (Ω) and f ∈ H 1 (Ω) there exists a mesh independent constant C such that, u -u h h ≤ C h[ u 2 + f 1 ]. (33)
Proof: By standard results applying to the RT 0 method, and since Π h [u -u h ] is obviously bounded above by a mesh independent constant times h|u| 1 , for a suitable constant C I independent of h it holds,

inf w∈V h u -w h ≤ C I h[ u 2 +|f | 1 ]. (34)
On the other hand we have

|[a * h -a h ](u * h , v)| = |([w -w 1 h ] • ∇ h u * h , Π h [v])|. Hence for a mesh inde- pendent constant C * W such that w -w 1 h 0,∞ ≤ C * W h|w| 1,∞ we derive, |[a * h -a h ](u * h , v)| ≤ C * W h|w| 1,∞ ∇ h u * h v h . (35)
Taking into account ( 14), ( 30)-( 34)-(35) readily yield (33), C being a mesh independent constant.

Next we give a fundamental result of this work:

Theorem 4.3 If Ω is convex, w ∈ [W 2,4 (Ω)] N , u ∈ H 2 (Ω), f ∈ H 1 (Ω)
and h is sufficiently small, there exists a constant C independent of h such that,

u -u h ≤ C h 2 ( u 2 + f 1 ). ( 36
)
Proof: The proof is a non-trivial extension of the proof of Theorem 2.3 in [START_REF] Ruas | Hermite finite elements for second order boundary value problems with sharp gradient discontinuities[END_REF], where the quadratic convergence was shown for the pure diffusion case. We first observe that w ∈ [W 1,∞ (Ω)] N according to the Sobolev Embedding Theorem [START_REF] Adams | Sobolev Spaces[END_REF]. Moreover using the definitions of a h and Π h , together with the continuity of the normal components of

K∇v h + w h Π h [v h ] on ∂T for v h ∈ V h , we easily obtain, a h (u -u h , v h ) + ((w 1 h -w) • ∇u, Π h [v h ]) = 0, (37)
for all v h ∈ V h . Similarly, owing to the continuity of the normal components of K∇u h on ∂T (cf. [START_REF] Ruas | Hermite finite elements for second order boundary value problems with sharp gradient discontinuities[END_REF]):

(u-u h , ∇•K∇v) = a h (u-u h , v)+(∇•K∇(u-u h ), v-Π h [v])+((w 1 h -w h )•∇ h (u-u h ), Π h [v]), (38)
for all v ∈ {v|v ∈ H 1 0 (Ω), ∇ • K∇v ∈ L 2 (Ω)}. By using the Aubin-Nitsche trick and (38) we can write [START_REF] Evans | Partial Differential Equations[END_REF]). Thus it is easy to see that due to (39) it holds,

u -u h = sup v∈D(Ω)\{0} (u -u h , ∇ • K∇v) ∇ • K∇v = sup v∈D(Ω)\{0} a h (u -u h , v) + (∇ • K∇(u -u h ), v -Π h [v]) + ((w 1 h -w h ) • ∇ h (u -u h ), Π h [v]) ∇ • K∇v (39) where D(Ω) := V ∩ H 1 0 (Ω). Let now u = u h (if u = u h then (36) trivially holds) and B D (0, 1) := {v|v ∈ D(Ω), ∇ • K∇v = 1}. We know that there exists v 0 ∈ B D (0, 1) such that ∇ • K∇v 0 = u -u h u -u h (cf.
u-u h = a h (u-u h , v 0 )+(∇•K∇(u-u h ), v 0 -Π h [v 0 ])+((w 1 h -w h )•∇ h (u-u h ), Π h [v 0 ]). ( 40
)
By combining (37) and (40) we further get for any [START_REF] Evans | Partial Differential Equations[END_REF]), we can define the standard interpolate

v h ∈ V h u -u h = a h (u -u h , v 0 -v h ) + (∇ • K∇(u -u h ), v 0 -Π h [v 0 ]) (41) +((w 1 h -w h ) • ∇ h (u -u h ), Π h [v 0 ]) -((w 1 h -w) • ∇u, Π h [v h ]). Since D(Ω) ⊂ H 2 (Ω) in case Ω is convex (cf.
I h v ∈ V h of every v ∈ D(Ω), based on the degrees of freedom of V h . Taking v h = I h v 0 in (41), we get, u -u h = a h (u -u h , v 0 -I h v 0 ) + (∇ • K∇(u -u h ), v 0 -Π h [v 0 ]) (42) +((w 1 h -w h ) • ∇ h (u -u h ), Π h [v 0 ]) -((w 1 h -w) • ∇u, Π h [I h v 0 ]).
By using the continuity (10) of a h and the definition of v 0 , together with the approximation properties of I h for h sufficiently small (see [START_REF] Ruas | Hermite finite elements for second order boundary value problems with sharp gradient discontinuities[END_REF], p. 239 for details), we have for a suitable h-independent constant C,

a h (u -u h , v 0 -I h v 0 ) ≤ 3 4 u -u h +C u -u h h [ v 0 -I h v 0 + ∇ h (v 0 -I h v 0 ) ]. (43)
Applying to (42) and (43) standard results to estimate

v 0 -Π h [v 0 ] together with v 0 -I h v 0 + ∇ h (v 0 -I h v 0 )
, recalling (33) we easily conclude that there exists another constant C independent of h such that,

u -u h ≤ Ch 2 (|u| 2 + |f | 1 ) + 4[|((w 1 h -w h ) • ∇ h (u -u h ), Π h [v 0 ])| + |((w 1 h -w) • ∇u, Π h [I h v 0 ])|].
(44) We proceed by estimating the two terms in brackets on the right hand side of (44) denoted by T 1 and T 2 . First we note that from the Sobolev Embedding Theorem and the convexity of Ω (cf. [START_REF] Grisvard | Elliptic Problems in Non Smooth Domains[END_REF]), there exist constants C ∞ and C ∞ depending only on Ω such that

Π h [v 0 ] 0,∞ ≤ v 0 0,∞ ≤ C ∞ v 0 2 ≤ C ∞ .
Thus for a mesh independent constant C 2 we have:

T 1 ≤ w 1 h -w h ∇ h (u -u h ) Π h [v 0 ] 0,∞ ≤ C 2 h 2 |u| 2 . ( 45
)
For deriving the estimate (45) we used the fact that both w 1 h and w h are interpolates of w, the result (33) and the boundedness of Π h [v 0 ] 0,∞ . C 2 is the product of C ∞ with another constant not depending on h and the semi-norm of w in [H 1 (Ω)] N . The term T 2 can be estimated in a similar way. Using now the fact that w 1 h is a piecewise linear interpolate of w, the regularity of the solution u and standard properties of Π h and I h , there holds

T 2 ≤ w 1 h -w 0,4 ∇u 0,4 Π h [I h v 0 ] ≤ C 3 h 2 u 2 , (46)
where C 3 equals an h-independent constant times |w| 2,4 . Putting together (44)-( 45)-( 46), we obtain the quadratic convergence (36).

A variant for the equations in divergence form

Like in [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF] it is possible to consider a variant of the method described in Section 3 applying to the case where the normal component of the total flux -K∇u + wu is continuous across the element interfaces. In the case of the mixed formulation this corresponds to introducing the auxiliary variable p given by the above expression, and write the C-D equation equation ( 1) in divergence form, namely Find u satisfying u = 0 on Γ and p such that

∇ • p -∇ • w u = f in Ω, p + K∇u -wu = 0 in Ω. ( 47
)
Recalling the space H(div;

Ω) := {q | q ∈ [L 2 (Ω)] N , ∇ • q ∈ L 2 (Ω)
}, a natural weak (variational) formulation equivalent to system (47) is given in [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF], that is, Find u ∈ L 2 (Ω) and p ∈ H(div; Ω) such that for all v ∈ L 2 (Ω), and for all q ∈ H(div; Ω),

(∇ • p, v) -(∇ • w u, v) = (f, v) (K -1 p, q) -(u, ∇ • q) -(K -1 w u, q) = 0. (48)
The extension of RT 0 to the C-D equation considered in [START_REF] Douglas | Mixed finite element methods for second order elliptic problems[END_REF] consists of using the Raviart-Thomas interpolation of the lowest order to represent p and q -i.e. to approximate H(div; Ω) -, and the space of constant functions in each element of the partition T h to represent u and v. In contrast, here we shall mimic (48) by resorting to the space U h , after adding up both relations in (48). More specifically we take in each element T ∈ T h , q |T = K∇v |T for v ∈ U h . Now u h will be searched for in a space W h defined hereafter. First we have to construct field a wh to replace w h (cf. Section 3), in order to preserve optimality of the approximation of u. In this aim it suffices that wh be of the form cx + d in each T ∈ T h for suitable real number c and real vector d = [d 1 , . . . , d N ] t . This representation is compatible with the requirement that the normal component of the flux variable p = -K∇ + wu be continuous across the mesh edges at discrete level. The natural choice of wh is certainly the interpolate of w in the Raviart-Thomas (RT 0 ) space. Now we define the Hermite finite element space W h : W h is the space of functions v of the form x t {K -1 [ax/2 + b]} + d in every T ∈ T h , such that the mean normal flux F (-K∇v + wh Π h [v]) • n F ds/meas(F ) is continuous across all the inner edges or faces F of the partition. This is about all that is needed to complete the definition of W h . Indeed using a procedure very similar to the one in Section 3 (cf. ( 5)) it is possible to uniquely determine the N +2 local basis functions η T i for each N -simplex T ∈ T h related to the above set of degrees of freedom completed with the function mean value in T , defining W h locally. More precisely, setting [ wh ] |T = a T w x + b T w , since by definition Π T η T i = 0 for i = 1, . . . , N + 1 and Π T η T N +2 = 1, recalling (5) we have:

η T i = ϕ T i for i = 1, . . . , N + 1, η T N +2 = x t {K -1 [a T w x/2 + b T w ]} + 1 -T x t {K -1 [a T w x/2 + b T w ]}dx/meas(T ). (49)
Now we replace in (48) :

• u with Π h [u h ];
• w with wh ;

• p with -K∇ h u h + wh Π h [u h ] (taking u h ∈ W h ); • q with -K∇ h v (taking v ∈ U h ); • f with Π h [f ].
This leads to the following equation:

T ∈T h [(∇ • {K∇u h -wh Π T [u h ]}, v) T + (∇ • wh Π T [u h ], v) T + (K∇u h -wh Π T [u h ], ∇v) T + (Π T [u h ], ∇ • K∇v) T + ( wh Π T [u h ], ∇v) T ] = -(Π h [f ], v) ∀v ∈ U h . (50)
After straightforward simplifications, and taking into account that

(Π h [f ], v) = (f, Π h [v]
), we come up with the following Hermite finite element counterpart of (1):

Find

u h ∈ W h such that for all v ∈ U h ãh (u h , v) = L h (v), (51)
where

∀u ∈ V + W h and ∀v ∈ V + U h ,    ãh (u, v) := T ∈T h [(∇ • K∇u, v) T + (∇u, K∇v) T + (u, ∇ • K∇v) T ] L h (v) := -(f, Π h [v]). (52)
At a first glance (52) seems to indicate that the velocity w does not appear in formulation (51). Nonetheless w remains implicit therein through the definition of space W h .

The fact that problem (51) has a unique solution can be established quite similarly to problem [START_REF] Burman | Bubble stabilized discontinuous Galerkin method for parabolic and elliptic problems[END_REF]. The convergence results that hold for this method can be proved very much like in the case of the method defined in Section 3. The main difference is that it is necessary to require a little more regularity of ∇•w, namely, that this function lies in W 1,∞ (Ω). Apart from this assumption, the results are qualitatively equivalent, in the sense that a priori error estimates completely analogous to those of Theorem 4.2 and Theorem 4.3 apply to problem (51) as well. As far as this work is concerned, resulting properties among others we have not formally established here, are illustrated by means of numerical examples given in the following section.

Numerical experiments

In this section we present some numerical results obtained with the methods described in Sections 3 and 5 for two test problems, which particularly highlight their behavior. The following nomenclature is used for the different numerical methods being experimented:

• Method A -Douglas & Roberts version in non divergence form of mixed method RT 0 ;

• Method HA -Hermite analog of Method A (cf. Section 3);

• Method B -Douglas & Roberts version in divergence form of mixed method RT 0 ;

• Method HB -Hermite analog of Method B (cf. Section 5).

In this section we will denote by e h the error function u -u h , and by ∆ h g the function defined by [∆ h g] |T = ∆[g |T ] ∀T ∈ T h , for any function g whose laplacian is well defined in the interior of every mesh element. Moreover the expression pseudo maximum norm will stand for the maximum absolute value of a function at the centroids of the mesh elements.

Test-problem 1: In these experiments Ω is the unit square and a manufactured solution u is given by u

(x 1 , x 2 ) = (x 1 -x 2 1 )(x 2 -x 2
2 )/4. This together with the choice K = I and w =Pé[x 2 1 , x 2 2 ] t / √ 2 where Pé is the Péclet number, produces a right hand side datum f . A sequence of uniform meshes was employed with 2L 2 triangles, for L = 8, 16, 32, 64, constructed by first subdividing Ω into L 2 equal squares and then each one of these squares into two triangles by means of their diagonals parallel to the line x 1 = x 2 . Quite abusively we denote by h the spacing 1/L.

In Tables 1,2, 3 and 4 we display the absolute errors in four different respects for increasing values of L, of the approximate solutions obtained with methods A, HA, B and HB for Pé= 1. More precisely the absolute errors of u, ∇u and ∆u = ∇ • K∇u measured in the norm of L 2 (Ω), and in the pseudo maximum norm, are shown in Tables 1,2, 3, 4, respectively. In Tables 5,6 1/16 0.54815217 x 10 -3 0.54824708 x 10 -3 0.44953727 x 10 -4 0.45890562 x 10 -4

1/32 0.27439727 x 10 -3 0.27440932 x 10 -3 0.11286733 x 10 -4 0.11531497 x 10 -4

1/64 0.13723841 x 10 -3 0.13723992 x 10 -3 0.28250216 x 10 -5 0.28866263 x 10 -5 1/32 0.11636060 x 10 -2 0.11951094 x 10 -2 0.11636981 x 10 -2 0.11919738 x 10 -2

1/64 0.58218263 x 10 -3 0.59798255 x 10 -3 0.58219418 x 10 -3 0.59640799 x 10 -3 • Methods A and HA are fairly equivalent to Methods B and HB in all respects for a low Péclet number.

• Methods A and HA are superior to Methods B and HB in all respects when the Péclet number is not low. 1/16 0.29448711 x 10 -4 0.47404901 x 10 -4 0.44100102 x 10 -4 0.66030394 x 10 -4

1/32 0.78792148 x 10 -5 0.12849987 x 10 -4 0.11178619 x 10 -4 0.17190429 x 10 -4

1/64 0.20428256 x 10 -5 0.33363926 x 10 -5 0.28130033 x 10 -5 0.43809703 x 10 -5 1/32 0.27441391 x 10 -3 0.41669269 x 10 -3 0.10350476 x 10 -4 0.37772184 x 10 -3

1/64 0.13724039 x 10 -3 0.15861798 x 10 -3 0.25386722 x 10 -5 0.95399515 x 10 -4 1/16 0.25481792 x 10 -2 0.59439183 x 10 -1 0.28082657 x 10 -2 0.56152848 x 10 -1

1/32 0.11934791 x 10 -2 0.27928892 x 10 -1 0.12388051 x 10 -2 0.26448385 x 10 -1

1/64 0.58595099 x 10 -3 0.13745079 x 10 -1 0.59256341 x 10 -3 0.13029605 x 10 -1 1/16 0.10035771 x 10 +0 0.12162465 x 10 +0 0.59883623 x 10 -1 0.12162465 x 10 +0

1/32 0.51102598 x 10 -1 0.60964240 x 10 -1 0.42990111 x 10 -1 0.60964240 x 10 -1

1/64 0.15249263 x 10 -2 0.15339429 x 10 -2 0.15249297 x 10 -2 0.15339429 x 10 -2 1/16 0.88757842 x 10 -4 0.20869472 x 10 -2 0.11114041 x 10 -3 0.24367046 x 10 -2

1/32 0.14427021 x 10 -4 0.54988868 x 10 -3 0.21157560 x 10 -4 0.65466283 x 10 -3

1/64 0.26841141 x 10 -5 0.13998920 x 10 -3 0.37752993 x 10 -5 0.16575725 x 10 -3

Table 8: Test-problem 1 with Pé=100 -Maximum errors of u at centroids for methods A, B, HA, HB

• The theoretical results of Section 4 for Method HA were confirmed in the case of both a low and a moderate Péclet number.

• The numerical convergence rate in the pseudo maximum norm is approximately two for all the four methods for Pé=1;

• The numerical convergence rate for Pé=100 in the pseudo maximum norm is significantly greater than two for methods A and HA, but remains close to two for methods B and HB.

• For Pé= 100 the mixed methods are a little more accurate than their Hermite counterparts, as far as errors at the triangle centroids are concerned.

• Both ∇u and ∆u tend to be equally approximated by a mixed method and its Hermite counterpart.

The numerical results for Methods B and HB deteriorated substantially as we switched to higher Péclet numbers, which was partially the case of Method HA, while most of the results obtained with Method A remained quite reasonable. Taking L = 64 we illustrate the behaviour of Methods A and HA in Tables 9 and10, respectively, for increasing Péclet numbers. More precisely we took Pé= 10 2k for k = 0, 1, 2, 3, for which we display the absolute errors of u, ∇u, ∆u measured in the norm of L 2 (Ω) and the absolute error of u measured in the pseudo maximum norm. Test-problem 2: In order to observe the behaviour of the four methods being checked, in the presence of a curved boundary, in this test-problem the domain is a disk with unit radius. The manufactured solution u is given by u(x 1 , x 2 ) = (1 -x 2 1 -x 2 2 )/4. Taking again K = I, the right hand side function

Pé e

h ∇ h e h ∆ h e h max T |e h (x T )| 1
f = 1 -(x 2 1 + x 2
2 )/2 corresponds to a convective velocity w =Pé[x 1 , x 2 ] t . For symmetry reasons the computational domain Ω is only the quarter of disk given by x 1 > 0 and x 2 > 0. A sequence of quasi-uniform meshes with 2L 2 triangles was employed for L = 8, 16, 32, 64, constructed by mapping the meshes of Test-problem 1 into the actual meshes of Ω using the transformation of cartesian into polar coordinates in the way described in [START_REF] Ruas | Automatic generation of triangular finite element meshes[END_REF]. For this procedure we have h = 1/L. We denote by Ω h the approximation of Ω consisting of the union of the triangles in T h . In Tables 11,[START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF], 13, 14 we display the absolute errors in four different respects for increasing values of L, of the approximate solutions obtained with methods A, HA, B and HB for Pé= 1. The errors and corresponding notations are the same as in Tables 1 through 8. More precisely the absolute errors of u, ∇u and ∆u = ∇ • K∇u measured in the norm of L 2 (Ω h ) together with the pseudo maximum norm are shown in Tables 11,[START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF][START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF][START_REF] Grisvard | Elliptic Problems in Non Smooth Domains[END_REF]respectively. From these tables we infer that:

• Methods A and HA are superior to Methods B and HB in all respects, except in the approximation of (constant) ∆u, which is almost exactly approximated by all the four methods. 1/16 0.26406317 x 10 -7 0.26406317 x 10 -7 0.26406293 x 10 -7 0.26406317 x 10 -7

1/32 0.26410294 x 10 -7 0.26410294 x 10 -7 0.26410290 x 10 -7 0.26410294 x 10 -7

1/64 0.26411289 x 10 -7 0.26411289 x 10 -7 0.26411287 x 10 -7 0.26411289 x 10 -7 1/16 0.30221612 x 10 -4 0.24492318 x 10 -1 0.10040660 x 10 -3 0.25368398 x 10 -1

1/32 0.75493397 x 10 -5 0.25794826 x 10 -1 0.25126450 x 10 -4 0.26020645 x 10 -1

1/64 0.18777283 x 10 -5 0.26328682 x 10 -1 0.63028673 x 10 -5 0.26385490 x 10 -1

Table 14: Test-problem 2 with Pé=1 -Maximum errors of u at centroids for methods A, B, HA, HB

• Methods A and HA do not seem to be affected by the curved boundary approximation by polygons, while this seems to be case of Methods B and HB.

• Method A and HA approximate both ∇u and ∆u to machine precision; this is an expected behavior since both functions in this test-problem can be exactly represented by the same underlying incomplete linear and constant interpolation for both methods.

• The approximations of u by Method HA converge as an O(h 2 ) in L 2 (Ω h ), while those computed by Method A converge as an O(h), the best we can hope for.

• The numerical convergence rate in the pseudo maximum norm is approximately two for both Method A and Method HA with an advantage of the former over the latter in terms of accuracy.

Akin to Test-problem 1, we checked the behaviour of Methods A and HA as the Péclet number increases. Here again we took L = 64 and Pé= 10 2k for k = 0, 1, 2, 3. The resulting errors measured in the same manner as in Tables 9 and 10 are displayed in From Tables 15 and16 we observe that both Method A and Method HA are accurate to machine precision, irrespective of the Péclet number, as far as the approximations of ∇u and ∆u are concerned. The approximations of u in L 2 (Ω h ) and in the pseudo maximum norm do not seem to be affected by the Péclet number either in this test-problem for both methods. In the former sense Method HA is much more accurate than Method A as expected, while in the latter sense Method A is slightly more precise than Method HA. Notice that this test-problem is a little peculiar, since the exact solution is a quadratic function, whose gradient can be exactly represented by the gradient of the underlying interpolating functions. Actually this also happens to the approximation of the function itself by Method HA, but in this case other sources of errors came into play, such as numerical integration (see also Remark 3 hereafter).

Test-problem 3: Finally, we consider a test problem with a non-polynomial analytical solution. The computational domain Ω is the unit square with boundary layers along x = 1 and y = 1. We are solving equation ( 1) with homogeneous Dirichlet boundary conditions, w = (1, 1) T and K = vI for v = 1 100 . The right hand side is given by f (x, y) := s v (x)g(y) + s v (y)g(x) + sin(πx) + sin(πy),

where

s v : [0, 1] → R is the solution of the equation -v[s v ] + [s v ] = 1 on [0, 1] and satisfying s v (0) = s v (1) = 0. The function g : [0, 1] → R is defined through g(z) := π[cos(πz) + πvsin(πz)].
In this setting, equation (1) admits the analytical solution u(x, y) = s v (x)sin(πy) + s v (y)sin(πx). We remark that for n = 25 and n = 50 method A is more accurate than its Hermite variant HA in the L 2 sense. However, for n = 100 the opposite occurs, and the expected second order convergence starts playing a role. The method A converges only linearly (in the L ∞ sense this is less clear up to this level of refinement). The H 1 errors are very big owing to the sharp boundary layers. The conclusion is that, for a moderate Péclet number the Hermite method HA is the best option, provided an affordable fine mesh is used. Of course if the Péclet number increases some stabilization is required, like for any other method, but addressing this issue is beyond the aim of this paper. Notice that the results for n = 25 indicate that the assumption on the magnitude of h for stability of method HA is not superfluous.

Concluding remarks

We summarize this work with a few conclusions and remarks.

Remark 1 The Hermite method described in Section 3 works as well as the corresponding Douglas and Roberts extension of the RT 0 element, as far as the fluxes are concerned. On the other hand the former behaves much better in terms of the error of the primal variable in L 2 (Ω).

Remark 2 According to the theoretical results derived in this work, numerical convergence in case the Péclet number is high could only be observed for meshes much finer than those used in Test-problem 1 and in [START_REF] Ruas | A Hermite finite element method for convection-diffusion equations[END_REF]. However running tests with such meshes may become unrealistic. Therefore the authors intend to study modifications of the variational formulations employed in this work, in order to obtain stable solutions within acceptable accuracy, even in the case of high Péclet numbers, without resorting to excessive mesh refinement. The work of Park and Kim [START_REF] Kim | A posteriori error estimators for the upstream weighting mixed methods for convection-diffusion problems[END_REF] for RT 0 discretizations could be a inspiring one in this connection. Remark 3 The bilinear forms and the linear form L h considered in this work do not really reduce to those in [START_REF] Ruas | Hermite finite elements for second order boundary value problems with sharp gradient discontinuities[END_REF] in the case where w ≡ 0. This is because somehow we wanted to incorporate numerical quadrature to the variational formulations in use, which is mandatory if f is not easy to integrate. However, in this case we should rather take L h (v) := (f h , v) for a suitable f h defined through point values of f only. Assuming for instance that f ∈ H 2 (Ω), f h can be chosen to be a piecewise linear interpolate of f in every T ∈ T h . Second order convergence results in L 2 (Ω) can still be proven to hold for such a choice, using the well-known analysis of variational crimes [START_REF] Strang | An Analysis of the Finite Element Method[END_REF] and [START_REF] Cuminato | Unification of distance inequalities for linear variational problems[END_REF]. The same qualitative results can also be obtained by using a suitable quadrature formula to compute the integral of the function g := f v in every element of the mesh. For more details we also refer to [START_REF] Strang | An Analysis of the Finite Element Method[END_REF], or to many other text books on the finite element method. Remark 4 The Hermite methods studied in this work can be viewed as a technique to improve the accuracy of the primal variable computations with mixed element RT 0 without resorting to post-processing (see e.g. [START_REF] Burman | Bubble stabilized discontinuous Galerkin method for parabolic and elliptic problems[END_REF], [START_REF] Lovadina | A posteriori error analysis of mixed finite element methods for second order elliptic equations[END_REF]) or hybridization (see e.g. [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] for the diffusion equation and [START_REF] Radu | Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study[END_REF] for convectiondiffusion problems). Incidentally the method proposed in [START_REF] Radu | Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study[END_REF] can be applied also to BDM elements to obtain optimal estimates [START_REF] Brunner | Optimal order convergence of a modified BDM1 mixed finite element scheme for reactive transport in porous media[END_REF] improving in this way the classical BDM method for convection-diffusion equations [START_REF] Demlow | Suboptimal and Optimal Convergence in Mixed Finite Element Methods[END_REF] or for stabilization purposes [START_REF] Radu | Accuracy of numerical simulations of contaminant transport in heterogeneous aquifers: a comparative study[END_REF][START_REF] Brunner | Analysis of an upwind-mixed hybrid finite element method for transport problems[END_REF] like in several previous work on the subject. Notice that our method allows to achieve better accuracy directly from the numerical solution procedure, at negligible additional cost. Thus it seems worthwhile searching for Hermite analogs of BDM methods as well in the future.

Remark 5 A comparison between the methods described in Section 3 (Method HA) and in Section 5 (Method HB), whose study was the main object of this work, advocates in favor of the former in all respects.

Remark 6 As a by-product of the numerical experimentation presented in Section 6 the two versions of the Douglas & Roberts extension of the RT 0 method were also compared to each other, in the framework of the solution of convection-diffusion equations at different Péclet numbers. We observed that the method in non divergence form (Method A) is substantially more accurate and reliable than the method in divergence form (Method B). To the best of authors' knowledge this kind of numerical study had not been carried out before.

  and b T ∈ N . Then the flux variable K∇v |T is of the form a T x + b T , and from a well-known property of the lowest order Raviart-Thomas mixed element a T and b T can be uniquely determined for prescribed V F (v) (resp U F (v)), ∀F ⊂ ∂T . Indeed by construction the flux variable for the Hermite element is locally defined by functions of the same form as for the lowest order Raviart-Thomas element. Once a T and b T are known, we determine the value of the additive constant d T to complete the expression of v |T , by enforcing the condition Π T [v] = 0. As for the basis function corresponding to the degree of freedom V T and U T , the values of a T and b T are obtained as follows: a T = 0 for both V T and U T while b T = -Kw h for V T and b T = 0 for U T . Then the value of d T is adjusted in such a way that the mean value of the corresponding quadratic function equals one. More precisely, for V T we have d T = Kw h • x T + 1 and for U T we have d T = 1.

  x 10 -3 0.15656954 x 10 -3 0.16950256 x 10 -3 0.24286177 x 10 -3

  0.13723841 x 10 -3 0.58218263 x 10 -3 0.15249263 x 10 -2 0.20428256 x 10 -5 100 0.13724039 x 10 -3 0.58595099 x 10 -3 0.25587661 x 10 -1 0.26841141 x 10 -5 10000 0.18370239 x 10 -2 0.51526514 x 10 +0 0.15093095 x 10 +3 0.52745011 x 10 -1 1000000 0.20738981 x 10 -3 0.57979017 x 10 -1 0.22150639 x 10 +2 0.12137294 x 10 -2

  , 7, 8 the same kind of results are displayed for Pé= 100.

	h	Method A	Method B	Method HA	Method HB
	1/8	0.10909896 x 10 -2 0.10917040 x 10 -2 0.17683545 x 10 -3 0.17991224 x 10 -3

Table 1 :

 1 Test-problem 1 with Pé=1 -L 2 errors of u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.45950033 x 10 -2 0.47137052 x 10 -2 0.46006556 x 10 -2 0.47021526 x 10 -2
	1/16 0.23211704 x 10 -2 0.23834158 x 10 -2 0.23219012 x 10 -2 0.23772504 x 10 -2

Table 2 :

 2 Test-problem 1 with Pé=1 -L 2 errors of ∇u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.12129518 x 10 -1 0.12201600 x 10 -1 0.12131302 x 10 -1 0.12201600 x 10 -1
	1/16 0.60914089 x 10 -2 0.61274572 x 10 -2 0.60916246 x 10 -2 0.61274572 x 10 -2
	1/32 0.30490236 x 10 -2 0.30670547 x 10 -2 0.30490504 x 10 -2 0.30670547 x 10 -2
	1/64 0.15249263 x 10 -2 0.15339429 x 10 -2 0.15249297 x 10 -2 0.15339429 x 10 -2

Table 3 :

 3 Test-problem 1 with Pé=1 -L 2 errors of ∆u for methods A, B, HA, HB From Tables 1 through 8 one can infer that:

Table 4 :

 4 Test-problem 1 with Pé=1 -Maximum errors of u at centroids for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.10964223 x 10 -2 0.37372813 x 10 -2 0.27608031 x 10 -3 0.45450302 x 10 -2
	1/16 0.54833181 x 10 -3 0.12935252 x 10 -2 0.46012246 x 10 -4 0.14346673 x 10 -2

Table 5 :

 5 Test-problem 1 with Pé=100 -L 2 errors of u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.59939985 x 10 -2 0.12889435 x 10 +0 0.77493931 x 10 -2 0.11984783 x 10 +0

Table 6 :

 6 Test-problem 1 with Pé=100 -L 2 errors of ∇u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.17192369 x 10 +0 0.24082822 x 10 +0 0.67758963 x 10 -1 0.24082822 x 10 +0

Table 7 :

 7 Test-problem 1 with Pé=100 -L 2 errors of ∆u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.54600580 x 10 -3 0.70291040 x 10 -2 0.55050954 x 10 -3 0.81134902 x 10 -2

Table 9 :

 9 Absolute errors for Test-problem 1 solved by Method A with h = 1/64

	Pé	e h	∇ h e h	∆ h e h	max T	|e h (x T )|
	1	0.28250216 x 10 -5 0.58219418 x 10 -3 0.15249297 x 10 -2 0.28130033 x 10 -5
	100	0.25386722 x 10 -5 0.59256341 x 10 -3 0.24402972 x 10 -1 0.37752993 x 10 -5
	10000	0.10089341 x 10 +4 0.26546637 x 10 +6 0.75214576 x 10 +7 0.57326794 x 10 +4
	1000000 0.13681695 x 10 +0 0.26569231 x 10 +2 0.66544379 x 10 +2 0.27070200 x 10 +1

Table 10 :

 10 Absolute errors for Test-problem 1 solved by Method HA with h = 1/64 Tables 9 and 10 indicate that Method A is more stable than its Hermite counterpart, as Pé increases.

  -2 0.12088608 x 10 -1 0.35534315 x 10 -3 0.82905447 x 10 -2 1/16 0.46131717 x 10 -2 0.91517998 x 10 -2 0.88950447 x 10 -4 0.80213268 x 10 -2 1/32 0.23075591 x 10 -2 0.82633266 x 10 -2 0.22251853 x 10 -4 0.79642825 x 10 -2 1/64 0.11539009 x 10 -2 0.80263174 x 10 -2 0.55709298 x 10 -5 0.79504439 x 10 -2

	h	Method A	Method B	Method HA	Method HB
	1/8	0.92106378 x 10

Table 11 :

 11 Test-problem 2 with Pé=1 -L 2 errors of u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.93229273 x 10 -8 0.49035715 x 10 -1 0.93226950 x 10 -8 0.49803379 x 10 -1
	1/16 0.93341685 x 10 -8 0.48942597 x 10 -1 0.93341123 x 10 -8 0.49096187 x 10 -1
	1/32 0.93369805 x 10 -8 0.48911400 x 10 -1 0.93369665 x 10 -8 0.48942095 x 10 -1
	1/64 0.93376835 x 10 -8 0.48902954 x 10 -1 0.93376800 x 10 -8 0.48909457 x 10 -1

Table 12 :

 12 Test-problem 2 with Pé=1 -L 2 errors of ∇u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.26390408 x 10 -7 0.26390408 x 10 -7 0.26390250 x 10 -7 0.26390408 x 10 -7

Table 13 :

 13 Test-problem 2 with Pé=1 -L 2 errors of ∆u for methods A, B, HA, HB

	h	Method A	Method B	Method HA	Method HB
	1/8	0.12068400 x 10 -3 0.22261595 x 10 -1 0.40130052 x 10 -3 0.25015968 x 10 -1

Table 15 for

 15 Method A and in Table 16 for Method HA. Pé e h Ω h ∇ h e h Ω h ∆ h e h Ω h max

	1	0.11539009 x 10 -2 0.93376835 x 10 -8 0.26411289 x 10 -7 0.18777283 x 10 -5
	100	0.11539009 x 10 -2 0.93376819 x 10 -8 0.26411285 x 10 -7 0.18777283 x 10 -5
	10000	0.11539009 x 10 -2 0.93377017 x 10 -8 0.26411346 x 10 -7 0.18777283 x 10 -5
	1000000 0.11539009 x 10 -2 0.96367921 x 10 -8 0.34526072 x 10 -7 0.18774905 x 10 -5

T |e h (x T )|

Table 15 :

 15 Absolute errors for Test-problem 2 solved by Method A with h = 1/64

	Pé	e h Ω h	∇ h e h Ω h	∆ h e h Ω h	max T	|e h (x T )|
	1	0.55709298 x 10 -5 0.93376800 x 10 -8 0.26411287 x 10 -7 0.63028673 x 10 -5
	100	0.55708924 x 10 -5 0.92319279 x 10 -8 0.26157933 x 10 -7 0.63028012 x 10 -5
	10000	0.55694727 x 10 -5 0.52231986 x 10 -8 0.15104352 x 10 -7 0.62996183 x 10 -5
	1000000 0.55709889 x 10 -5 0.93376916 x 10 -8 0.26406856 x 10 -7 0.63029895 x 10 -5

Table 16 :

 16 Absolute errors for Test-problem 2 solved by Method HA with h = 1/64

  The computations were done with uniform 2 × n × n meshes, for n = 25, n = 50 and n = 100,

	corresponding to a mesh diameter h =	√	2/25,	√ 2/50 and	√	2/100. The results are presented in Tables
	17-18 below.
	h 2/25 √ √ 2/50 √ 2/100 0.18041177 x 10 -1 0.20250530 x 10 -1 e h ∇ h (Ke h ) 0.88378563 x 10 -1 0.77159915 x 10 -1 0.37831058 x 10 -1 0.39786417 x 10 -1	∇ h • K∇ h e h 0.78473709 x 10 1 0.50624144 x 10 1 0.27814490 x 10 1	max T 0.38808847 x 10 0 |e h (x T )| 0.13186931 x 10 0 0.41014836 x 10 -1

Table 17 :

 17 Absolute errors for Test-problem 3 solved by Method A

	h 2/25 √ √ 2/50 √ 2/100 0.10370033 x 10 -1 0.23350942 x 10 -1 e h ∇ h (Ke h ) 0.27252781x 10 7 0.24839776 x 10 7 0.54136232 x 10 -1 0.73013269 x 10 -1	∇ h • K∇ h e h 0.21605994 x 10 9 0.91471443 x 10 1 0.45183509 x 10 0 max |e h (x T )| T 0.25252501 x 10 8 0.31967471 x 10 1 0.10196511 x 10 0

Table 18 :

 18 Absolute errors for Test-problem 3 solved by Method HA

Acknowledgment: The first author is thankful for the financial support provided by CNPq through grant 307996/2008-5 and the second author gratefully acknowledges the support of Statoil through the Akademia agreement.