Pierre Cantin

Alexandre Ern

An edge-based scheme on polyhedral meshes for vector advection-reaction equations

Keywords: classifications. 65N12, 65N15, 76R99, 76D07, 76W05 Keywords. Vector advection-reaction problems, polyhedral meshes, Friedrichs' assumptions, quasi-optimal a priori error estimates

come

An edge-based scheme on polyhedral meshes for vector advection-reaction equations

Introduction

Let Ω be a polyhedral domain of R d with d = 3 and consider a polyhedral mesh of Ω. We use boldface fonts for R d or R d×d -valued quantities. The purpose of this paper is to devise an approximation, using scalar degrees of freedom (dofs) attached to the edges of a mesh, of the R d -valued function u solving the vector advection-reaction problem: ∇(β•u) + (∇×u)×β + µu = s a.e. in Ω, (1a)

u = u D a.e. on ∂Ω -. (1b)
The R d -valued advective field β is assumed to be Lipschitz continuous in Ω and the R d×d -valued reaction tensor µ is assumed to be bounded in Ω. The subset ∂Ω -⊂ ∂Ω denotes the inflow part of the boundary where β•n < 0 with n the unit outward normal to Ω. The model problem [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF] is encountered in various situations. For example, it models the static advection of a magnetic field (u here) by a moving plasma of velocity β and of anisotropic conductivity µ. In the context of differential geometry, the operator ∇(β•u) + (∇×u)×β is the proxy of the Lie derivative of a differential 1-form (also called circulation) in R 3 (see Abraham et al. [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF] or Heumann [START_REF] Heumann | Eulerian and Semi-Lagrangian Methods for Advection-Diffusion of Differential Forms[END_REF]). The Lie derivative describes more generally the advection along the vector field β of a differential form on a manifold. The model problem [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF] is also relevant to study, in the advection-dominant regime, the advection-diffusion of a R d -valued field, which is one the building blocks of the Oseen problem or of the magneto-hydrodynamic problem. Using vector calculus rules, we observe that

∇(β•u) = (∇β) t u + (∇u) t β, (2a)
(∇×u)×β = (∇u)β -(∇u) t β, (2b)
where we have denoted ∇v the Jacobian matrix of v : Ω → R 3 such that its (i, j)-th component is ∂ j v i . As a result, combining the two above equations yields ∇(β•u) + (∇×u)×β = (∇u)β + (∇β) t u, so that the particular choice µ = -(∇β)

t yields the pure advection problem (with the more usual writing (∇u)β = (β•∇)u in this context):

(β•∇)u = s a.e. in Ω, (3a)

u = u D a.e. on ∂Ω -. (3b)
Edge-based schemes, that is, schemes using one scalar degree of freedom (dof) per mesh edge, are rarely addressed in the literature despite the fact that they are the natural way to discretize differential 1-forms, such as the electric field in electromagnetism or the flow velocity in fluid mechanics. For the Maxwell and the Stokes problem respectively, we mention for example the work of Zaglmayr [START_REF] Zaglmayr | High Order Finite Element Methods for Electromagnetic Field Computation[END_REF] and that of Girault [START_REF] Girault | The Navier-Stokes Equations Theory and Numerical Methods[END_REF] using Nédélec edge elements. In the context of our problem (1), Heumann and Hiptmair proposed in [START_REF] Heumann | Stabilized Galerkin for Transient Advection of Differential Forms[END_REF] an H(curl; Ω)conforming discretization of arbitrary order using Nédélec edge elements on simplicial meshes with a stabilization term in the spirit of the discontinuous Galerkin method (see Lesaint & Raviart [22], or Johnson & Pitkäranta [START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF]). In a different context and motivated by the discretization of the Lie derivative of a 1-form, we mention the Ph.D. thesis of Palha [START_REF] Palha | High order mimetic discretization[END_REF] approximating on square meshes a problem similar to [START_REF] Abraham | Manifolds, tensor analysis, and applications[END_REF] with the spectral element method. Based on the work of Bossavit [START_REF] Bossavit | Extrusion, contraction: their discretization via Whitney forms[END_REF], Mullen et al. also studied in [START_REF] Mullen | Discrete Lie advection of differential forms[END_REF] an approximation of (1) by extruding the edges of a simplicial mesh along the vector field β. All of the above schemes are devised on either simplicial or tensor-product meshes.

The first salient contribution of this work is to devise an edge-based scheme to approximate the model problem (1) on polyhedral meshes. The advantage of considering polyhedral meshes is multifold; it allows for more flexibility when meshing a complex geometry, it provides a natural framework to handle nonmatching mesh refinement and mesh coarsening by cell agglomeration, and it may even yield lower computational costs and better accuracy compared to the case of the simplicial meshes (see Bonelle's Ph.D. thesis [START_REF] Bonelle | Compatible Discrete Operator schemes on polyhedral meshes for elliptic and Stokes equations[END_REF]). The analysis framework for our scheme hinges on the notions of reduction and reconstruction maps as, e.g., in the mimetic approach of Kreeft et al. in [START_REF] Kreeft | Mimetic framework on curvilinear quadrilaterals of arbitrary order[END_REF], see also Gerritsma [START_REF] Gerritsma | An introduction to a compatible spectral discretization method[END_REF], or the Compatible Discrete Operator (CDO) approach of Bonelle & Ern [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF][START_REF] Bonelle | Analysis of compatible discrete operator schemes for stokes problems on polyhedral meshes[END_REF]. In particular, we consider a reconstruction map defining piece-wise constant vector-valued functions on an edge-based diamond partition of each mesh cell. This map has been introduced by Codecasa et al. in [START_REF] Codecasa | A new set of basis functions for the discrete geometric approach[END_REF] and has been recently revisited in the context of CDO schemes in [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF]. The novelty here is to perform the stability analysis in L q -spaces for q ∈ [1, ∞) and to prove a quasi-local consistency result by composing the reconstruction map on the right with a novel reduction map à la Clément that is stable for all integrable functions on a macro cell collecting all diamonds attached to the cell edges. This technique is key to establish an O(h q ) convergence rate as soon as the weak solution belongs to W 1,q (Ω) with q ∈ 3 2 , 2 without invoking a more stringent regularity assumption.

The second salient contribution of this work is to extend the well-posedness analysis at the discrete level to the non-coercive case. Specifically, we introduce an extended hypothesis on the problem coefficients (the fields β and µ) that allows one to go beyond the classical (and somewhat restrictive) assumption à la Friedrichs requiring the positivity of the minimal eigenvalue of the symmetric tensor

σ β,µ = (∇β + ∇β t ) -(∇•β)Id + (µ + µ t ) : Ω → R 3×3 .
Under this hypothesis, the well-posedness of the discrete problem classically hinges on a coercivity argument. However, this assumption is somehow restrictive; e.g. , the basic case of a constant vector field β with no reaction term does not fulfill this hypothesis. Motivated by our recent work [START_REF] Cantin | Vertex-based compatible discrete operator schemes on polyhedral meshes for advection-diffusion equations[END_REF] related to scalar advection-diffusion problem (see also the work of Deuring et al. in [START_REF] Deuring | L 2 -stability independent of diffusion for a Finite Element-Finite Volume discretization of a linear convection-diffusion equation[END_REF] for face-based finite volume schemes), we propose to extend the analysis to the non-coercive case where the minimal eigenvalue λ can take null or slightly negative values. Even if our analysis is presented here for our scheme, we emphasize that the main idea can be adapted to other schemes, such as Nédélec edge elements. We denote λ the minimal eigenvalue of σ β,µ over the domain Ω, i.e.

(Ω) = {v ∈ L 2 (Ω) | (β•∇)v ∈ L 2
(Ω)} if the fields β and µ satisfy the following hypothesis:

(H1) λ > 0. We define the reference time τ = λ -1 .

(H2) -C λ < λ ≤ 0, where C λ > 0 is a constant independent the mesh size, and there exists a potential

ζ ∈ W 1,∞ (Ω) satisfying ζ ≥ 1 and ess inf Ω (-β•∇ζ) > 0. We define the reference time τ = (ess inf Ω (-β•∇ζ)) -1 .
In the case of a continuously differentiable vector field β ∈ C 1 (Ω), the existence of the potential ζ is proved by Devinatz et al. in [START_REF] Devinatz | The asymptotic behavior of the first real eigenvalue of second order elliptic operators with a small parameter in the highest derivatives[END_REF]Lemma 2.3] by considering the Cauchy problem d t x(t) = β(x(t)), x(0) = x 0 ∈ Ω when the solution remains in the domain Ω for a finite time only. As a result, the hypothesis (H2) is satisfied if the vector field β has no closed curves and no stationary points in Ω. The analysis of our polyhedral edge-based scheme under this second hypothesis (H2) differs since the key idea is now to bound, at first-order in the mesh size, the commutator between the reconstruction map and the multiplication by the potential ζ. Using this technique, we can prove inf-sup stability (and infer the same convergence rates as above) as soon as the mesh size is smaller than a reference length that linearly depends on the quantity

||∇β t + µ|| -1 L ∞ (Ω) .
In particular, for the advective problem (3) (where µ ≡ -∇β t ), inf-sup stability holds with no restriction on the mesh size. This paper is organized as follows. In Section 2, we introduce the notation and the analysis tools on polyhedral meshes. In Section 3, we introduce the edge-based reconstruction map and we present the numerical scheme with dofs attached to edges. In Section 4, we state the main analytic results, namely, stability under hypothesis (H1) or (H2), boundedness and a priori error estimates delivering quasi-optimal decay rates for solutions in W 1,q (Ω) with q ∈ 3 2 , 2 . The proofs are postponed to Section 6 to facilitate the reading. Finally, we present in Section 5 numerical results on three-dimensionnal polyhedral meshes. A natural perspective for this work is to use the present scheme to discretize the advective operator in the Oseen (and Navier-Stokes) equations, while using the CDO scheme of [START_REF] Bonelle | Analysis of compatible discrete operator schemes for stokes problems on polyhedral meshes[END_REF] to discretize the Stokes operator in curl-curl formulation.

Notation and analysis tools on polyhedral meshes

We consider a general mesh M of Ω ⊂ R d with d = 3, composed of polyhedral cells c ∈ C (3-cells), planar faces f ∈ F (2-cells), straight edges e ∈ E (1cells), and vertices v ∈ V (0-cells). We collect the interior faces in the set F • = {f = ∂c ∩ ∂c | c = c and c, c ∈ C}, and we define F ∂ = F\F • the set collecting boundary faces. For any A, X ∈ {V, E, F, C}, we define the subset X a with a ∈ A as {x ∈ X | a ⊂ ∂x} if the dimension of a is smaller than that of the elements of X and as X a = {x ∈ X | x ⊂ ∂a} otherwise. For example, the set C e = {c ∈ C | e ⊂ ∂c} collects all the mesh cells containing the edge e, whereas the set E c = {e ∈ E | e ⊂ ∂c} collects all the mesh edges contained in the cell c, and so on. For any geometric entity x, we denote |x| its Hausdorff measure of appropriate dimension. In this paper, we assume mesh regularity in the sense that

• The mesh M := {V, E, F, C} defines a cellular complex (see Christiansen [START_REF] Christiansen | A construction of spaces of compatible differential forms on cellular complexes[END_REF]), i.e. the boundary of any k-cell, 1 ≤ k ≤ d (recall d = 3), is composed of a uniformly finite number of (k -1)-cells in M.

• Faces and cells are star-shaped with respect to their barycenters.

• Let x v denote the coordinates of v ∈ V in R d . Let x f and x c denote the coordinates of the barycenters of f ∈ F and c ∈ C, respectively, in R d . Then, the simplicial sub-mesh composed of the tetrahedra [x v , x v , x f , x c ] (where [x 1 , ..., x k+1 ] is the convex hull of the set {x 1 , ..., x k+1 }) for all c ∈ C, all f ∈ F c and all e ∈ E f with e = [x v , x v ] (see Figure 1, left panel) is shape-regular in the usual sense of Ciarlet.

For every cell c ∈ C, we introduce the edge-based diamond partition P c which plays a central role in our analysis. We define P c = ∪{p e,c ; e ∈ E c } where the diamond p e,c is defined by

p e,c = f ∈Fc∩Fe [x v , x v , x f , x c ] with e = [x v , x v ] ,
see Fig. 1, right panel. Note that P c is composed of #E c diamonds and that each diamond p e,c is composed of two tetrahedra, since #(F e ∩ F c ) = 2, with # faces are oriented by a fixed unit normal vector n f . For all f = ∂p e,c ∩∂p e ,c ∈ F c with e, e ∈ E c and n f pointing from p e,c to p e ,c , we define the jump and the average, respectively, as

• x v ù x e ◊ x f ¶ x c • x v • x v Õ ù x e ◊ x f ¶ x c • • ù x e ◊ ◊ ¶ x c 2 • x v • x v Õ ù x e ◊ x f ¶ x c • • ù x e ◊ ◊ ¶
• • ù x e ù x e Õ ◊ ◊ ◊ ¶ x c • • ù x e ◊ x f ◊ ¶ x c • • ù x e ◊ x f ◊ ¶ x c
[[v]] = v |pe,c -v |p e ,c and { {v} } := 1 2 v |pe,c + v |p e ,c .
Similarly, for all f = ∂p e,c ∩ ∂p e,c ∈ F f with c, c ∈ C f , e ∈ E f , and n f pointing from p e,c to p e,c , we define

[[v]] := v |pe,c -v |p e,c and { {v} } := 1 2 v |pe,c + v |p e,c .
We denote |•| 2 the Euclidean and the Frobenius norm on R d and R d×d , respectively. For every set ω ⊂ Ω, we denote

L q (ω) with q ∈ [1, ∞] the Banach space of R d or R d×d -valued functions v such that ||v|| L q (ω) := || |v| 2 || L q (ω) < ∞.
Lemma 2.1 (Mutliplicative trace inequality). There exists C T > 0 such that

||v|| L q (f) ≤ C T ||v|| 1-1 q L q (p) h -1 q c ||v|| 1 q L q (p) + |v| 1 q W 1,q (p) , (4)
for all c ∈ C with h c the diameter of c, all p ∈ P c , all f ∈ ∂p and all v ∈ W 1,q (p)

with q ∈ [1, ∞].
Proof. Observing that p ⊂ P c is composed of two tetrahedra connected by a sub-face f ∈ F c , this result follows proceeding as in Ern & Guermond [START_REF] Ern | Finite element quasi-interpolation and best approximation[END_REF].

3 Discrete Scheme

Degrees of freedom

We consider an approximation of the continuous problem (1) with scalar dofs attached to edges. We denote E ≡ R #E the linear space collecting these dofs and we denote v e the entry of v ∈ E attached to the edge e ∈ E. We additionally introduce the linear space E c collecting the dofs attached to the subset E c for all c ∈ C. We denote v a generic element of E or E c .

Reconstruction map

The global reconstruction map L E is defined locally, so that L E (v) |c = L Ec (v), for all c ∈ C. The local reconstruction map L Ec : E c → P 0 (P c ; R d ), where P 0 (P c ; R d ) is composed of piece-wise constant R d -valued polynomials over the diamond partition P c , is such that

L Ec (v)(x) := e∈Ec v e e,c (x), ∀v ∈ E c , ∀x ∈ c, (5)
where for all e ∈ E c , the basis function e,c ∈ P 0 (P c ; R d ), is defined by

and δ e,e is the Kronecker symbol equal to 1 if e = e and 0 otherwise. Moreover, for all e ∈ E, t e is a fixed unit tangent vector to e, such that e = |e|t e , and fc (e) = fc(e) n fc(e) where the dual face f (e) is composed of two elementary

triangles fc (e) = f ∈Fc∩Fe [x e , x f , x c ] ,
see Fig. 3, and where n fc(e) is the unit normal vector to fc (e) satisfying n fc(e) •t e ≥ 0. The basis functions e,c were first considered in the context of the Discrete Geometric Approach by Codecasa et al. [START_REF] Codecasa | A new set of basis functions for the discrete geometric approach[END_REF] and were recently revisited by Bonelle & Ern in [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF][START_REF] Bonelle | Analysis of compatible discrete operator schemes for stokes problems on polyhedral meshes[END_REF] to build Hodge operators within the CDO framework. They satisfy the following properties:

• • ù x e ◊ ◊ ¶ x c • x v ù ù ù ◊ ◊ ◊ ¶ x c • • ù x e ◊ ◊ ¶ x c • • • x v ◊ x f ¶ x c 3

Discrete scheme

The discrete scheme is formulated using the global bilinear form

A β,µ : E×E → R such that A β,µ (u, v) = A β,µ (u, v) + A ∂ (β•n) -(u, v), (7)
where A β,µ approximates (1a) and A ∂ (β•n) -weakly enforces the boundary condition (1b). The bilinear form A β,µ : E×E → R is composed of three bilinear forms also defined on E×E:

A β,µ (u, v) := g β,µ (u, v) + n β (u, v) + s β (u, v). (8)
The bilinear form g β,µ is assembled cell-wise as

g β,µ (u, v) = c∈C g β,µ;c (u, v), (9)
and each local bilinear form g β,µ;c results from the standard Galerkin approximation of (1a) in c using the reconstruction map L Ec :

g β,µ;c (u, v) = p∈Pc p (∇(β•L Ec (u)) + (∇×L Ec (u))×β)•L Ec (v)+ c µL Ec (u)•L Ec (v).
(10) Using the identities (2) and since L Ec (v) is piece-wise constant, we can reformulate this expression as

g β,µ;c (u, v) = c ∇β t + µ L Ec (u) • L Ec (v). ( 11
)
Because L Ec (v) jumps across inter-cell and intra-cell sub-faces, we also consider the bilinear form n β such that

n β (u, v) = c∈C n β;c (u, v) + f ∈F • n β;f (u, v), (12)
where the local bilinear forms n β;x with x = f or x = c are defined as

n β;x (u, v) = - f∈Fx f (β•n f )[[L E (u)]]•{ {L E (v)} }, (13)
and the stabilization bilinear form s β such that

s β (u, v) = c∈C s β;c (u, v) + f ∈F • s β;f (u, v), (14)
where the local bilinear forms s β;x with x = f or x = c are defined as

s β;x (u, v) = f∈Fx f |β•n f |[[L E (u)]]•[[L E (v)]]. (15)
The bilinear forms n β and s β are devised similarly to the discontinuous Galerkin method; n β corresponds to centered fluxes and n β +s β to upwind fluxes. Finally, the Dirichlet boundary condition is weakly enforced by means of the bilinear form

A ∂ α : E × E → R (with α = (β•n) -) such that A ∂ α (u, v) = f ∈F ∂ A ∂ α;f (u, v). (16)
The local bilinear form A ∂ α;f is defined as

A ∂ α;f (u, v) = f α L Ec f (u)•L Ec f (v), (17)
with c f is the unique cell containing the boundary face f . The discrete scheme consists in finding u ∈ E such that

A β,µ (u, v) = Σ(s, u D ; v), ∀v ∈ E, (18)
with the right-hand side form Σ(s,

u D ; •) : E → R such that Σ(s, u D ; v) := c∈C c s•L Ec (v) + f ∈F ∂ f (β•n) -u D •L Ec f (v). ( 19
)
4 Stability and error analysis

|||v||| q,c ≤ ||L Ec (v)|| L q (c) ≤ C |||v||| q,c , for all c ∈ C, all v ∈ E c , all q ∈ [1, ∞)
and where

|||v||| q,c = e∈Ec |p e,c | |e| q |v e | q 1 q . ( 20
)
Remark 4.2 (Alternative definition). In lieu of (20), we could also consider the simpler discrete L q -norm given by |||v||| q q,c = h d-q c e∈Ec |v e | q . Owing to mesh regularity, this definition is equivalent to [START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF] up to a uniform constant with respect to the mesh-size. We prefer to use [START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF] since it simplifies the proof of Proposition 4.1.

We introduce the reduction map

R E : L 1 (Ω) → E such that R E (v)| e := 1 |p e | pe v•e , ∀e ∈ E, ( 21
)
where p e = ∪{p e,c ; c ∈ C e } is the diamond volume surrounding the edge e and ĉ is the local diamond patch ĉ = ∪{p e ; e ∈ E c } surrounding the cell c; notice that c ĉ. We also define the local reduction map R Ec : L 1 (ĉ) → E c from definition [START_REF] Kreeft | Mimetic framework on curvilinear quadrilaterals of arbitrary order[END_REF] for all e ∈ E c . e v•e for every e ∈ E can be used as well, provided that v ∈ H 1+ (Ω) [START_REF] Ern | Finite element quasi-interpolation and best approximation[END_REF]

or v ∈ {w ∈ H 1 2 + (Ω), ∇×w ∈ L 2+ (Ω)} [2] with > 0.
For each cell c ∈ C, we denote I Ec the local interpolation operator obtained by composing the local reconstruction map with the local reduction map, i.e.

I Ec = L Ec • R Ec , so that I Ec : L 1 (ĉ) → P 0 (P c ; R d ).
Proposition 4.4 (Consistency). For all c ∈ C and all U ∈ P 0 (ĉ; R d ) (so that U is a constant function in ĉ), we have I Ec (U ) = U |c .

Lemma 4.5 (Interpolation error).

There exists C Int > 0 such that for all c ∈ C and all v ∈ W 1,q (ĉ)

with q ∈ [1, ∞), ||v -I Ec (v)|| L q (c) ≤ C Int h c |v| W 1,q (ĉ) , (22)
and for all p ∈ P c ,

||v -I Ec (v)|| L q (∂p) ≤ C Int h 1-1 q c |v| W 1,q (ĉ) . (23)

Well-posedness under (H1)

We consider the following stability norm on the edge dof space E:

|||v||| := τ -1 |||v||| 2 2 + |v| 2 ∂ + |v| 2 s 1 2 , ( 24
)
where the reference time τ > 0 is defined by assumption (H1) or (H2),

|||•||| 2 2 = c∈C |||•||| 2 2,c is the discrete L 2 -norm with |||•||| 2,c defined by (20), |•| 2 ∂ = A ∂ |β•n| (•, •) is the semi-norm
1 2 |||v||| 2 ≤ A β,µ (v, v), ∀v ∈ E.
Consequently, the discrete problem (18) is well-posed.

Well-posedness under (H2)

In this section, we address the stability of the bilinear form A β,µ under the hypothesis (H2). We consider the reference length

h -1 0 = 4C 2 L ζ τ ||µ+∇β t || L ∞ (Ω)
, where C results from Proposition 4.1 and

L ζ = |ζ| W 1,∞ (Ω) is the Lipschitz con- stant of ζ. If µ = -∇β t ,
we conventionally set h 0 = +∞. Recalling that λ denotes the smallest eigenvalue of the tensor σ β,µ over the domain Ω, we assume that 1 + 2ϑτ λ > 0 and h < h

0 (1 + 2ϑτ λ ) , (25)
where ϑ > 0 is a non-dimensional constant that linearly depends on ||ζ|| A β,µ (v, w), ∀v ∈ E.

L ∞ (Ω) + C T C L ζ max(|Ω| 1 d , ||β|| L ∞ (Ω) τ ). By convention, the second condition in (25) is void if µ = -∇β t .
Consequently, the discrete problem (18) is well-posed.

In the proof of Proposition 4.7, the idea is to introduce a discrete test function ζv ∈ E defined as (ζv) e = ζ(x e )v e for all v ∈ E and for all e ∈ E. The key argument to obtain the well-posedness of the discrete problem [START_REF] Heumann | Eulerian and Semi-Lagrangian Methods for Advection-Diffusion of Differential Forms[END_REF] under hypothesis (H2) is then to bound the commutator δ defined as

δ(v) |c = L Ec (ζv) -ζL Ec (v), ∀v ∈ E c , ∀c ∈ C. ( 26
)
Lemma 4.8 (Bounds on δ). For all c ∈ C, we have

||δ(v)|| L 2 (c) ≤ 2C L ζ h c |||v||| 2,c , ∀v ∈ E c . ( 27a
)
and for all f ∈ F c ,

||δ(v)|| L 2 (f ) ≤ 2C T C L ζ h 1 2 c |||v||| 2,c , ∀v ∈ E c . (27b)
Table 1 recapitulates the different situations where the discrete problem ( 18) is well-posed.

λ > 0 -1 2ϑτ < λ ≤ 0 (H1) (H2) µ = -∇β t µ = -∇β t h ∈ R >0 h ∈ R >0 h ∈ (0, h 0 (1 + 2ϑτ λ ))
Table 1: Stability of the discrete problem [START_REF] Heumann | Eulerian and Semi-Lagrangian Methods for Advection-Diffusion of Differential Forms[END_REF] with respect to λ and the meshsize h.

Bound on consistency error and a priori error estimate

In this section, we derive an a priori error estimate by bounding the consistency error

E(u) = sup v∈E,|||v|||=1 |A β,µ (R E (u), v) -Σ(s, u D ; v)| .
In what follows, the notation A B stands for A ≤ CB where C is a positive constant uniform with respect to the mesh size and the model parameters.

Lemma 4.9 (Bound on consistency error). Assume that the exact solution satisfies u ∈ W 1,q (Ω) with q ∈ [1, 2]. Then, the following holds:

E(u) c∈C ||∇β + µ t -(∇•β)Id || q L ∞ (c) τ q 2 h d 2 (q-2) c ||u -I Ec (u)|| q L q (c) 1 q +   c∈C p∈Pc ||β|| q 2 L ∞ (c) h (d-1) 2 (q-2) c ||u -I Ec (u)|| q L q (∂p)   1 q .
We can now state the main result of this paper which follows from Lemmata 4.5 and 4.9.

Theorem 4.10 (A priori estimate). Assume that the assumptions stated in Table 1 hold. Assume that the exact solution of (1) satisfies u ∈ W 1,q (Ω) with q ∈ 2d d+1 , 2 . Then, we have

|||u -R E (u)||| c∈C ||∇β + µ t -(∇•β)Id || q L ∞ (c) τ q 2 h d+2 2 (q-2d d+2 ) c |u| q W 1,q (ĉ) 1 q +   c∈C p∈Pc ||β|| q 2 L ∞ (c) h d+1 2 (q-2d d+1 ) c |u| q W 1,q (ĉ)   1 q . For d = 3, it follows that |||u -R E (u)||| = O h 2-3 q for all q ∈ 3 2 , 2 .

Numerical results

We investigate numerically the edge-based scheme (18) on four sequences of three-dimensional polyhedral meshes. Each mesh is obtained as a uniform refinement of an initial mesh. Meshes from the first sequence, denoted H, are composed of hexahedra, those from the second one, denoted PrT, are composed of prisms with a triangular basis, those from the third one, denoted PrG, are composed of prisms with a hexagonal basis, and those of the last one, denoted CB, are composed of hexahedra with non-matching interfaces; see Figure 4. The domain is the unit cube Ω := [0, 1] 3 . The exact solution corresponds to a Taylor-Green velocity field, the advective vector field β is affine (see Figure 5, left panel) and the reaction tensor µ is diagonal and constant: u =   sin(πx) cos(πy/2) cos(πz/2) cos(πx/2) sin(πy) cos(πz/2) cos(πx/2) cos(πy/2) sin(πz)

  , β = 1 2   (x -2y)/2 (y -2x)/2 -z   , µ = 1 2
Id. Note that ∇•β = 0 and that the eigenvalues of the tensor σ β,µ are {0, 1 2 , 5 2 }, so that the discrete scheme ( 18) is well-posed owing to Proposition 4.7 if the mesh size is small enough.

We perform a convergence study by computing the relative discrete L 2 -error attached to edge dofs, denoted Er E (u), and defined by

Er E (u) = |||u -R E (u)||| 2 |||R E (u)||| 2 ,
with the norm |||•||| 2 on every cell of the mesh by ( 20). The convergence rates, shown in the right panel of 5, lie between 1 2 and 1 for the PrT and PrG mesh sequences and are closer to 1 for the H and CB mesh sequences. Note that the considered meshes being quasi-uniform, we have h ∼ (#E) -1/3 ; the reference slopes indicated in Figure 5 are based on this scaling, i.e., are with respect to h. Table 6 provides additional information on the computational costs by reporting the size of the linear system (#E), the mean stencil St, the values of the discrete error Er E (u), and the ratios #E/#V and #E/#C, indicating that the present scheme may involve less dofs than traditional Finite Volume schemes placing R d -valued unknowns at mesh vertices or at mesh cells. Note that owing to the Euler-Poincaré characteristic formula (in dimension d = 3; see e.g. [START_REF] Ziegler | Lectures on polytopes[END_REF]Chapter 8]),

#V #E + #F #E -#C #E = 2 #E -1 ≈ 1. #E St Er E (u) #E #V #E #C
3.0e+02 21

Recalling that ||•||

L q (c) = |||•| 2 || L q , we infer that ||L Ec (v)|| q L q (c) = e∈Ec |v e a e + b e | 2 q L q (pe,c)
.

Using the Property ( 1), we observe that a e •b e ≡ 0 on p e,c , so that |v e a e + b e | 2 ≥ |v e a e | 2 , whence

||L Ec (v)|| q L q (c) ≥ e∈Ec |v e a e | 2 q L q (pe,c) = e∈Ec |v e | q ||a e || q L q (pe,c) .
Hence, the expected lower bound follows from ||a e || q L q (pe,c) = |pe,c| |e| q . (ii) Upper bound. The discrete Hölder inequality yields

||L Ec (v)|| q L q (c) ≤ (#E c ) q-1 e∈Ec |v e | q || e,c || q L q (c) .
Since || e,c || q L q (c) ≤ |c||| e,c || q L ∞ (c) , we have || e,c || q L q (c) ≤ C q (#E c ) 1-q |pe,c| |e| q with the constant

C = (#E c ) 1-1 q max e∈Ec |c| |p e,c | 1 q |e||| e,c || L ∞ (c) ,
that is uniformly bounded owing to mesh regularity, yielding the expected upper bound. Specifically, a straightforward calculation shows that Proof of Lemma 4.5. Let c ∈ C and let v ∈ W 1,q (ĉ) with q ∈ [1, ∞). Owing to the triangle inequality and the P 0 -consistency of the reconstruction map from Proposition 4.4, we infer that

||v -I Ec (v)|| L q (c) ≤ ||v -v ĉ|| L q (c) + ||I Ec (v -v ĉ)|| L q (c) with v ĉ = |ĉ| -1 ĉ v.
In addition, we observe that, for all w ∈ L q (ĉ),

|||R Ec (w)||| q q,c = e∈Ec |p e,c | |e| q 1 |p e | pe w•e q ≤ e∈Ec |p e,c | |p e | q ||w|| q L 1 (pe) ≤ e∈Ec 1 |p e | q-1 ||w|| q L 1 (pe) ,
where we have used that |p e,c | ≤ |p e | to infer the last inequality. Owing to the Hölder inequality, it then follows that ||w|| q L 1 (pe) ≤ ||w|| q L q (pe) ||1|| q L q (pe) with

1 q + 1 q = 1. Since ||1|| q L q (pe) = |p e | q-1 , we infer that |||R Ec (w)||| q q,c ≤ ||w|| q L q (ĉ) .
Using this estimate and the upper bound from Proposition 4.1, we obtain

||I Ec (v -v ĉ)|| L q (c) ≤ C |||R Ec (v -v ĉ)||| q,c ≤ C ||v -v ĉ|| L q (c) , so that ||v -I Ec (v)|| L q (c) ≤ (1 + C )||v -v ĉ|| L q (c) . Hence, ||v -I Ec (v)|| L q (c) ≤ (1 + C )φ ĉ,q h c |v| W 1,q (ĉ) with φ ĉ,q = sup w∈W 1,q (ĉ) ||w -w ĉ|| L q (ĉ) h c |w| W 1,q (ĉ) .
Finally, we observe that the diamond ĉ can be decomposed as

ĉ = e∈Ec p e = e∈Ec c∈Ce p e,c ,
where p e,c consists of two tetrahedra, so that ĉ is composed of 2 e∈Ec #C e tetrahedra connected through elements of F c and F f with f ∈ F c . Then, proceeding as Ern & Guermond in [15, Lemma 5.7], we infer that the quantity φ ĉ,q is uniformly bounded for all c ∈ C and all q ∈ [1, ∞).

Well-posedness under (H1)

Proof of Proposition 4.6. Let c ∈ C and consider v, w ∈ E c . The definition of the bilinear form g β,µ;c together with the definition of the tensor σ β,µ yield

g β,µ;c (v, w)+g β,µ;c (w, v) = c L Ec (v)•σ β,µ •L Ec (w)+ p∈Pc p ∇•(β L Ec (v)•L Ec (w)).
(28) Choosing w = v in this relation leads to

1 2 p∈Pc p ∇•(β|L Ec (v)| 2 2 )- f∈Fc f (β•n f )[[L E (v)]]•{ {L E (v)} } = 1 2 f ∈Fc f (β•n c )|L Ec (v)| 2 2 ,
with n c the unit outward normal vector to c, so that recalling the definition (13) of n β;c , we infer that

c∈C (g β,µ;c (v, v) + n β;c (v, v)) = 1 2 c∈C c L Ec (v)•σ β,µ •L Ec (v)+ 1 2 c∈C f ∈Fc f (β•n c )|L Ec (v)| 2 2 .
The above rightmost term is reformulated as

1 2 c∈C f ∈Fc f (β•n c )|L E (v)| 2 2 = 1 2 f ∈F ∂ f (β•n)|L Ec f (v)| 2 2 + f ∈F • f∈F f f (β•n f )[[L E (v)]]•{ {L E (v)} } = 1 2 f ∈F ∂ f (β•n)|L Ec f (v)| 2 2 + f ∈F • n β;f (v, v)
so that, using the definition [START_REF] Devinatz | The asymptotic behavior of the first real eigenvalue of second order elliptic operators with a small parameter in the highest derivatives[END_REF] of n β , we arrive at

g β,µ (v, v) + n β (v, v) = 1 2 c∈C c L Ec (v)•σ β,µ •L Ec (v) + 1 2 f ∈F ∂ f (β•n)|L Ec f (v)| 2 2 .
Recalling the definition (7) of A β,µ and combining the above relation with the bilinear forms s β and A ∂ (β•n) -, defined by ( 14) and ( 16) respectively, we obtain

A β,µ (v, v) = 1 2 c∈C c L Ec (v)•σ β,µ •L Ec (v) + 1 2 A ∂ |β•n| (v, v) + 1 2 s β (v, v). (29)
The expected result is inferred from (H1) .

Well-posedness under (H2)

Proof of Lemma 4.8. Let v ∈ E and let c ∈ C.

(i) Proof of (27a). Let ζ c be the mean-value of ζ over c given by

ζ c = |c| -1 c ζ. Since L Ec (ζ c v) = ζ c L Ec (v) because ζ c is constant, we have δ(v) |c = (ζ- ζ c )L Ec (v) -L Ec ((ζ -ζ c )v)
, so that the triangle inequality, the Hölder inequality and the upper bound in Proposition 4.1 yield

||δ(v)|| L 2 (c) ≤ ||ζ -ζ c || L ∞ (c) ||L Ec (v)|| L 2 (c) + ||L Ec ((ζ -ζ c )v)|| L 2 (c) ≤ C ||ζ -ζ c || L ∞ (c) |||v||| 2,c + C |||(ζ -ζ c )v||| 2,c ≤ 2C ||ζ -ζ c || L ∞ (c) |||v||| 2,c . Observing that ||ζ -ζ c || L ∞ (c) ≤ L ζ h c , the expected result follows.
(ii). Proof of (27b). Let p ∈ P c and let f ⊂ ∂p. Owing to the multiplicative trace inequality (4), we have

||δ(v)|| L 2 (f) ≤ C T ||δ(v)|| 1 2 L 2 (p) h -1 2 c ||δ(v)|| 1 2 L 2 (p) + |δ(v)| 1 2 H 1 (p) . Observe that |δ(v)| H 1 (p) = |||∇δ(v)| 2 || L 2 (p) where |∇δ(v)| 2 2 = d i,j |∂ j δ(v) i | 2
in the Cartesian basis of R d and where ∂ i is the weak derivative in the direction i. Since L Ec (v) is piece-wise constant on P c , it then follows that |∇δ(v)|

2 2 = d i,j |L Ec (v) i ∂ j ζ| 2 = |L Ec (v)| 2 2 |∇ζ| 2 2 . As a result, |δ(v)| H 1 (p) ≤ L ζ ||L Ec (v)|| L 2 (p) . Moreover, proceeding as in (i), we infer that ||δ(v)|| L 2 (p) ≤ 2L ζ h c |||L Ec (v)||| L 2 (p) .
Collecting these bounds, we infer that

||δ(v)|| L 2 (f) ≤ 2C T L ζ h 1 2 c ||L Ec (v)|| L 2 (p) .
Then, summing over F f and using the upper bound of Proposition 4.1 yield the expected result.

In what follows, we consider the non-dimensional number

ω ζ = L ζ max(|Ω| 1 d , ||β|| L ∞ (Ω) τ ).
|||ζv||| 2 = c∈C τ -1 |||ζv||| 2 2,c + c∈C s β,c (ζv, ζv) + f ∈F ∂ A |β•n|;f (ζv, ζv) + f ∈F • s β,f (ζv, ζv) = T 1 + T 2 + T 3 + T 4 .
We want to use the Lipschitz regularity of ζ to bound separately these terms by |||v||| 2 . We recall the notation ζ c = |c| -1 c ζ from the proof of Lemma 4.8.

(i) Bound on T 1 . First, the triangle inequality implies that

1 2 T 1 ≤ c∈C τ -1 |||ζ c v||| 2 2,c + c∈C τ -1 |||(ζ -ζ c )v||| 2 2,c = T 1,1 + T 1,2 .
Since

|ζ c | ≤ ||ζ|| L ∞ (c) , we infer that T 1,1 ≤ c∈C τ -1 ||ζ|| 2 L ∞ (c) |||v||| 2 2,c ≤ ||ζ|| 2 L ∞ (Ω) |||v||| 2 . The bound T 1,2 easily follows from the Lipschitz regularity of ζ since T 1,2 ≤ c∈C τ -1 L 2 ζ h 2 c |||v||| 2 2,c ≤ L 2 ζ h 2 c |||v||| 2 . Combining these two bounds with h c ≤ |Ω| 1 d
and the definition of ω ζ yields

T 1 ≤ 2 ω 2 ζ + ||ζ|| 2 L ∞ (Ω) |||v||| 2 .
(ii) Bound on T 2 . Since the bilinear form s β;c is symmetric and positive, we infer that

1 2 T 2 ≤ c∈C s β;c (ζ c v, ζ c v) + c∈C s β;c ((ζ -ζ c )v, (ζ -ζ c )v) = T 2,1 + T 2,2 ,
and we have directly that

T 2,1 ≤ c∈C ||ζ|| 2 L ∞ (c) s β,c (v, v) ≤ ||ζ|| 2 L ∞ (Ω) |||v||| 2 .
To bound T 2,2 , we use the multiplicative trace inequality (4) and that L Ec is piecewise constant to infer that

s β;c ((ζ-ζ c )v, (ζ-ζ c )v) = f∈Fc f |β•n f | |[[L Ec ((ζ -ζ c )v)]]| 2 2 ≤ 2C 2 T ||β|| L ∞ (c) f∈Fc p∈Pc∩P f h -1 c ||L Ec ((ζ-ζ c )v)|| 2 L 2 (p)
where P f = {p ∈ P | f ⊂ ∂p}. Observing that the boundary of each diamond p e,c is composed of 4 sub-faces in F c , exchanging the sums yields

s β;c ((ζ-ζ c )v, (ζ-ζ c )v) ≤ 8C 2 T ||β|| L ∞ (c) p∈Pc h -1 c ||L Ec ((ζ-ζ c )v)|| 2 L 2 (p) = 8C 2 T ||β|| L ∞ (c) h -1 c ||L Ec ((ζ-ζ c )v)|| 2 L 2 (c) ,
Owing to upper bound from Proposition 4.1, the Lipschitz regularity of ζ, and the definition of ω ζ , we infer that

s β;c ((ζ-ζ c )v, (ζ-ζ c )v) ≤ 8C 2 T C 2 ||β|| L ∞ (c) h -1 c |||(ζ-ζ c )v)||| 2 2,c ≤ 8C 2 T C 2 ω 2 ζ τ -1 |||v||| 2 2,c .
Finally, collecting these two bounds leads to

T 2 ≤ 2 ||ζ|| 2 L ∞ (Ω) + 8C 2 T C 2 ω 2 ζ |||v||| 2 .
(iii) Bound on T 3 . We proceed as in the previous step (ii) to infer that

T 3 ≤ 2 ||ζ|| 2 L ∞ (Ω) + n F,∂ C 2 T C 2 ω 2 ζ |||v||| 2 ,
where n F,∂ = (max c∈C #(F c ∩ F ∂ )) is the naximal number of boundary faces that a mesh cell can have.

(iv) Bound on T 4 . To bound this last term, we use a different decomposition, namely

T 4 = f ∈F • s ζ 2 β;f (v, v) + f ∈F • ∆ f (v) = T 4,1 + T 4,2 , with ∆ f (v) = s β;f (ζv, ζv)-s ζ 2 β;f (v, v). Observing that s ζ 2 β;f (v, v) ≤ ||ζ|| 2 L ∞ (f ) s β;f (v, v) for all f ∈ F • , it follows that T 4,1 ≤ ||ζ|| 2 L ∞ (Ω) |||v||| 2 .
To bound the second term T 4,2 , we recall the quantity δ(v) defined by [START_REF] Ziegler | Lectures on polytopes[END_REF] and we obtain

∆ f (v) = f |β•n f | [[δ(v)]] + ζ[[L E (v)]] 2 -ζ 2 [[L E (v)]] 2 .
Then, applying Young's inequality and the trace inequality (27b) yields

|∆ f (v)| ≤ 2 f |β•n f | [[δ(v)]] 2 + f |β•n f | ζ 2 [[L E (v)]] 2 ≤ 4(2C T C L ζ ) 2 ||β|| L ∞ (f ) c∈C f h c |||v||| 2 2,c + ||ζ|| 2 L ∞ (f ) s β;f (v, v).
As a result, since #C f = 2 for all f ∈ F • and introducing ω ζ , we infer that

T 4,2 ≤ 32C 2 T C 2 ω 2 ζ τ -1 |||v||| 2 2 +||ζ|| 2 L ∞ (Ω) f ∈F • s β;f (v, v) ≤ 32C 2 T C 2 ω 2 ζ + ||ζ|| 2 L ∞ (Ω) |||v||| 2 , whence T 4 ≤ 2 16C 2 T C 2 ω 2 ζ L ζ h + ||ζ|| 2 L ∞ (Ω) |||v||| 2 .
(v) Conclusion. The expected inequality then follows from the above four bounds.

Proof of Proposition 4.7. Let v ∈ E and define

S = sup w∈E\{0} A β,µ (v, w) |||w||| .
Let us take w = ζv+θv with θ > 0 to be chosen below. We infer from Lemma 6.1 that

A β,µ (v, w) ≤ S|||w||| ≤ S θ + C ζ ||ζ|| L ∞ (Ω) + ω ζ |||v|||,
so that it remains to prove that A β,µ (v, w) |||v||| 2 . First, we split A β,µ as follows:

A β,µ (v, w) = A β,-∇β t + 1 2 (∇•β)Id (v, w) + H µ+∇β t -1 2 (∇•β)Id (v, w) = T 1 + T 2 , where the bilinear form H α is defined on E×E by H α (v, w) = c∈C c αL Ec (v) • L Ec (w),
for all α ∈ L ∞ (Ω). Let us bound from below the two terms T 1 and T 2 .

(i) Bound on T 1 . We bound from below this term by considering the following decomposition

T 1 = A β,-∇β t + 1 2 (∇•β)Id (v, w) = A ζβ,-∇(ζβ) t + 1 2 ζ(∇•β)Id (v, v) + A β,-∇β t (v, ζv) -A ζβ,-∇(ζβ) t (v, v) + H 1 2 (∇•β)Id (v, ζv) -H 1 2 ζ(∇•β)Id (v, v) + θA β,-∇β t + 1 2 (∇•β)Id (v, v) = T 1,1 + T 1,2 + T 1,3 + T 1,4 . Regarding T 1,1 , we use the relation (29) to infer that T 1,1 = 1 2 c∈C c L Ec (v)•σ ζβ,-∇(ζβ) t + 1 2 ζ(∇•β)Id •L Ec (v) + A ∂ |ζβ•n| (v, v) + s ζβ (v, v) ≥ 1 2 c∈C c L Ec (v)•σ ζβ,-∇(ζβ) t + 1 2 ζ(∇•β)Id •L Ec (v) + A ∂ |β•n| (v, v) + s β (v, v) , since ζ ≥ 1. Then, observing that σ ζβ,-∇(ζβ) t + 1 2 ζ(∇•β)Id = -β•∇ζ
Id and using hypothesis (H2) together with the lower bound from Proposition (4.1), we infer that T 1,1 ≥ 1 2 |||v||| 2 . The next step consists in bounding the perturbation term T 1,2 . To do so, we recall the identity (11) for g β,µ;c , and we observe that g β,-∇β t ;c ≡ 0 and g ζβ,-∇(ζβ) t ;c ≡ 0, so that T 1,2 solely consists of surfacic terms:

T 1,2 = n β (v, ζv) -n ζβ (v, v) + s β (v, ζv) -s ζβ (v, v) + A ∂ (β•n) -(v, ζv) -A ∂ (ζβ•n) -(v, v) .
Now, introducing the function δ(v) locally defined by [START_REF] Ziegler | Lectures on polytopes[END_REF] and recalling that

β ∈ W 1,∞ (Ω), ζ ≥ 1, and ζ ∈ W 1,∞ (Ω), , so that ζ{ {L Ec (v)} } = { {ζ L Ec (v)} }, we observe that n β;x (v, ζv) -n ζβ;x (v, v) = f∈Fx f (β•n f )[[L E (v)]]•{ {δ(v)} }, s β;x (v, ζv) -s ζβ;x (v, v) = f∈Fx f |β•n f |[[L E (v)]]•[[δ(v)]],
for all x ∈ F • or x ∈ C, and

A ∂ (β•n) -,f (v, ζv) -A ∂ (ζβ•n) -,f (v, v) = f (β•n) -L E (v)•δ(v),
for all f ∈ F ∂ . Then, applying the Cauchy-Schwarz inequality to these three terms yields

T 1,2 ≤ 6 |v| 2 ∂ + |v| 2 s 1 2   2 c∈C ||β|| L ∞ (c) p∈Pc ||δ(v)|| 2 L 2 (∂p)   1 2
.

In addition, observing that σ β,-∇β t + 1 2 (∇•β)Id ≡ 0 and using the identity (29), we have

A β,-∇β t + 1 2 (∇•β)Id (v, v) = 1 2 |v| 2 ∂ + |v| 2 s ,
so that combining this expression with the above estimate yields

T 1,2 ≤ 12 A β,-∇β t + 1 2 (∇•β)Id (v, v) 1 2   c∈C ||β|| L ∞ (c) p∈Pc ||δ(v)|| 2 L 2 (∂p)   1 2
.

Finally, we use the inequalities (27a)-(27b) together with the definition of ω ζ , to infer that

T 1,2 ≤ C δ ω ζ A β,-∇β t + 1 2 (∇•β)Id (v, v) 1 2 τ -1 |||v||| 2 2 1 2
where C δ > 0 depends exclusively on the numerical constants C T and C . Now, we collect the bounds on T 1,1 and T 1,2 and we apply Young's inequality to obtain

A β,-∇β t + 1 2 (∇•β)Id (v, w) ≥ 1 4 |||v||| 2 + (θ -C 2 δ ω 2 ζ )A β,-∇β t + 1 2 (∇•β)Id (v, v) + T 1,3 . As a result, choosing θ = C 2 δ ω 2 ζ yields A β,-∇β t + 1 2 (∇•β)Id (v, w) ≥ 1 4 |||v||| 2 + T 1,3 . (30)
(ii) Bound on T 2 . First, we rewrite this term as:

T 2 = θH µ+∇β t -1 2 (∇•β)Id (v, v) + H ζ(µ+∇β t -1 2 (∇•β)Id) (v, v) + H µ+∇β t -1 2 (∇•β)Id (v, ζv) -H ζ(µ+∇β t -1 2 (∇•β)Id) (v, v) = T 2,1 + T 2,2 + T 2,3 .
Concerning T 2,1 , we have

T 2,1 = θ 2 c∈C c σ β,µ L Ec (v)•L Ec (v) ≥ θλ 2 c∈C ||L Ec (v)|| 2 L 2 (c) ≥ C 2 λ 2 θ|||v||| 2 2 ,
where we have used hypothesis (H2) (recall that λ ≤ 0) and the upper bound from Proposition 4.1. The second term T 2,2 is treated similarly:

T 2,2 = 1 2 c∈C c ζσ β,µ L Ec (v)•L Ec (v) ≥ C 2 λ 2 ||ζ|| L ∞ (Ω) |||v||| 2 2 .
Collecting these bounds yields

T 2 ≥ λ 2 ϑ|||u||| 2 2 + T 2,3 , (31)
with ϑ = C 2 (θ + ||ζ|| L ∞ (Ω) ).
(iii) Bound on T 1 + T 2 . Collecting the estimates (30) and (31), we obtain

A β,µ (v, w) ≥ 1 4 |||v||| 2 + τ λ 2 ϑτ -1 |||v||| 2 2 + T 1,3 + T 2,3 . 20
We observe that

T 1,3 + T 2,3 = H 1 2 (∇•β)Id (v, ζv) -H 1 2 ζ(∇•β)Id (v, v) + H µ+∇β t -1 2 (∇•β)Id (v, ζv) -H ζ(µ+∇β t -1 2 (∇•β)Id) (v, v) = H µ+∇β t (v, ζv) -H ζ(µ+∇β t ) (v, v) = c∈C c µ + ∇β t L Ec (v)•δ(v).
Applying successively the Hölder inequality, the inequality (27a) and the upper bound from Proposition 4.1, we infer that

|T 1,3 +T 2,3 | ≤ c∈C ||µ+∇β t || L ∞ (c) ||L E (v)|| L 2 (c) ||δ(v)|| L 2 (c) ≤ ||µ+∇β t || L ∞ (Ω) C 2 L ζ h|||v||| 2 2 .
As a result, we obtain

A β,µ (v, w) ≥ 1 4 |||v||| 2 + τ λ 2 ϑ - h 4h 0 τ -1 |||v||| 2 2 ,
with the reference length

h 0 = 4C 2 ||µ + ∇β t || L ∞ (Ω) τ L ζ -1
. Hence, there exists > 0 such that A β,µ (v, w) ≥ |||v||| 2 , as soon as λ and h satisfy (25). Indeed, the first term T 1 is obtained using the definition (11) of g β,µ;c together with the following integration by part formula (28) and Let q ∈ [1, 2] and denote q ≥ 2 its conjugate number, i.e. , 1 = 1/q + 1/q . From the Hölder inequality, we infer that L q (c)

Bound on consistency error and a priori estimate

1 q c∈C ||L Ec (v)|| q L 2 (c) 1 q
.

Moreover, recalling that q ≥ 2 so that |•| q ≤ |•| 2 , and using the upper bound in Proposition 4.1 leads to

|T 1 | = c∈C c (∇β + µ t -(∇•β)Id) L Ec (v)•y c∈C N q ∞ h θq c ||y|| q L q (c) 1 q |||v||| 2 .
To bound the two terms T 2 and T 3 , we consider a sub-face f ∈ F x for all x ∈ X with X ∈ {F • , C}. As above, the Hölder inequality yields

f (β•n f )[[L E (v)]]•{ {y} } ≤ ||β|| 1 2 L ∞ (f) ||{ {y} }|| L q (f) |||β•n f | 1 2 [[L E (v)]]|| L q (f) ,
so that using a local inverse inequality, we obtain

f (β•n f )[[L E (v)]]•{ {y} } ≤ h θ f ||β|| 1 2 L ∞ (f) ||{ {y} }|| L q (f) |||(β•n f )| 1 2 [[L E (v)]]|| L 2 (f) ,
with θ = (d -1) 1 q -1 2 . Hence, denoting f = X∈{F • ,C} x∈X f∈Fx , it follows from the triangle inequality, the Hölder inequality and q ≥ 2 that f f

(β•n f )[[L E (v)]]•{ {y} } ≤   f h θ q f ||β|| q 2 L ∞ (f) ||{ {y} }|| q L q (f)   1 q   f |||β•n f | 1 2 [[L E (v)]]|| 2 L 2 (f)   1 2
.

Next, owing to the definitions ( 12) and ( 14) of n β and s β respectively, the mesh regularity and recalling the inequality |a ± b| q ≤ 2 q-1 (|a| q + |b| q ), we infer that

|T 2 + T 3 |   c∈C p∈Pc h θ q c ||β|| q 2 L ∞ (c) ||y|| q L q (∂p)   1 q s β (v, v) 1 2
.

Finally, proceeding similarly, we also infer that

|T 4 |   f ∈F ∂ h θ q c f ||β|| q 2 L ∞ (f ) ||y|| q L q (f )   1 q A ∂ |β•n| (v, v) 1 2
, and the expected result follows from the above bounds.

x c 2 Figure 1 :

 21 Figure 1: Left panel: tetrahedron [x v , x v , x f , x c ]. Right panel: local diamond p e,c .

Figure 2 :

 2 Figure 2: In blue. Left: intra-cell sub-face f = ∂p e,c ∩ ∂p e ,c ∈ F c . Right: inter-cell sub-face f = ∂p e,c ∩ ∂p e,c ∈ F f .

  e,c | p e ,c = Idfc (e )⊗e d|p e ,c | fc (e) |c| + fc (e) d|p e,c | δ e,e , ∀e ∈ E c ,

Figure 3 :

 3 Figure 3: Local dual face fc (e).

4. 1

 1 Properties of the reconstruction map

Proposition 4 . 1 (

 41 Stability). There exists C > 0 such that

Remark 4 . 3 (

 43 De Rham's map). Requiring more regularity, the usual de Rham's reduction map defined by R E (v)| e = |e| -1

  induced by the bilinear form A |β•n| defined by (17), and |•| 2 s := s β (•, •) is the semi-norm induced by the bilinear form s β defined by (14).

Proposition 4 . 6 (

 46 Coercivity). Assume that (H1) holds. Then,

Proposition 4 . 7 (

 47 Inf-sup stability). Assume that (H2) and (25) hold. Then, there exists > 0 such that |||v||| ≤ sup w∈E,|||w|||=1

Figure 4 :

 4 Figure4: Examples of meshes from the four sequences. From left to right: hexahedral mesh (H), prismatic mesh with triangular basis (PrT), prismatic mesh with hexagonal basis (PrG), and Checkerboard mesh with non-matching interfaces (CB).

1 Figure 5 :

 15 Figure 5: Left panel: inflow boundary ∂Ω -in blue and some streamlines of the vector field β. Right panel: Discrete errors on H (

6 Proofs 6 . 1

 661 |p e,c | δ e,e -e • fc (e) |p e ,c | fc (e ) |c| Stabilization term , ∀e ∈ E c . Numerical experiments show that it is possible to replace the parameter d -1 by a positive value that is reasonably close to d -1 ; however, in the stability analysis, this modification impacts the property ( 1) which is used to obtain the lower bound in Proposition 4.1. Properties of the reconstruction map Proof of Proposition 4.1. Let c ∈ C, let v ∈ E c and let q ∈ [1, ∞). (i) Lower bound. Owing to the definition (5) of L Ec , we have for all e ∈ E c , L Ec (v) |pe,c = v e a e + b e with a e = e |e| 2 and b e = e,c -e |e| 2 v e + e ∈Ec\{e} v e e ,c .

2 , 2 4 .

 224 e,c | p e,c 2 ≤ | fc (e)| |c| |c| d|p e,c | and e,c | p e ,c 2 ≤ | fc (e)| |c| 1 + 1 cos 2 (t e , n fc(e ) ) 1 leading to |e||| e,c || L ∞ (c) ≤ |e|| fc (e)| |c| max (t e , n fc(e ) ) Let c ∈ C and let e ∈ E c . The consistency property relies on the property ( 2). Indeed, given U ∈ P 0 (ĉ; R d ), we infer that, for all x ∈ p e ,c , L Ec R Ec (U )(x) = e∈Ec R Ec (U ) |e e,c (x) = e∈Ec (U •e) e,c (x) = e∈Ec e,c (x)⊗e U = U .

Lemma 6 . 1 (

 61 Multiplicative stability). There exists C ζ > 0 independent of the mesh size and the model parameters such that |||ζv||| ≤ C ζ ||ζ|| L ∞ (Ω) + ω ζ |||v|||, ∀v ∈ E. Proof. Let v ∈ E and let us rewrite |||ζv||| 2 as

Proof of Lemma 4 . 9 .

 49 Let y |c = (u -I E (u)) |c for all c ∈ C. Note that y |∂p ∈ L q (∂p) for all p ∈ P c . Let v ∈ E. Owing to the definitions of A β,µ and Σ, it follows that Σ(s,u D ; v) -A β,µ (R E (u), v) = T 1 + T 2 + T 3 + T 4 , with T 1 := c∈C c (∇β + µ t -(∇•β)Id) L Ec (v)•y, T 2 := X∈{F • ,C} x∈X f∈Fx f (β•n)[[L E (v)]]•{ {y} }, T 3 := X∈{F • ,C} x∈X f∈Fx f |β•n| [[L E (v)]]•[[y]] and T 4 := f ∈F ∂ f∈F ∂ f f (β•n) -L Ec f (v)•y.

  •∇)y) •L Ec (v) = -p∈Pc p ((β•∇)L Ec (v)) •yc (∇•β) L Ec (v)•y+ p∈Pc p ∇• (β y•L Ec (v)) ,holding for all c ∈ C and all v ∈ E c . The terms T 2 and T 3 result from the rightmost term of the relation (28) and the fact that (β•n)[[u]] |f ≡ 0 for all f ∈ F x . Finally, the term T 4 is inferred observing that u D = u |∂Ω . It remains to bound these four terms. First, let us consider T 1 .

c

  (∇β + µ t -(∇•β)Id) L Ec (v)•y ≤ N ∞ ||y|| L q (c) ||L Ec (v)|| L q (c) , with N ∞ = ||∇β + µ t -∇•βId|| L ∞ (c) .Then, using a local inverse inequality (see[START_REF] Ern | Theory and practice of finite elements[END_REF] Lemma 1.138]), we infer thatc (∇β + µ t -(∇•β)Id) L Ec (v)•y ≤ N ∞ h θ c ||y|| L q (c) ||L Ec (v)|| L 2 (c) ,with θ = d 1 2 -1 q , so that the Hölder inequality yields c∈C c(∇β + µ t -(∇•β)Id) L Ec (v)•y ≤ c∈C N q ∞ h θq c ||y|| q

Table 5 :

 5 Features of prism meshes

		PrT20 PrT30 29 791 114 390 138 600 9 261 34 860 41 600	16 000 54 000	(a) PrT10 Mesh
		PrT40 68 921 267 320 326 400 128 000
	mes for Stokes Problem: Principles and First Results le Discrete Operator Schemes for Stokes Problem: Principles and First Results H-I83-2013-03326-EN Version 1.0	H-I83-2013-03326-EN Version 1.0
							(b) PrT20 Mesh
	# F # V	# C # E M	# F # V	# C # E	# F	# C
	240 1 728 13 056 1 331 9 261 29 791 114 390 138 600 64 512 4 730 5 400 34 860 41 600 PrG10 3 080 PrG20 20 160 4 096 PrG30 63 240 154 200 119 791 28 830 (a) H4 Mesh 2 000 16 000 7 200 5 331 1 210 (a) PrT10 Mesh 48 600 37 261 8 820 54 000 101 376 32 768 68 921 267 320 326 400 128 000 PrG40 144 320 354 000 276 921 67 240
	rtesian meshes ble 5: Features of prism meshes
					(b) H8 Mesh		(b) PrT20 Mesh
	# F # V	# C # E	# F	# C
	120 864 6 528 50 688 399 360 196 608 48 384 144 320 354 000 276 921 67 240 24 576 63 240 154 200 119 791 28 830 3 072 (a) TU3 Mesh 3 080 7 200 5 331 1 210 20 160 48 600 37 261 8 820		(a) PrG10 Mesh
	tetrahedral meshes ures of prism meshes with polygonal basis
					(b) TU4 Mesh		(b) PrG20 Mesh
	#F	#C
	500 # V	215	# E	# F	# C
	4 308 8 248 6 148 edral meshes 2 003 3 898 7: Features of checkerboard meshes 4 988 61 052 254 977 700 416 592 896 147 456 2 787 30 480 33 025 89 088 74 496 18 432 1 691 15 266 4 417 11 520 9 408 2 304 7 711 (c) T2 Mesh 97 216 156 36 625 1 536 1 200 288		(a) CB4 Mesh
					(d) T3 Mesh		(b) CB8 Mesh
	Page 29 sur 54			©EDF SA 2013
	nte			Page 30 sur 54		©EDF SA 2013

Table 6 :

 6 Features of prism meshes with polygonal basis

	(a) PrG10 Mesh
	(b) PrG20 Mesh

Table 7 :

 7 Features of checkerboard meshes

	(a) CB4 Mesh
	(b) CB8 Mesh

Table 6 :

 6 Features of prism meshes with polygonal basis

	(b) PrG20 Mesh

Table 7 :

 7 Features of checkerboard meshes

	(a) CB4 Mesh
	(b) CB8 Mesh

  Mean stencil St and discrete error Er E (u) for the H (upper left panel), PrT (upper right panel), PrG (lower left panel), and the CB (lower right panel) mesh sequences.

					#E	St Er E (u) #E #V	#E #C
			3.9e-01 2.40 4.69 4.7e+03 38 2.4e-01 3.55 2.37
	1.9e+03 25 1.8e-01 2.67 3.80 3.5e+04 46 1.5e-01 3.76 2.18
	1.4e+04 28 9.4e-02 2.82 3.39 1.1e+05 48 1.1e-01 3.84 2.12
	1.0e+05 30 4.9e-02 2.91 3.19 2.7e+05 49 9.1e-02 3.88 2.09
	#E	St Er E (u) #E #V	#E #C	#E	St Er E (u) #E #V	#E #C
	7.2e+03 83	2.2e-01 2.34 5.95 1.5e+03 112 3.6e-01 2.46 5.33
	4.9e+04 110 1.4e-01 2.41 5.51 1.2e+04 144 1.8e-01 2.61 5.00
	1.5e+05 120 1.1e-01 2.44 5.35 8.9e+04 162 9.8e-02 2.70 4.83
	3.5e+05 125 8.5e-02 2.45 5.26 7.0e+05 180 5.1e-02 2.75 4.75
	Figure 6:

Remark 5.1 (Stabilization parameter). As observed in Bonelle et al.

[START_REF] Bonelle | Low-order reconstruction operators on polyhedral meshes: application to compatible discrete operator schemes[END_REF]

, one can reformulate the basis functions e,c as a consistent term plus a stabilization term:


Erik Burman

Alexandre Ern

A nonlinear consistent penalty method weakly enforcing positivity in the finite element approximation of the transport equation

Keywords: Stabilized finite element method, consistent penalty, positivity preserving, transport equation, discrete maximum principle

We devise and analyze a new stabilized finite element method to solve the first-order transport (or advection-reaction) equation. The method combines the usual Galerkin/Least-Squares approach to achieve stability with a nonlinear consistent penalty term inspired by recent discretizations of contact problems to weakly enforce a positivity condition on the discrete solution. We prove the existence and the uniqueness of the discrete solution. Then we establish quasi-optimal error estimates for smooth solutions bounding the usual error terms in the Galerkin/Least-Squares error analysis together with the violation of the maximum principle by the discrete solution. Numerical examples are presented to illustrate the performances of the method.

Introduction

The design of robust and accurate finite element methods for first-order transport (or advection-reaction) equations or for advection-dominated advectiondiffusion equations remains an active field of research. Indeed, the task of designing a numerical scheme that is of higher order than one in the zone where the exact solution is smooth, but preserves the monotonicity properties of the exact solution on the discrete level, is nontrivial. Since it is known that such a scheme necessarily must be nonlinear even for linear equations, one typical strategy adopted when working with stabilized finite element methods is to add an additional nonlinear shock-capturing term, designed to make the method satisfy a discrete maximum principle; see, e.g., [START_REF] Mizukami | A Petrov-Galerkin finite element method for convection-dominated flows: an accurate upwinding technique for satisfying the maximum principle[END_REF][START_REF] Burman | Stabilized Galerkin approximation of convectiondiffusion-reaction equations: discrete maximum principle and convergence[END_REF]. These methods, however, often result in ill-conditioned nonlinear equations and include parameters that may be difficult to tune and that depend on the mesh geometry. Another approach is the so-called flux-corrected finite element method [START_REF] Löhner | Finite element flux-corrected transport (FEM-FCT) for the Euler and Navier-Stokes equations[END_REF][START_REF] Kuzmin | Flux correction tools for finite elements[END_REF]. In this scheme, the system matrix is manipulated so that it becomes a so called M-matrix, the inverse of which has positive coefficients which yields a maximum principle. Such a scheme is monotonicity preserving, but of first order. In order to improve the accuracy, anti-diffusive mechanisms, or flux-limiter techniques, have been proposed so as to reduce the amount of dissipation in the smooth regions by blending a low-and a high-order approximation [START_REF] Kuzmin | Flux correction tools for finite elements[END_REF][START_REF] Guermond | A second-order maximum principle preserving Lagrange finite element technique for nonlinear scalar conservation equations[END_REF].

In this paper, we consider a method that follows a completely different approach. The starting observation is that, if the problem satisfies a maximum principle of the form u ≥ 0, [START_REF] Mizukami | A Petrov-Galerkin finite element method for convection-dominated flows: an accurate upwinding technique for satisfying the maximum principle[END_REF] then this constraint can be added to the problem without any perturbation. On the discrete level however, the condition (1) is not necessarily satisfied, unless enforced by the numerical method, which is the purpose of all the methods discussed above. One may argue that one can solve the problem with any method under the constraint (1). This was proposed in [START_REF] Evans | Enforcement of constraints and maximum principles in the variational multiscale method[END_REF]. The resulting method, in the form of a variational inequality, is unwieldy, with the need of Lagrange multipliers to impose the constraint and associated stability and solver issues. In the present work, we instead draw on recent advances in the field of contact problems [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Burman | Galerkin least squares finite element method for the obstacle problem[END_REF], where the variational inequality instead is discretized by means of a nonlinear consistent penalty method. Note however that in the present context the formulation cannot be associated with an augmented Lagrangian method, due to the lack of symmetry of the formulation.

Our method combines the well-known Galerkin/Least Squares (GaLS) discretization of the transport equation (see [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF][START_REF] Ern | Theory and Practice of Finite Elements[END_REF]) with a nonlinear switch inspired from [START_REF] Burman | Galerkin least squares finite element method for the obstacle problem[END_REF] and that changes the equations in the zones where ( 1) is violated to a least-squares penalty on this inequality (more precisely, the negative part of the discrete solution) together with a least-squares penalty on the residual. Our method is not meant to enforce strictly a discrete maximum principle, but to blend asymptotically the satisfaction of the GaLS approximation of the PDE with the satisfaction of the discrete maximum principle, in the same spirit of the above-mentioned methods for contact problems. We first prove the existence and uniqueness of the discrete solution. Then, our main result is Theorem 3.1 where we establish a quasi-optimal error estimate bounding at the same time the error measured in the usual GaLS norm (combining the L 2 -norm on the solution, its boundary values, and the advective derivative weighted by the local mesh size) and the violation of the positivity condition (1) measured by the weighted L 2 -norm of the negative part of the discrete solution. These convergence results hold for all polynomial orders k ≥ 1. Another salient feature of our method is its flexibility in incorporating a priori lower and upper bounds on the discrete solution by simply adding the corresponding consistent penalty term to the discrete formulation. Finally, we report some numerical experiments illustrating that accurate solutions with mild and asymtotically vanishing violations of the discrete maximum principle can be obtained at moderate computational costs.

Model problem and GaLS discretization

Let Ω be an open, bounded, Lipschitz set in R d , d ∈ {2, 3}, let β ∈ W 1,∞ (Ω; R d ) be a given advection velocity, and let σ ∈ L ∞ (Ω; R) be a given reaction coefficient. We assume that β and σ satisfy the following (classical) positivity assumption: There exists σ 0 > 0 such that

0 < σ 0 ≤ σ - 1 2 ∇•β, a.e. in Ω. (2)
We split the boundary ∂Ω of Ω as

∂Ω = ∂Ω -∪ ∂Ω 0 ∪ ∂Ω + with ∂Ω -= {x ∈ ∂Ω | (β•n)(x) < 0} (inflow boundary), ∂Ω 0 = {x ∈ ∂Ω | (β•n)(x) = 0} (
u : Ω → R such that A(u) := β•∇u + σu = f in Ω, (3a)
u = g on ∂Ω -. (3b)
We assume that f ∈ L 2 (Ω; R) and g ∈ L 2 |β•n| (∂Ω -; R) and look for a weak solution in the graph space

V := {v ∈ L 2 (Ω; R) | β•∇v ∈ L 2 (Ω; R)}.
Assuming that dist(∂Ω -, ∂Ω + ) > 0, one can show [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] that functions in the graph space admit a trace in the weighted space L 2 |β•n| (∂Ω; R) and that there exists one and only one weak solution in the graph space V to the model problem [START_REF] Löhner | Finite element flux-corrected transport (FEM-FCT) for the Euler and Navier-Stokes equations[END_REF]. In particular, we observe that the following positivity condition holds:

σ 0 v 2 Ω + 1 2 |β•n| 1 2 v 2 ∂Ω ≤ (β•∇v + σv, v) Ω -(β•n)v, v ∂Ω -, ∀v ∈ V. (4)
Several stabilized H 1 -conforming finite element methods are available in the literature to discretize [START_REF] Löhner | Finite element flux-corrected transport (FEM-FCT) for the Euler and Navier-Stokes equations[END_REF]. We focus on the Galerkin/Least-Squares method (GaLS). Let T h be a mesh from a shape-regular mesh sequence. We assume for simplicity that Ω is a polytope (polygon/polyhedron) so that T h can cover Ω exactly. Let k ≥ 1 be the polynomial degree and consider the H 1 -conforming finite element space

V k h := {v h ∈ C 0 (Ω; R) | v h | T ∈ P k (T ; R), ∀T ∈ T h }, (5)
where P k (T ; R) denotes the space composed of R-valued functions that are the restriction to T of d-variate polynomials of degree at most k. We consider the following discrete problem: Find

u h ∈ V k h such that a τ h (u h , w h ) = ℓ τ h (w h ), ∀w h ∈ V k h , (6)
with the following bilinear and linear forms:

a τ h (v h , w h ) := (A(v h ), w h + τ A(w h )) Ω -(β•n)v h , w h ∂Ω -, (7a) ℓ τ h (w h ) := (f, w h + τ A(w h )) Ω -(β•n)g, w h ∂Ω -. (7b)
The stabilization parameter τ is piecewise constant on T h and is of the order of the local mesh size h T for all T ∈ T h ; more precisely, a suitable choice is

τ | T = min(σ -1 0 , β -1 T h T ) with β T = β L ∞ (T ;R d ) .
By construction, the discrete bilinear form a τ h is coercive with respect to the norm:

a τ h (v, v) ≥ |||v||| 2 := σ 0 v 2 Ω + |β•n| 1 2 v 2 ∂Ω + τ 1 2 A(v) 2 Ω , ∀v ∈ V. (8)
Moreover, exact consistency holds, and the following quasi-optimal error estimates can be established [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF][START_REF] Ern | Theory and Practice of Finite Elements[END_REF]: There exists C, uniform, such that 9)

|||u-u h ||| ≤ C inf v h ∈V k h ( τ -1 2 (u-v h ) Ω + |β•n| 1 2 (u-v h ) ∂Ω + τ 1 2 A(u-v h ) Ω ), (
and if u ∈ H k+1 (Ω), |||u -u h ||| ≤ C( T ∈T h φ T h 2k+1 T |u| 2 H k+1 (T ;R) ) 1 2 with φ T = max(β T , σ 0 h T ).

The consistent penalty method

The model problem (3) has a maximum principle; for instance, if f ≥ 0 and g ≥ 0, then u ≥ 0 in Ω. Unfortunately, this property rarely carries over to finite element discretizations. Our goal is to modify the GaLS finite element approximation (6) by using a consistent penalty method.

Let γ > 0 be a penalty parameter. For any function v ∈ V , let us define the function

ξ γ : Ω → R such that ξ γ (v) := [v -γ(A(v) -f )] -, (10)
where x -= 1 2 (x -|x|) denotes the negative part of the real number x. Note that ξ γ (u) = 0 in Ω for the weak solution u since A(u) = f and u -= 0. Let us consider the following discrete problem: Find

u h ∈ V k h such that a τ γ h (u h ; w h ) = ℓ τ h (w h ), ∀w h ∈ V k h , (11)
with

a τ γ h (v h ; w h ) := a τ h (v h , w h ) + (γ -1 ξ γ (v h ), w h ) Ω . (12)
Since ξ γ (u) vanishes identically in Ω, exact consistency still holds for [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF]. The discrete problem [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] remains meaningful and exactly consistent if the penalty parameter γ is replaced by a function taking uniformly positive values in Ω.

We will take γ to be piecewise constant on the mesh T h since the error analysis below will reveal that quasi-optimal error estimates are obtained by taking γ to be locally of the order of h T (on quasi-uniform mesh sequences, a constant function γ can be considered).

Rationale of the consistent penalty method

Before embarking on the analysis of the method, let us briefly discuss the design principle behind the approach. First, we observe that if [u h -γ(A(u h )-f )] -= 0, then the formulation coincides with the standard GaLS discretization. Assume now that [u hγ(A(u h )f )] -= 0 everywhere in the macroelement Ω i := supp(ϕ i ) where ϕ i is an interior nodal (or hat) basis function. Then, since x -= x if x -= 0, we observe that the standard Galerkin part is eliminated by the second term in the penalty term, so that [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] with

w h = ϕ i becomes (γ -1 u h , ϕ i ) Ωi + (τ (A(u h ) -f ), A(ϕ i )) Ωi = 0. ( 13
)
This shows that the nonlinear penalty term changes the discrete equation locally to the sum of two least-squares contributions, one on the violation of positivity by the discrete solution and one on the PDE residual. By choosing γ small, one can expect that the violation of the maximum principle is reduced. This is indeed one of the main conclusions of the error analysis below, where we additionally prove that quasi-optimal error estimates of the form (9) also hold for the consistent penalty method.

Remark 3.1 (Variant) Variants are possible in [START_REF] Temam | Navier-Stokes equations[END_REF] for the penalty term, such as

a τ γ h (v h ; w h ) := a τ h (v h , w h ) + (γ -1 ξ γ (v h ), w h + τ A(w h )) Ω . ( 14
)
This variant can be analyzed using the same arguments as below. In particular, considering an interior nodal basis function, we now obtain (γ

-1 u h , ϕ i + τ A(ϕ i )) Ωi = 0.
Comparing with [START_REF] Kuzmin | Goal-oriented mesh adaptation for flux-limited approximations to steady hyperbolic problems[END_REF], only the penalty on the violation of the positivity remains, but the term is no longer symmetric.

Well-posedness and convergence

Let us first establish that a τ γ h has reasonable monotonicity properties. Lemma 3.1 (Monotonicity) Assume that 0 < γ ≤ τ . Then the following holds for all u 1 , u 2 ∈ V :

1 2 (|||u 1 -u 2 ||| 2 + γ -1 2 (ξ γ (u 1 ) -ξ γ (u 2 )) 2 Ω ) ≤ a τ γ h (u 1 ; u 1 -u 2 ) -a τ γ h (u 2 ; u 1 -u 2 ), (15a) 1 4
(|||u

1 ||| 2 + γ -1 2 ξ γ (u 1 ) 2 Ω ) ≤ a τ γ h (u 1 ; u 1 ) + γ 1 2 f 2 Ω . (15b)
Proof. Let us prove (15a). We observe that

a τ γ h (u 1 ;u 1 -u 2 ) -a τ γ h (u 2 ; u 1 -u 2 ) = a τ h (u 1 -u 2 , u 1 -u 2 ) + (γ -1 (ξ γ (u 1 ) -ξ γ (u 2 )), u 1 -u 2 ) Ω ≥ |||u 1 -u 2 ||| 2 + (γ -1 (ξ γ (u 1 ) -ξ γ (u 2 )), u 1 -u 2 ) Ω ,
where we have used [START_REF] Burman | Galerkin least squares finite element method for the obstacle problem[END_REF]. Moreover, we have

(γ -1 (ξ γ (u 1 ) -ξ γ (u 2 )), u 1 -u 2 ) Ω = (γ -1 (ξ γ (u 1 ) -ξ γ (u 2 )), u 1 -γ(A(u 1 ) -f ) -(u 2 -γ(A(u 2 ) -f ))) Ω + (ξ γ (u 1 ) -ξ γ (u 2 ), A(u 1 -u 2 )) Ω ≥ γ -1 2 (ξ γ (u 1 ) -ξ γ (u 2 )) 2 Ω + (ξ γ (u 1 ) -ξ γ (u 2 ), A(u 1 -u 2 )
) Ω , where we have used the fact that

|x --y -| 2 ≤ (x --y -)(x -y), ∀x, y ∈ R. ( 16
)
Using Young's inequality and the fact that γ ≤ τ , we infer that

((ξ γ (u 1 )-ξ γ (u 2 )), A(u 1 -u 2 )) Ω ≥ - 1 2 γ -1 2 (ξ γ (u 1 )-ξ γ (u 2 )) 2 Ω - 1 2 τ 1 2 A(u 1 -u 2 ) 2 Ω .
Putting everything together shows that (15a) holds. Finally, the proof of (15b) follows from (15a) by taking u 2 = 0 and using the fact that

1 2 γ -1 2 (ξ γ (u 1 ) - ξ γ (0)) 2 Ω ≤ γ -1 2 ξ γ (u 1 ) 2 Ω + γ -1 2 ξ γ (0) 2 Ω ≤ γ -1 2 ξ γ (u 1 ) 2 Ω + γ 1 2 f 2 Ω .
We can now prove that the discrete problem [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] is well-posed. Proposition 3.1 (Well-posedness) Assume that 0 < γ ≤ τ . Then the discrete problem [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] admits one and only one solution.

Proof. Uniqueness follows from (15a). To prove existence, let N := dimV k h and let G : R N → R N be the map defined by (G(U ), V ) R N := a τ γ h (u h ; v h )ℓ τ h (v h ), where U, V ∈ R N are the component vectors associated with the functions u h , v h in the Lagrange basis of V k h . It is readily seen that G is a continuous map (observe in particular that |x --y -| ≤ |x-y| for all x, y ∈ R). Furthermore, since Cauchy-Schwarz inequalities and τ ≤ σ -1

0 show that |ℓ τ h (v h )| ≤ K|||v h ||| with K = (2σ -1 2 0 f Ω + |β•n| 1 2 g ∂Ω -), we infer using (15b) that (G(U ), U ) R N = a τ γ h (u h ; u h ) -ℓ τ h (u h ) ≥ 1 4 (|||u h ||| 2 + γ -1 2 ξ γ (u h ) 2 Ω ) -γ 1 2 f 2 Ω -K|||u h |||.
This proves that there is a real number, say K ′ , so that (G(U ), U ) R N > 0 for all U ∈ R N with U R N ≥ K ′ . Indeed, using norm equivalence on discrete spaces, we infer that there exists

C N > 0 such that C N U R N ≤ |||u h ||| for all U ∈ R N with associated discrete function u h ∈ V k h . This leads to (G(U ), U ) R N ≥ 1 8 |||u h ||| 2 -γ 1 2 f 2 Ω -2K 2 ≥ 1 8 C 2 N U 2 R N -γ 1 2 f 2 Ω -2K 2 ,
and we conclude that the expected inequality holds with

K ′ = 8 C N γ 1 2 f 2 Ω + 2K 2 + 1.
Existence then follows using well-known arguments (see, for instance, [12, Lemma 1.4, Chapter 2]). The next theorem is the main result of this paper. It shows that the GaLS finite element method with penalty has essentially the same behavior as that without penalty when approximating smooth solutions. Theorem 3.1 (Error estimate) Let u ∈ V be the solution to (3) and let u h ∈ V k h be the solution to [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF]. Assume that 0 < γ ≤ τ . Then there exists C > 0, uniform, such that

|||u -u h ||| + γ -1 2 [u h ] -Ω ≤ C inf v h ∈V k h ( τ 1 2 A(u -v h ) Ω + |β•n| 1 2 (u -v h ) ∂Ω + γ -1 2 (u -v h ) Ω ). ( 17
)
Moreover, if u ∈ H k+1 (Ω), τ is chosen as in the GaLS method as τ | T = min(σ -1 0 , β -1 T h T ), and cτ | T ≤ γ| T for all T ∈ T h with c uniformly bounded from below away from zero, then

|||u -u h ||| + γ -1 2 [u h ] -Ω ≤ C T ∈T h φ T h 2k+1 T u 2 H k+1 (T ;R) 1 2 . ( 18
)
Proof. Let e = uu h . Then, using (8), we infer that

|||e||| 2 ≤ a τ h (e, e) = a τ h (e, u -v h ) + a τ h (e, v h -u h ).
Moreover, the exact consistency of the GaLS approximation and the definition of the discrete problem [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] imply that

a τ h (e, v h -u h ) = ℓ τ h (v h -u h ) -ℓ τ h (v h -u h ) + (γ -1 ξ γ (u h ), v h -u h ) Ω = (γ -1 ξ γ (u h ), v h -u h ) Ω = (γ -1 ξ γ (u h ), v h -u + γA(e)) Ω + (γ -1 ξ γ (u h ), u -u h -γA(e)) Ω .
Since ξ γ (u) = 0, using the monotonicity property (16), we infer that

(γ -1 ξ γ (u h ), u -u h -γA(e)) Ω ≤ -γ -1 2 ξ γ (u h ) 2 Ω .
As a result, we obtain

|||e||| 2 + γ -1 2 ξ γ (u h ) 2 Ω ≤ a τ h (e, u -v h ) + (γ -1 ξ γ (u h ), v h -u + γA(e)) Ω .
The boundedness properties of the GaLS approximation yield

a τ h (e, u -v h ) ≤ |||e|||( τ 1 2 A(u -v h ) Ω + |β•n| 1 2 (u -v h ) ∂Ω + τ -1 2 (u -v h ) Ω ).
Moreover, we have

(γ -1 ξ γ (u h ), v h -u + γA(e)) Ω ≤ γ -1 2 ξ γ (u h ) Ω ( γ -1 2 (u -v h ) Ω + γ 1 2 A(e) Ω ) ≤ 3 4 γ -1 2 ξ γ (u h ) 2 Ω + γ -1 2 (u -v h ) 2 Ω + 1 2 γ 1 2 A(e) 2 Ω .
Collecting these two bounds and using that γ ≤ τ , we infer that |||e|||+ γ -1 2 ξ γ (u h ) Ω is bounded by the right-hand side of (17). To conclude that (17) holds, it suffices to establish that

γ -1 2 [u h ] -Ω ≤ γ -1 2 ξ γ (u h ) Ω + |||e|||.
This inequality, in turn, follows from the elementary inequality |[x + y] -| ≤ |x -| + |y -| for arbitrary real numbers x and y, leading to

γ -1 2 [u h ] -Ω ≤ γ -1 2 ξ γ (u h ) Ω + γ 1 2 [(A(u h ) -f )] -Ω ,
together with the fact that A(u) = f and the assumption that γ ≤ τ so that γ

1 2 [(A(u h ) -f )] -Ω ≤ |||e|||.
Finally, (18) results from the approximation properties of finite elements and the assumptions on τ and γ.

The error estimates derived in Theorem 3.1 show that the present consistent penalty method delivers similar bounds on the error |||uu h ||| to those obtained with the usual GaLS discretization, while additionally controlling the violation of positivity by means of the measure γ -1 2 [u h ] -Ω (note that γ scales as the mesh size, so that the factor γ -1 2 in front of [u h ] -makes the bound even stronger.

Numerical example

In this section, we assess the proposed method on two test cases: the first one features a solution with inner layer, and the second one a solution with discontinuity.

Test case 1: solution with inner layer

We consider problem (3) in the domain Ω ⊂ R 2 shown in the left panel of Figure 1, with

β = 1 (x 2 + y 2 ) 1 2
(y, -x) T , σ = 0, f = 0.

The advection field β rotates clockwise, and the inflow boundary corresponds to the part of ∂Ω where x = 0. The exact solution given by u = 1 2 (tanh(((x 2 + y 2 ) 1 2 -0.5)/ǫ) + 1.0) is a consequence of inflow data imposed on the inflow boundary (see the contourlines in the right panel of Figure 1). The boundary data has a sharp layer of width ǫ at y = 0.5. This creates spurious underand overshoots that are transported downstream throughout the domain. We compute approximate solutions for both a mild layer (ǫ = 0.1) and a sharp layer (ǫ = 0.01) using either the nonlinear method [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF] with the bilinear form a τ γ h defined by [START_REF] Temam | Navier-Stokes equations[END_REF] or the standard (linear) GaLS method obtained by dropping the nonlinear term. Computations not reported here indicate that using the definition (14) leads to similar results. We consider affine (k = 1) and quadratic (k = 2) finite elements on a sequence of quasi-uniform unstructured meshes characterized by the mesh sizes h = 0.09 × 2 -l with l ∈ {0, 1, 2, 3, 4}. The nonlinear penalty term is evaluated using nodal quadrature when k = 1 and a quadrature that is exact for polynomials of order up to five when k = 2. Since the integrand of the nonlinear term is only H 1 in the interior of the elements, an exact quadrature requires a careful local analysis of where the nonlinearity is active. To assess the effect of possible under-integration of the nonlinear term, we also consider, in the case where k = 2, evaluating the nonlinear term with a hybrid quadrature rule obtained as a linear combination of two low-order terms, one using nodal quadrature and the other midpoint quadrature. The hybrid quadrature is insufficient to resolve the integration of the nonlinearity, but it gives control on the degrees of freedom of the quadratic polynomials.

We set the penalty parameter to γ = 0.0001h for k = 1 and to γ = 0.005h for k = 2 when consistent quadrature is used and to γ = 0.0001h for the loworder quadrature. Our results show that these choices are sufficient to reduce undershoots to less than one percent in all cases. Strengthening the penalty l e Ω β•∇e Ω e min e max Φ(u h ) 0 1.13e-2 (-) 9.13e-2 (-) -1.09e-2 4.59-3 1 2.33e-3 (2.3) 2.65e-2 (1.8) -1.93e-3 1.49e-4 2 9.57e-4 (1.3) 1.94e-2 (0.4) -9.20e-4 2.72e-5 3 2.54e-4 (1.9) 9.92e-3 (1.0) -3.90e-4 8.59e-6 4 6.35e-5 (2.0) 4.64e-3 (1.1) -1.6e-4 2.17e-6

Table 1: Test case 1, mild layer (ǫ = 0.1). Nonlinear method with consistent penalty,

k = 1, γ = 0.0001h, τ = h/2 l e Ω β•∇e Ω e min e max Φ(u h ) 0 1.18e-2 (-) 9.69e-2 (-) 2.37e-2 1.09e-2 3.44e-4 1 2.33e-3 (2.3) 2.65e-2 (1.9) - 1.93e-3 1.45e-4 2 9.57e-4 (1.3) 1.94e-2 (0.4) - 9.
20e-4 2.72e-5 3 2.54e-4 (1.9) 9.92e-3 (1.0) -3.90e-4 8.59e-6 4 6.35e-5 (2.0) 4.64e-3 (1.1) -1.60e-4 2.17e-6

Table 2: Test case 1, mild layer (ǫ = 0.1). Linear GaLS method, k = 1, τ = h/2 in the case of quadratic approximation does not improve the positivity, but increases the stiffness of the nonlinear problem. We present in Tables 1 to 5 the results for (i) the error e = uu h in the L 2 -norm and in the streamline derivative (experimental convergence orders are given between parenthesis), (ii) the violations of the maximum principle evaluated as e min :=min xi∈N (u h (x i )) - and e max := max xi∈N u(x i ) -1, where N denotes the set of nodes for the Lagrange basis functions (the symbol '-' means that the discrete maximum principle is actually satisfied), and (iii) the error on the global conservation property Φ(u h ) := | ∂Ω (β•n)u h ds| = 0. Note that the lack of exact global conservation for the linear problem is due to quadrature errors. Note also that we only impose weakly by the consistent penalty method the lower bound on the discrete solution; the upper bound could be imposed similarly.

For the above results on the nonlinear method with consistent penalty, the nonlinear system is solved using fixed-point iteration to an accuracy of T OL = 10 -6 on the increment u k hu k+1

h

Ω
where k is the index of the fixed-point iteration. This convergence criterion has been used for the sole purpose of numerical illustrations. In practice, a computationally-effective possibility is to prescribe the fixed-point convergence tolerance so that the error induced by the stopping criterion of the iteration is comparable to that of the a priori estimate. 7.43e-2 (-) 5.64e-1 (-) 9.94e-2 8.86e-2 5.68e-4 1 5.23e-2 (0.5) 5.04e-1 (0.2) 7.01e-2 1.13e-1 1.41e-4 2 2.69e-2 (1.0) 6.67e-1 (-0.4) 6.41e-2 6.77e-2 1.84e-5 3 1.59e-2 (0.8) 6.47e-1 (0.0) 6.01e-2 5.88e-2 1.20e-5 4 5.31e-3 (1.6) 2.95e-1 (1.1) 1.94e-2 2.0e-2 2.01e-6

, k = 1, γ = 0.0001h, τ = h/2
Assume that the fixed-point iteration is a contraction so that for some 0 < δ < 1, the approximation at iteration k of u h , say u k h , satisfies

u h -u k h Ω ≤ δ u h -u k-1 h Ω . Since u h -u k h Ω ≤ u h -u k+1 h Ω + u k h -u k+1 h Ω ≤ δ u h -u k h Ω + u k h -u k+1 h Ω , we infer that u h -u k h Ω ≤ (1 -δ) -1 u k h -u k+1 h Ω .
Then, we obtain the following error estimate for the error at the iteration k:

u k h -u Ω ≤ u h -u Ω + u h -u k h ≤ Ch k+ 1 2 + (1 -δ) -1 u k h -u k+1 h Ω .
Hence, if δ stays uniformly bounded away from 1 during the fixed-point iterations, we can reasonably stop the iterations whenever

u k h -u k+1 h Ω ∼ Ch k+ 1 2 ∼ e 0 /(2 l ) k+ 1 2
, where e 0 denotes the L 2 -norm error on the coarsest mesh and l denotes the level of mesh refinement. We call the resulting iterative method the balanced fixed-point iteration. We present the results using balanced fixedpoint iteration with T OL = 0.01/(2 l ) k+ 1 2 in Tables 8 and9. We observe that the results are of comparable or even better quality when the balanced fixedpoint iterations are used. In this context, it is also interesting to compare the behaviour of the present consistent penalty method with a method employing diffusion-based shock-capturing terms with artificial viscosity depending on the residual. We observe that for the present method, the computational cost is reduced as the accuracy of the solution improves and the violation of the DMP is isolated close to layers, whereas no such reduction is observed for nonlinear diffusion where the nonlinearity appears to have a much more global character. Indeed, it is known that the effects of the nonlinearity for shock-capturing can propagate into the zones where the solution is smooth in the form of gradient oscillations (this is sometimes called a terracing phenomenon). No such spurious gradient fluctuations were observed for the present consistent penalty method.

Test case 2: solution with discontinuity

In this section we apply the present consistent penalty method to approximate a solution with a discontinuity. The example that we consider is taken from [START_REF] Kuzmin | Goal-oriented mesh adaptation for flux-limited approximations to steady hyperbolic problems[END_REF]. The equation is similar to in the previous example, but this time we set

Figure 1 :

 1 Figure 1: Test case 1. Left: computational mesh. Right: contourlines of discrete solution for ǫ = 0.01

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure 2: Test case 2. Elevations of solutions using piecewise affine elements. Left: Linear GaLS, violation DMP 15%. Right: consistent penalty method with lumped quadrature and balanced fixed-point iteration, nodal violation DMP less than 4 • 10 -3 %

Table 3

 3

	: Test case 1, sharp layer (ǫ = 0.01). Nonlinear method with consistent
	penalty, k = 1, γ = 0.0001h, τ = h/2

Table 4 :

 4 Test case 1, sharp layer (ǫ = 0.01). Linear GaLS method, k = 1,

	τ = h/2
	l	e Ω	β•∇e Ω	e min	e max	Φ(u h )
	0	5.60e-2 (-)	6.47e-1 (-)	4.36e-3 7.58e-2 5.82e-3
	1 2.77e-2 (1.0) 7.17e-1 (-0.1) 5.25e-3 6.27e-2 2.47e-3
	2 1.05e-2 (1.4) 4.39e-1 (0.7) 1.81e-4 3.91e-2 1.06e-3
	3 2.86e-3 (1.9) 1.90e-1 (1.2) 4.07e-4 8.81e-3 1.10e-4
	4 3.96e-4 (2.9) 4.07e-2 (2.2) 2.31-19	-	3.78e-9

Table 5 :

 5 Test case 1, sharp layer (ǫ = 0.01). Nonlinear method with consistent penalty, fifth-order quadrature, k = 2, γ = 0.005h, τ = h/2

	l	e Ω	β•∇e Ω	e min	e max	Φ(u h )
	0	5.40e-2 (-)	6.27e-1 (-)	9.35e-6 7.66e-2 7.47e-3
	1 2.60e-2 (1.1) 6.68e-1 (-0.1) 8.87e-6 6.21e-2 2.83e-3
	2 1.03e-2 (1.3) 4.50e-1 (0.6)	-	3.80e-2 1.37e-3
	3 2.86e-3 (1.8) 1.90e-1 (1.2)	-	8.82e-3 1.20e-4
	4 3.96e-4 (2.9) 4.07e-2 (2.2)	-	-	3.78e-9
	Table 6: Test case 1, sharp layer (ǫ = 0.01). Nonlinear method with consistent
	penalty, hybrid quadrature, k = 2, γ = 0.0001h, τ = h/2
	l	e Ω	β•∇e Ω	e min	e max	Φ(u h )
	0	5.33e-2 (-)	6.58e-1 (-)	9.26e-2 7.68e-2 3.34e-4
	1 2.56e-2 (1.1) 7.29e-1 (-0.1) 7.82e-2 6.22e-2 5.84e-5
	2 1.05e-2 (1.3) 4.70e-1 (0.6)	4.17e-2 3.83e-2 2.07e-6
	3 2.93e-3 (1.8) 1.92e-1 (1.3)	1.03e-2 8.82e-3 5.31e-8
	4 3.96e-4 (2.9) 4.07e-2 (2.2) 4.76e-18	-	3.78e-9
	Table 7: Test case 1, sharp layer (ǫ = 0.01). Linear GaLS method, k = 2,
	τ = h/2

Table 8 :

 8 Test case 1, sharp layer (ǫ = 0.01). Nonlinear method, balanced fixedpoint iteration

l e Ω β•∇e Ω e min e max Φ(u h ) 0 5.30e-2 (-) 6.56e-1 (-) -7.68e-2 9.08e-3 1 2.60e-2 (1.0) 7.30e-1 (-0.15) -6.21e-2 1.83e-3 2 1.03e-2 (1.3) 4.63e-1 (0.7) -3.80e-2 1.46e-3 3 2.87e-3 (1.8)

1.98e-1 (1.2) -8.82e-3 1.20e-4 4 3.96e-4 (2.9) 4.07e-2 (2.3) --3.78e-9

Table 9: Test case 1, sharp layer (ǫ = 0.01). Nonlinear method, hybrid quadrature, balanced fixed-point iteration, k = 2, γ = 0.0001h, τ = h/2 Ω = (-1, 1) × (0, 1), take β = (y, -x) T (so that ∂Ω -= (-1, 0) × {0} ∪ {-1} × (0, 1) ∪ (0, 1) × {1}) and use the inflow data g = 1 on (-0.65, -0.35) × {0}, 0 elsewhere on ∂Ω -.

The corresponding exact solution reads

We compute the approximate solutions to this problem on structured meshes with mesh-size h = 0.1×2 -l , l ∈ {0, . . . , 4}. First, using the linear GaLS method and piecewise linear or quadratic finite elements, we recorded violations of positivity of more than 14% for the former case and more than 11% in the latter on all meshes. Using the present consistent penalty method with lumped quadrature for linear elements and the hybrid low-order quadrature (using vertices and midpoints) for quadratic elements resulted in strictly nodally positive solutions in all cases. Balanced fixed-point iterations were also used and we observed convergence after two iterations, yielding an accuracy similar to that of the linear method. In Figure 2, we report an illustration for the piecewise linear case with h = 1/20. In Figure 3 we report the contour plots of the nodal interpolant of the negative part of the solution. The magnitude of the DMP violation in the linear case is of the order 15% and in the nonlinear case of the order 4 • 10 -3%. The elevation of the linear GaLS method is presented in the left panel and that of the present method in the right one, here a term eliminating local overshoots has been added as well. In Figure 4 and 5, we report the same results, but this time for piecewise quadratic approximation and h = 1/10. The magnitude of the DMP violation in the linear case is of the order 20% and in the nonlinear case of the order 4 • 10 -3%. Observe that in spite of the small nodal DMP violation, the piecewise quadratic approximation is observed to violate the DMP with up to 10% due to the non-positivity of the quadratic basis functions.