J Ér

Ôme Droniou

Robert Eymard
email: robert.eymard@u-pem.fr

HIGH-ORDER MASS-LUMPED SCHEMES FOR NONLINEAR DEGENERATE ELLIPTIC EQUATIONS

Keywords: Stefan problem, porous medium equation, nonlinear degenerate elliptic equations, numerical scheme, mass-lumping, gradient discretisation method, error estimate, Finite Elements, Discontinuous Galerkin AMS subject classification: 65N12, 65N15, 65N30, 35J70

We present and analyse a numerical framework for the approximation of nonlinear degenerate elliptic equations of the Stefan or porous medium types. This framework is based on piecewise constant approximations for the functions, which we show are essentially necessary to obtain convergence and error estimates. Convergence is established without regularity assumption on the solution. A detailed analysis is then performed to understand the design properties that enable a scheme, despite these piecewise constant approximations and the degeneracy of the model, to satisfy high-order error estimates if the solution is piecewise smooth. Numerical tests, based on continuous and discontinuous approximation methods, are provided on a variety of 1D and 2D problems, showing the influence on the convergence rate of the nature of the degeneracy and of the design choices.

Introduction

The goal of numerical methods for partial differential equations is to approximate, as accurately as possible, the continuous solution. For mesh-based methods, it is well-known that when the problem is linear and the solution has sufficient regularity properties, for a fixed number of degrees of freedom high-order methods provide more accurate solutions than low-order methods. This result must however be questioned in the case of nonlinear problems for which, even if the solution is smooth enough, stability and high-order estimates might not be achievable without the proper structure of the chosen discretisation. We propose in this work to explore this question, considering the following nonlinear degenerate elliptic equation as the basis of our discussion:

β(u) -div(Λ∇ζ(u)) = f + div(F )
in Ω,

ζ(u) = 0 on ∂Ω, (1.1)
for which the corresponding weak formulation is Find u ∈ L 2 (Ω) such that ζ(u) ∈ H 1 0 (Ω) and

Ω β(u)v + Ω Λ∇ζ(u) • ∇v = Ω f v - Ω F • ∇v , ∀v ∈ H 1 0 (Ω).
(1.2)

Throughout the paper, we denote by • L 2 the norms in L 2 (Ω) or L 2 (Ω) d , and we make the following assumptions:

• Ω is an open bounded connected subset of R d (d ∈ N ) and T > 0, (1.3a)

• ζ : R → R is continuous and non-decreasing, ζ(0) = 0 and, for some M 0 , M 1 > 0, |ζ(s)| ≥ M 0 |s| -M 1 for all s ∈ R, (1.3b)

• β : R → R is continuous and non-decreasing, β(0) = 0 and, for some

K 0 , K 1 > 0, |β(s)| ≤ K 0 |s| + K 1 for all s ∈ R, (1.3c)
• β + ζ : R → R is strictly increasing, (1.3d) This work was supported by the Australian Government through the Australian Research Council's Discovery Projects funding scheme (project DP170100605), and by the French Government through the Agence Nationale de la Recherche (project CHARMS, ANR-16-CE06-0009).

• Λ : Ω → M d (R) is measurable and there exists λ ≥ λ > 0 such that, for a.e. x ∈ Ω, Λ(x) is symmetric with eigenvalues in [λ, λ].

(1.3e)

• f ∈ L 2 (Ω) and F ∈ L 2 (Ω) d .

(1.3f)

Remark 1.1 (On assumptions). Given that ζ(u) only appears through its gradient in (1.1), assuming ζ(0) = 0 is not restrictive. Likewise, upon changing f into f -β(0), the assumption β(0) = 0 is also not restrictive.

The theoretical study of (1.1) is covered by the pioneering paper [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF], extended in [START_REF] Andreianov | Well-posedness results for triply nonlinear degenerate parabolic equations[END_REF], on problems also including a nonlinear convection term; using techniques necessiting to multiply the equation by various functions of the unknown, existence and uniqueness of an entropy solution are obtained. An existence and uniqueness result for the simpler problem considered here is given by Theorem A.1 in Appendix A, without refering to entropy solutions.

The case ζ = Id fits into quasilinear second-order elliptic problems, the approximation of which is covered by a rather large literature, see e.g. [START_REF] Milner | Mixed finite element methods for quasilinear second-order elliptic problems[END_REF][START_REF] Chen | Expanded mixed finite element methods for quasilinear second order elliptic problems[END_REF][START_REF] Zeng | Error estimates of the lumped mass finite element method for semilinear elliptic problems[END_REF][START_REF] Bi | Global superconvergence and a posteriori error estimates of the finite element method for second-order quasilinear elliptic problems[END_REF]. The case ζ = Id, on which we focus in this paper, raises severe issues and is less often considered in the literature, especially when considering the question of high-order schemes. First, for such a problem, the solution can display discontinuities when ζ has plateaux. Moreover, the nonlinearities challenge the design of numerical methods that simultaneously (i) only require to compute integrals of polynomials (integrals that can be exactly computed in general), (ii) are amenable to error estimates (or, at the very least, proven to be convergent), and (iii) are of order higher than 1.

Extending the entropy method used in [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF] to the notion of entropy process solutions, the convergence of a Two-Point Flux Approximation (TPFA) finite volume method is proved in [START_REF] Eymard | Convergence of a finite volume scheme for nonlinear degenerate parabolic equations[END_REF] for a time-dependent version of (1.1) with Λ = Id. The entropy method requires to consider φ(u) as test functions for various nonlinear functions φ, a process that can only be reproduced at the discrete level for the TPFA scheme, see [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF]Section 7] and [START_REF] Eymard | TP or not TP, that is the question[END_REF]. Unfortunately, the TFPA scheme is only applicable on very specific grids, which usually forces Λ = Id, and only low-order error estimates can be expected from the application of the doubling variable technique as in [START_REF] Eymard | Error estimate for approximate solutions of a nonlinear convectiondiffusion problem[END_REF].

In the general case of an anisotropic heterogeneous field Λ, we need to consider more versatile schemes than the TPFA scheme, which will necessarily reduce the range of admissible test functions. Nevertheless, an important feature to preserve, if one wants to ensure the stability of the discretisation, is the capacity to choose appropriate test functions to simultaneously get diffusion estimates from the gradient terms, and a positive sign from the reaction term.

Let us first consider the case of conforming Galerkin methods. Given a subspace V h of H 1 0 (Ω), a conforming scheme for (1.2) is written Find u ∈ V h such that:

Ω β(u)v + Ω Λ∇ζ(u) • ∇v = Ω f v - Ω F • ∇v , ∀v ∈ V h .
(1.4)

If u ∈ V h we have ζ(u) ∈ H 1 0 (Ω) and the key of the convergence analysis is that the chain rule ∇ζ(u) = ζ (u)∇u enables us to take v = u ∈ V h as test function in the scheme. This choice creates from the diffusion term the quantity of interest ζ (u)|∇u| 2 , while the reaction term is non-negative since β(u)u ≥ 0. However, to deduce any sort of estimate from this choice of test function, we are forced to set F = 0, since the term Ω F • ∇u cannot in general be estimated using Ω ζ (u)|∇u| 2 . A better choice of test function to estimate the term resulting from the presence of F would be v = ζ(u) since the diffusion term would provide the quantity Ω |∇ζ(w)| 2 , which can be used to estimate Coming back to the general case of (1.2) with possibly F = 0, we consider numerical methods for which the chain rule does not hold at the discrete level (as is the case for the majority of nonconforming methods). The only reasonable test function to consider in order to get estimates is then v = ζ(u), which formally provides |∇ζ(u)| 2 from the diffusion term. More precisely, let us consider a scheme where the discrete unknowns z = (z i ) i∈I represent pointwise values of the solutions at certain nodes, and functions z h are reconstructed from these values and used in the weak formulation (this is the choice made in [START_REF] Elliott | Error analysis of the enthalpy method for the Stefan problem[END_REF][START_REF] Amiez | Error estimates for Euler forward scheme related to two-phase Stefan problems[END_REF] in the case of the transient problem, through the use of "Lagrange interpolation operators"). Then, for u = (u i ) i∈I , one can easily define v = ζ(u) pointwise, setting v i = ζ(u i ) for all i ∈ I. The weak formulation then involves ∇(ζ(u)) h • ∇v h and, taking v = ζ(u), this diffusion term generates the quantity

|∇(ζ(u)) h | 2 .
With this choice of v, the reaction term creates the quantity β(u h )(ζ(u)) h . This function is, at the considered nodes, equal to β(u i )ζ(u i ) ≥ 0 (see (1.3b)-(1.3c)). However, outside the nodes, no particular sign can be ensured for β(u h )(ζ(u)) h and it is not clear that this reaction term will indeed lead to proper estimates on the solution to the scheme.

The way to solve this conundrum is, in the reaction term, to use a different reconstruction of functions than the natural reconstruction (•) h used for the diffusion term. Utilising for example a piecewise constant reconstruction, in which the only values taken by the reconstruction are nodal values, ensures that the positivity of the reaction term -valid at these nodal values -extends to the entire domain (this is again done for low-order methods in [START_REF] Elliott | Error analysis of the enthalpy method for the Stefan problem[END_REF][START_REF] Amiez | Error estimates for Euler forward scheme related to two-phase Stefan problems[END_REF] to handle the accumulation terms issued from the time derivative). For linear models, using piecewise constant reconstructions for reaction/accumulation terms leads to what is called mass-lumped schemes. There is a large literature on the mass-lumping of Finite Element methods for second order problems, see e.g. [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF][START_REF] Geevers | New higher-order mass-lumped tetrahedral elements for wave propagation modelling[END_REF][START_REF] Rösch | Mass lumping for the optimal control of elliptic partial differential equations[END_REF][START_REF] Jund | Arbitrary high-order finite element schemes and high-order mass lumping[END_REF] and references therein. In most of these references, though, the construction of mass-lumped versions of high-order methods is justified by a need to reduce computational costs: for explicit discretisations of time-dependent linear problems, a mass-lumped scheme ensures a diagonal mass matrix which, unlike the standard mass matrix, is trivial to invert at each time step. This property of diagonal mass matrix has also been heavily used in schemes for eigenvalues problems related to linear elliptic operators (see for example [START_REF] Andreev | Some results in lumped mass finite-element approximation of eigenvalue problems using numerical quadrature formulas[END_REF] and references therein). On the contrary, for a nonlinear degenerate model as (1.1), as explained above the mass-lumping is not just a way to improve the method's efficiency, but appears as an imperative to establish convergence and error estimates -and thus rigorously ensure that the scheme has high-order approximation properties. Additionally, the usual interpretation of mass-lumping as a specific choice of quadrature rules for the mass matrix is mostly meaningful in the linear setting. For nonlinear models, the less standard interpretation based on piecewise constant reconstructions is more appropriate (even though, as we will see, there is still some link to exploit with local quadrature rules). Finally, let us notice that, to our best knowledge, mass-lumping techniques seem to only be considered in the literature on Finite Element methods, not in the literature covering other high-order polynomial-based methods such as Discontinuous Galerkin. This is understandable when the goal is to simplify the inversion of the mass matrix; mass-lumping is then not much useful to methods such as Discontinuous Galerkin schemes, for which the standard mass matrix is easy to invert due to its block diagonal structure (which can also easily be made fully diagonal by a simple choice of orthogonal local polynomial basis). However, when the primary objective of mass-lumping is to enable convergence and error estimates for nonlinear models, the question of designing mass-lumped Discontinuous Galerkin (or other methods based on local polynomials) is fully relevant.

Our goal in this paper is to design high-order mass-lumped schemes for the nonlinear degenerate model (1.1). Our main contributions can be summarised as follows:

• design of a general analysis framework that treats in a unified way many different methods, including Finite Elements and Discontinuous Galerkin methods (and others); • proof of error estimates in this general framework;

• identification of conditions on the mass-lumping to ensure high-order convergences (when the exact solution is piecewise smooth), despite the nonlinearities and degeneracy in the model; • extensive numerical tests, using both P k Finite Elements and Discontinuous Galerkin schemes, on realistic test cases (porous, Stefan) to validate the theoretical analysis.

Let us describe the organisation of this paper. We first provide in Section 2 a general formulation of numerical schemes, based on schemes written in fully discrete form: approximate functions and gradients are reconstructed without direct relation, and the approximate functions are piecewise constant. This construction is performed in the gradient discretisation method [START_REF] Droniou | The gradient discretisation method[END_REF], a framework that provides efficient notations and notions for the design and analysis of such schemes. After proving a first convergence result (Theorem 2.8) in Section 2.1, we establish in Section 2.2 error estimates on the approximation of ζ(u) when using mass-lumped schemes (Theorem 2.11 and Corollary 2.14). As demonstrated in Section 2.3, this general error estimate yields a high-order convergence rate (Theorem 2.23) for piecewise smooth solutions to (1.2), provided the mass-lumping is performed in a way that corresponds to sufficiently high-order local quadrature rules. These conditions on the local quadrature rules are similar to those highlighted for Finite Elements in [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF][START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF] but, interestingly, they appear here from the need of estimating quite different error terms than in the case of linear models as in these references. Extensive numerical tests are presented in Section 3, both on porous medium equations and on Stefan problems, using mass-lumped Finite Element and Discontinuous Galerkin schemes; the results confirm that high-order approximations are obtained only if the aforementioned local quadrature rules hold, even if the theoretical assumptions are not fully satisfied (e.g. the solution is not piecewise smooth). The paper is completed with a short conclusion (Section 4) and two appendices. In Appendix A, the properties of the continuous problem are analysed, and Appendix B sketches the study of the conforming scheme (1.4) with F = 0, and highlights its limitations compared to the method in Section 2: strong convergence of the gradients only under some regularity assumption on the continuous solution, no error estimate, no uniqueness of the discrete solution.

Schemes with piecewise constant approximation

To present the discretisation of (1.1), we use the gradient discretisation method (GDM) [START_REF] Droniou | The gradient discretisation method[END_REF], a generic numerical analysis framework for diffusion equations that encompasses many different discretisations: finite element, finite volumes, etc. Using this framework enables a unified treatment of all these different schemes, and also gives efficient setting and tools to deal with them, including the notion of mass-lumping that will be essential to design a scheme for which an error estimate can be established.

The principle of the GDM is to introduce discrete elements -a finite dimensional space, an operator that reconstructs functions, and an operator that reconstructs gradients -together called a gradient discretisation (GD), and to replace the continuous counterparts in the weak formulation (1.2) with these discrete elements, leading to a gradient scheme (GS) for (1.1).

Definition 2.1 (Gradient discretisation). A gradient discretisation is

D = (X D ,0 , Π D , ∇ D , Q D ) such that • X D ,0 a finite-dimensional space. • Π D : X D ,0 → L 2 (Ω) and ∇ D : X D ,0 → L 2 (Ω) d are
linear operators reconstructing, respectively, a function and a gradient; ∇ D must be chosen such that

• D := ∇ D • L 2 is a norm on X D ,0 . • Q D : L 2 (Ω) → L 2 (Ω) is a quadrature operator.
Remark 2.2 (Quadrature operator). Quadrature rules for source terms are usually not accounted for in the definition and analysis of gradient schemes. In the context of mass-lumped schemes, however, accounting for quadrature rules is essential to establishing optimal high-order error estimates.

Note that Q D is not assumed to be linear, which enables different choices of quadrature rules depending on the regularity subspace of the considered functions (for example, if f is continuous/piecewise continuous functions, Q D f could be defined using pointwise values of f , whereas for more irregular f the quadrature function Q D f could involve average values of f ).

To deal with the nonlinearity in the derivatives in (1.1) we need the following notion. Definition 2.3 (Piecewise constant reconstruction). Let D be a gradient discretisation such that, for some finite sets I and I ∂Ω ⊂ I, it holds

X D ,0 = {v = (v i ) i∈I : v i ∈ R ∀i ∈ I , v i = 0 ∀i ∈ I ∂Ω }.
We say that the reconstruction Π D is piecewise constant if there exists a partition U = (U i ) i∈I of Ω (some of the U i can be empty) such that

∀v = (v i ) i∈I ∈ X D ,0 , Π D v = i∈I v i 1 Ui ,
where 1 Ui is the characteristic function of U i . In other words, (Π D v) |Ui = v i for all i ∈ I.

In the setting of this definition, if g : R → R is a function satisfying g(0) = 0, we define (with an abuse of notation) g : X D ,0 → X D ,0 by applying g coefficient by coefficient:

∀v = (v i ) i∈I , g(v) = (g(v i )) i∈I . (2.1)
We note that this definition actually depends on the choice of the basis of X D ,0 . In practice, this basis being canonical and chosen once and for all, we do not make explicit the dependency of g(v) with respect to this basis. If Π D is a piecewise constant reconstruction, then (2.1) leads to

∀v ∈ X D ,0 , Π D g(v) = g(Π D v). (2.2)
The accuracy properties of a GD are assessed through the following quantities. The first one measures a discrete Poincaré constant of D, the second one is an interpolation error, whilst the last one measures the conformity defect of the method (how well a discrete divergence formula holds).

C D := max v∈X D ,0 \{0} Π D v L 2 v D , (2.3
)

∀ϕ ∈ H 1 0 (Ω) , S D (ϕ) = min v∈X D ,0 ∇ D v -∇ϕ L 2 + Π D v -ϕ L 2 , (2.4
)

∀ψ ∈ H div (Ω) , W D (ψ) := max v∈X D ,0 \{0} 1 v D Ω ∇ D v • ψ + Π D vdivψ . (2.5)
In the following, unless otherwise specified, the notation a b means that a ≤ Cb with C > 0 depending only on the data in Assumption (1.3) and on an upper bound of C D .

Given a gradient discretisation D = (X D ,0 , Π D , ∇ D , Q D ) with piecewise constant reconstruction as in Definition 2.3, the gradient scheme for (1.1) is (compare with the weak formulation (1.2)):

Find u ∈ X D ,0 such that Ω β(Π D u)Π D v + Ω Λ∇ D ζ(u) • ∇ D v = Ω Q D f Π D v - Ω F • ∇ D v , ∀v ∈ X D ,0 . (2.6)
Remark 2.4 (Quadrature for div(F )). A quadrature operator could also be introduced for F (for example, considering Q D component-wise, or selecting a different quadrature operator more appropriate to the structure of the gradient reconstruction). For simplicity of the presentation we decide not to include it in the analysis.

2.1. Convergence analysis. We first prove an a priori estimate on the solution to the gradient scheme. This estimate is used to prove the existence of this solution, its convergence, and the error estimate (2.10).

Lemma 2.5 (Bounds on the solution to the GS). Let D be a GD with piecewise constant reconstruction as in Definition 2.3, and let u ∈ X D ,0 be a solution to the gradient scheme (2.6). Then

Π D u L 2 + Π D β(u) L 2 + ζ(u) D Q D f L 2 + F L 2 .
(2.7)

Proof. Letting v = ζ(u) in (2.6) we get Ω β(Π D u)ζ(Π D u) + Ω Λ∇ D ζ(u) • ∇ D ζ(u) = Ω Q D f Π D ζ(u) - Ω F • ∇ D ζ(u),
where we have used

λ ∇ D ζ(u) 2 L 2 ≤ Q D f L 2 Π D ζ(u) L 2 + F L 2 ∇ D ζ(u) L 2 ≤ C D Q D f L 2 + F L 2 ∇ D ζ(u) L 2 . Recalling that ζ(u) D = ∇ D ζ(u) L 2 ,
ζ(u) L 2 Q D f L 2 + F L 2 . By (2.
2) this gives an L 2 (Ω)estimate on ζ(Π D u) and, using assumption (1.3b), translates into the bound on Π D u in (2.7). The estimate on β(Π D u) follows from the sub-linearity of β stated in Assumption (1.3c).

Lemma 2.6 (Existence and uniqueness for the GS). Assume (1.3) and let D be a gradient discretisation with a piecewise constant reconstruction as in Definition 2.3. Then there exists a solution to the gradient scheme (2.6) and, if (u 1 , u 2 ) are two solutions to this scheme, then ζ(u 1 ) = ζ(u 2 ) and

Π D u 1 = Π D u 2 .
Remark 2.7 (Counter-example to u 1 = u 2 ). In general, we cannot claim that u 1 = u 2 , as the following counter-example shows. Consider β(s) = s and ζ : R → R such that ζ(s) = 0 for all s ∈ [0, 1]. Take F = 0 and f ∈ L 2 (Ω) such that 0 ≤ f ≤ 1 almost everywhere, and consider an HMM gradient scheme [START_REF] Droniou | The gradient discretisation method[END_REF]Chapter 13] on a polytopal mesh of Ω. Denoting by M and F, respectively, the sets of cells and faces of this mesh, the corresponding gradient discretisation satisfies

X D ,0 = {v = ((v K ) K∈M , (v σ ) σ∈F ) : v K ∈ R , v σ ∈ R , v σ = 0 if σ ⊂ ∂Ω} and (Π D v) |K = v K for all K ∈ M. We select Q D = Id,
and the precise expression of ∇ D v is irrelevant to our counter-example. Then u = ((u K ) K∈M , (u σ ) σ∈F ) ∈ X D ,0 defined by

u K = 1 |K| K f ∀K ∈ M , u σ ∈ [0, 1] ∀σ ∈ F , u σ = 0 if σ ⊂ ∂Ω
is solution to the gradient scheme (2.6). Indeed, all the components of such a vector belong to [0, 1] and thus ζ(u) = 0. The scheme equation on u thus reduces to

Ω Π D uΠ D v = Ω f Π D v , ∀v ∈ X D ,0 ,
which holds given the choice of the cell values (u K ) K∈M .

Proof of Lemma 2.5. The existence is obtained via a topological degree argument; we refer the reader to [START_REF] Deimling | Nonlinear functional analysis[END_REF] for the definition and properties of this degree. Fix an arbitrary Euclidean structure, with inner product •, • , on the finite dimensional space X D ,0 . For a

∈ [0, 1] let ζ a (s) = aζ(u) + (1 -a)u. Define F : [0, 1] × X D ,0 → X D ,0 the following way: for a ∈ [0, 1] and u ∈ X D ,0 , F(a, u) is the unique element of X D ,0 such that, for all v ∈ X D ,0 , F(a, u), v = Ω aβ(Π D u)Π D v + Ω Λ∇ D ζ a (u) • ∇ D v - Ω aQ D f Π D v - Ω aF • ∇ D v.
We note that u is a solution to the gradient scheme (2.6) if and only if F(1, u) = 0. By continuity of β and ζ, and the finite dimension of X D ,0 , the mapping F is clearly continuous. Assume that F(a, u) = 0 for some a ∈ [0, 1]. The arguments in the proof of Lemma 2.5, using v = ζ a (u) as a test function, show that ζ a (u) D ≤ C 1 with C 1 not depending on a; by equivalence of norms on the finite dimensional space X D ,0 , this shows that ζ a (u) ∞ ≤ C 2 with C 2 still independent on a and • ∞ the supremum norm in X D ,0 on an arbitrary basis. The mapping ζ a satisfies (1.3b) with M 0 , M 1 independent of a. As a consequence, the bound on ζ a (u) ∞ shows that u ∞ < R with R independent of t.

Hence, any solution to F(a, u) = 0 lies in the open ball B R of X D ,0 , centered at 0 and of radius R in the • ∞ norm. This ball being independent of a, the topological degree theory ensures that deg

(F(1, •), B R , 0) = deg(F(0, •), B R , 0). The mapping F(0, •) : X D ,0 → X D ,
0 is linear and the estimate obtained on the solutions to F(0, u) = 0 shows that F(0, •) has a trivial kernel, and is therefore invertible. This implies deg(F(0, •), B R , 0) = 0 and thus deg(F(1, •), B R , 0) = 0, which proves that the equation F(1, u) = 0 has a solution u ∈ B R .

We now consider the uniqueness of the solution to the scheme. Subtracting the equations satisfied by u 1 and u 2 and taking v = ζ(u 1 ) -ζ(u 2 ) ∈ X D ,0 as a test function, we have

Ω (β(Π D u 1 ) -β(Π D u 2 ))Π D (ζ(u 1 ) -ζ(u 2 )) + Ω Λ∇ D (ζ(u 1 ) -ζ(u 2 )) • ∇ D (ζ(u 1 ) -ζ(u 2 )) = 0.
Property (2.2) and the monotonicity of β and ζ show that

(β(Π D u 1 ) -β(Π D u 2 ))Π D (ζ(u 1 ) -ζ(u 2 )) = (β(Π D u 1 ) -β(Π D u 2 ))(ζ(Π D u 1 ) -ζ(Π D u 2 )) ≥ 0. Hence, ∇ D (ζ(u 1 ) -ζ(u 2 )) L 2 = 0 which, by property of ∇ D , ensures that ζ(u 1 ) = ζ(u 2 ).
We now come back to the equations satisfied by u 1 and u 2 , subtract them and take

v = β(u 1 ) - β(u 2 ) ∈ X D ,0 as a test function to get Ω (β(Π D u 1 ) -β(Π D u 2 )) 2 + Ω ∇ D (ζ(u 1 ) -ζ(u 2 )) • ∇ D (β(u 1 ) -β(u 2 )) = 0. Since ζ(u 1 ) = ζ(u 2 ), we infer that β(Π D u 1 )-β(Π D u 2 ) = 0. Owing to Hypothesis (1.3d), we conclude that Π D u 1 = Π D u 2 from β(Π D u 1 ) + ζ(Π D u 1 ) = β(Π D u 2 ) + ζ(Π D u 2 ).
The next theorem is our first main convergence result. It states the strong convergence of the solution to the gradient scheme without assuming any regularity property on the continuous solution.

Theorem 2.8 (Convergence of the scheme). Let (D m ) m∈N be a sequence of GDs, with piecewise constant reconstructions as in Definition 2.3, that satisfies the following properties:

• (Coercivity) The sequence (C D m ) m∈N is bounded, where C D m is defined by (2.3) for D = D m . • (Consistency) Recalling the definition (2.4), there holds ∀ϕ ∈ H 1 0 (Ω), lim m→∞ S D m (ϕ) = 0 and lim m→∞ Q D m f -f L 2 = 0. (2.8)
• (Limit-conformity) Recalling the definition (2.5), there holds

∀ψ ∈ H div (Ω) , lim m→∞ W D m (ψ) = 0.
(2.9)

• (Compactness) For any (v m ) m∈N such that v m ∈ X D m ,0 for all m ∈ N and (∇ D m v m ) m∈N is bounded in L 2 (Ω) d , the set {Π D m v m : m ∈ N} is relatively compact in L 2 (Ω).
For any m ∈ N, let u m be a solution of Scheme (2.6). Then there exists u solution to

(1.2) such that, as m → ∞, Π D m ζ(u m ) → ζ(u) strongly in L 2 (Ω), ∇ D m ζ(u m ) → ∇ζ(u) strongly in L 2 (Ω) d , and Π D m β(u m ) → β(u) weakly in L 2 (Ω).
Proof. Using Estimate (2.7) and (2.8) (which shows that Q D m f L 2 is bounded), the compactness and the limit-conformity of (D m ) m∈N , [START_REF] Droniou | The gradient discretisation method[END_REF]Lemma 2.15] gives Z ∈ H 1 0 (Ω) and B ∈ L 2 (Ω) such that, up to a subsequence (not made explicit in the following),

Π D m ζ(u m ) → Z strongly in L 2 (Ω), ∇ D m ζ(u m ) → ∇Z weakly in L 2 (Ω) d and β(Π D m u m ) → B weakly in L 2 (Ω). By weak/strong convergence we infer that lim m→∞ Ω Π D m β(u m )Π D m ζ(u m ) = Ω ZB.
The monotonicity properties β and ζ then enable us to apply [START_REF] Droniou | The gradient discretisation method[END_REF]Lemma D.10] (a Minty's trick) to get u ∈ L 2 (Ω) such that Z = ζ(u) and B = β(u).

We now show that u solves (1.2). Let ϕ ∈ H 1 0 (Ω) and let v m ∈ X D m ,0 be an element that realises the minimum defining S D m (ϕ). By (2.8) (2.8), enable us to take the limit m → ∞ of the scheme to see that u is a solution to (1.2). The uniqueness of u (see Theorem A.1) shows that the convergence properties holds for the whole sequence (u m ) m∈N instead of just along the subsequence previously extracted.

we have Π D m v m → ϕ in L 2 (Ω) and ∇ D m v m → ∇ϕ in L 2 (Ω) d . Use v m as a test function in GS (2.6) satisfied by u m . The convergence properties of β(Π D m u m ) and ∇ D m ζ(u m ) towards B = β(u) and ∇Z = ∇ζ(u), together with the convergence Q D m f → f in L 2 (Ω) stated in
It remains to establish the strong convergence of

∇ D m ζ(u m ). We now let m → +∞ in the GS (2.6) with D = D m and v = ζ(u m ), that is, Ω β(Π D m u m )Π D m ζ(u m ) + Ω Λ∇ D m ζ(u m ) • ∇ D m ζ(u m ) = Ω Q D m f Π D m ζ(u m ) - Ω F • ∇ D ζ(u m ). This yields lim m→∞ Ω Λ∇ D m ζ(u m ) • ∇ D m ζ(u m ) = Ω (f ζ(u) -F • ∇ζ(u) -β(u)ζ(u)) = Ω Λ∇ζ(u) • ∇ζ(u),
where the conclusion follows using v = u in (1.2). Since (ξ, η) → Ω Λξ • η is an inner product on L 2 (Ω), this relation and the weak convergence of (∇

D m ζ(u m )) m∈N imply the strong convergence of ∇ D m ζ(u m ) to ∇ζ(u) in L 2 .
2.2. Error estimate. The analysis above pinpoints the required structure on a numerical scheme to ensure proper bounds and convergence of the solution -namely, the piecewise constant reconstruction property. We now want to establish error estimates to better assess this convergence. In practice, one usually starts from a given numerical method and would like to apply it to the model under consideration. Following our discussion above, if the given method does not have a piecewise constant function reconstruction, it has to be modified into a method that has such a reconstruction. This process is called the mass-lumping of the original scheme. In the context of the GDM, this notion is translated in the following definition. The following theorem states a general error estimate on the gradient scheme (2.6).

Definition 2.9 (Mass-lumped GD). Let D * = (X D ,0 , Π D * , ∇ D , Q D ) be a gradient discretisation. A gradient discretisation D is a mass-lumped version of D * if it only differs from D * through the function reconstruction (that is, D = (X D ,0 , Π D , ∇ D , Q D )),
Theorem 2.11 (Error estimate for the GS). Assume (1.3) and let D be a mass-lumped version of a D * in the sense of Definition 2.9. Let u be a solution to the gradient scheme (2.6), and let u be the solution to (1.2) (see Theorem A.1). Then, for any

I D ζ(u) ∈ X D ,0 , there holds ∇ D [I D ζ(u) -ζ(u)] L 2 W D * (Λ∇ζ(u) + F ) + ∇ D I D ζ(u) -∇ζ(u) L 2 + R D ,D * (u, f ) + T D (u, u), (2.10) where R D ,D * (u, f ) = max v∈X D ,0 \{0} 1 ∇ D v L 2 Ω Π D v[β(Q D u) -Q D f ] -Π D * v[β(u) -f ] , (2.11)
and

T D (u, u) = max Ω β(Q D u) -β(Π D u) Π D I D ζ(u) -ζ(Q D u) ; 0 1/2 . ( 2
∇ D ζ(u) -∇ζ(u) L 2 ≤ ∇ D [ζ(u) -I D ζ(u)] L 2 + ∇ D I D ζ(u) -∇ζ(u) L 2 .
Similarly, introducing ±Π D * I D ζ(u), using the triangle inequality and the definition of C D * , we have Proof of Theorem 2.11. Since u is the solution to (1.2), we have div(Λ∇ζ(u

Π D * ζ(u) -ζ(u) L 2 ≤ Π D * [ζ(u) -I D ζ(u)] L 2 + Π D * I D ζ(u) -ζ(u) L 2 ≤ C D * ∇ D [ζ(u) -I D ζ(u)] L 2 + Π D * I D ζ(u) -ζ(u) L 2 . An estimate on ∇ D [ζ(u) -I D ζ(u)]
)+F ) = β(u)-f ∈ L 2 (Ω). Hence, by definition (2.5) of W D applied to D * , for any v ∈ X D ,0 , v D W D * (Λ∇ζ(u) + F ) ≥ Ω ∇ D v • (Λ∇ζ(u) + F ) + Π D * vdiv(Λ∇ζ(u) + F ) = Ω ∇ D v • Λ∇ζ(u) + F • ∇ D v + Π D * v[β(u) -f ].
Substituting the term involving F using (2.6), we get

Ω Λ∇ D v • (∇ζ(u) -∇ D ζ(u)) + (Q D f -β(Π D u))Π D v + Π D * v[β(u) -f ] ≤ v D W D * (Λ∇ζ(u) + F ). (2.13) Introducing ±Λ∇ D v • ∇ D I D ζ(u) and ±β(Q D u)Π D v
in the left-hand side, using the Cauchy-Schwarz inequality and recalling that v D = ∇ D v L 2 , we infer

Ω Λ∇ D v • (∇ D I D ζ(u) -∇ D ζ(u)) + (β(Q D u) -β(Π D u))Π D v + Ω Π D * v[β(u) -f ] -Π D v[β(Q D u) -Q D f ] ∇ D v L 2 W D * (Λ∇ζ(u) + F ) + ∇ D I D ζ(u) -∇ζ(u) L 2 . (2.14) Choose v = I D ζ(u) -ζ(u). Introducing ±ζ(Q D u) and using the monotonicity of ζ and β (which yields [β(b) -β(a)][ζ(b) -ζ(a)] ≥ 0 for all a, b ∈ R) together with (2.2), we have (β(Q D u) -β(Π D u))Π D v = [β(Q D u) -β(Π D u)] [Π D I D ζ(u) -ζ(Q D u)] + [β(Q D u) -β(Π D u)] [ζ(Q D u) -ζ(Π D u)] ≥ [β(Q D u) -β(Π D u)] [Π D I D ζ(u) -ζ(Q D u)].
Plugging this into (2.14) and using (1.3e) leads to

∇ D [I D ζ(u) -ζ(u)] 2 L 2 ∇ D v L 2 W D * (Λ∇ζ(u) + F ) + ∇ D I D ζ(u) -∇ζ(u) L 2 + Ω Π D v[β(Q D u) -Q D f ] -Π D * v[β(u) -f ] + Ω [β(Q D u) -β(Π D u)] [Π D I D ζ(u) -ζ(Q D u)] ∇ D v L 2 W D * (Λ∇ζ(u) + F ) + ∇ D I D ζ(u) -∇ζ(u) L 2 + R D ,D * (u, f ) + T D (u, u) 2 .
Using the Young inequality on the first term in the right-hand side and recalling that

v = I D ζ(u)-ζ(u) leads to ∇ D [I D ζ(u) -ζ(u)] 2 L 2 W D * (Λ∇ζ(u) + F ) + ∇ D I D ζ(u) -∇ζ(u) L 2 + R D ,D * (u, f ) 2 + T D (u, u) 2 .
The proof of (2.10) is complete taking the square root of this estimate and using

√ a 2 + b 2 ≤ a + b for all a, b ≥ 0.
From the general estimate (2.10) we deduce the following bound on the error, which is often leads to (low order) rates of convergence as noted in Remark 2.15. This estimate will be improved, for situations corresponding to classical mass-lumping versions of schemes with nodal interpolants, in Section 2.3.

Corollary 2.14. Under the assumptions of Theorem 2.11, define

α D ,D * = max v∈X D ,0 \{0} Π D v -Π D * v L 2 ∇ D v L 2
and let I D ζ(u) be given by

I D ζ(u) = argmin v∈X D ,0 ∇ D v -∇ζ(u) L 2 + Π D v -ζ(u) L 2 . Then, ∇ D [I D ζ(u) -ζ(u)] L 2 W D * (Λ∇ζ(u) + F ) + S D (ζ(u)) + α D ,D * + β(Q D u) -β(u) L 2 + Q D f -f L 2 + S D (ζ(u)) + ζ(u) -ζ(Q D u) L 2 1 2 , (2.15)
where the hidden multiplicative constant in additionally depends on

β(Q D u) L 2 and Q D f L 2 . Remark 2.
(Q D u) -β(u) L 2 , Q D f -f L 2 and ζ(u) -ζ(Q D u) L 2 .
The estimate (2.15) can thus be expected, most of the time, to provide an

O(h 1
2 ) rate of convergence, the limiting factor in the right-hand side of (2.15) being the last one, coming from T D (u, u). We will see in Section 2.3 that this estimate is however very pessimistic and, in many cases, can be improved to higher powers of h (see Remark 2.26).

Proof. We estimate each term, except the first one, in the right-hand side of (2.10). By choice of

I D ζ(u) and definition (2.4) of S D , ∇ D I D ζ(u) -∇ζ(u) L 2 + Π D I D ζ(u) -ζ(u) L 2 = S D (ζ(u)).
Hence the term

∇ D I D ζ(u) -∇ζ(u) L 2 in (2.10) is bounded above by S D (ζ(u)).
Using the definition of α D ,D * and of C D , we have

R D ,D * (u, f ) ≤ α D ,D * β(u) -f L 2 + max v∈X D ,0 \{0} 1 ∇ D v L 2 Ω Π D v[β(Q D u) -Q D f ] -Π D v[β(u) -f ] α D ,D * + C D [β(Q D u) -Q D f ] -[β(u) -f ] L 2 α D ,D * + β(Q D u) -β(u) L 2 + Q D f -f L 2 .
This gives the third, fourth and fifth terms in the right-hand side of (2.15). For the last term in this estimate, we write, using the Cauchy-Schwarz inequality and the a priori bound (2.7) on β(Π D u),

T D (u, u) 2 ≤ β(Q D u) -β(Π D u) L 2 Π D I D ζ(u) -ζ(Q D u) L 2 ( β(Q D u) L 2 + Q D f L 2 + F L 2 ) Π D I D ζ(u) -ζ(u) L 2 + ζ(u) -ζ(Q D u) L 2 .
The proof is complete taking the square root and recalling that Π

D I D ζ(u)-ζ(u) L 2 ≤ S D (ζ(u)).
2.3. Suitable quadrature rules lead to high-order estimates. Let us first make the following broken regularity assumption on the data and solution.

Assumption 2.16 (Data and exact solution). F = 0 and, u being the solution to (1.2) and s ≥ 1 being an integer, f and β(u) belong to the broken Sobolev space

W s,∞ (M) := g ∈ L ∞ (Ω) : g |K ∈ W s,∞ (K) ∀K ∈ M .
This space is endowed with the norm g W s,∞ (M) := max K∈M g W s,∞ (K) .

Remark 2.17 (Piecewise continuity and local smoothness). W s,∞ (M) is a subspace of

C(M) := {g ∈ L ∞ (Ω) : g |K ∈ C(K) ∀K ∈ M}. ( 2
β(u) ∈ C(M). By (1.3b)-(1.3d), β + ζ is an homeomorphism of R. Hence, u = (β + ζ) -1 (β(u) + ζ(u))
and thus u ∈ C(M).

In the rest of this section, we consider a slightly more precise setting than in Section 2.2. We assume that D * has a piecewise polynomial function reconstruction (possibly of high-order) and unknowns associated to nodes on the domain, and that specific local quadrature rules can be chosen. Typically, P k or Q k finite elements and Symmetric Interior Penalty discontinuous Galerkin (SIPG) schemes, with mass-lumping constructed using dual meshes around the nodes, fit into this setting. In what follows, h X denotes the diameter of a set X ⊂ R d . (1) (Mesh) M be a polytopal mesh of Ω ⊂ R d with d ≤ 3, in the sense of [START_REF] Droniou | The gradient discretisation method[END_REF]Definition 7.2] (this definition actually represents the mesh as a quadruple of sets of cells, faces, points and vertices, that will not be useful to our purpose; we therefore confuse the mesh with the set of cells). The mesh size is h = max K∈M h K . (2) (Space) There is a finite set I, partitioned into I Ω and I ∂Ω , such that

X D ,0 = {v = (v i ) i∈I : v i ∈ R , ∀i ∈ I , v i = 0 , ∀i ∈ I ∂Ω } . (3) (Local polynomial reconstructions) There is a polynomial degree k ≥ 1 such that, for all K ∈ M and all v ∈ X D ,0 , (Π D * v) |K ∈ P k . (4) (Broken gradient bound) There is C ∇ ≥ 0 such that, for all v ∈ X D ,0 , ∇ h (Π D * v) L 2 ≤ C ∇ ∇ D v L 2
, where ∇ h is the usual broken gradient on M. (5) (Nodes) There is a family (x i ) i∈I of points in Ω, and subsets (I K ) K∈M of I, such that I = (∪ K∈M I K ) ∪ I ∂Ω and, for all v = (v i ) i∈I ∈ X D ,0 , all K ∈ M and all i ∈ I K , we have

x i ∈ K and v i = (Π D * v) |K (x i ). Additionally, x i ∈ ∂Ω whenever i ∈ I ∂Ω . (6) (Mass-lumping) D = (X D ,0 , Π D , ∇ D , Q D ) is a mass-lumped version of the gradient discretisation D * = (X D ,0 , Π D * , ∇ D , Q D )
in the sense of Definition 2.9, which means that Π D is piecewise constant on a partition U = (U i ) i∈I in the sense of Definition 2.3. We further assume that

U i ∩ K = ∅ only if i ∈ I K . (7) (Interpolant) I D ζ(u) ∈ X D ,0 is given by the nodal values of ζ(u), that is, (I D ζ(u)) i = ζ(u)(x i ) for all i ∈ I. This is well-defined since ζ(u) ∈ C(Ω) (see Remark 2.17). (8) (Quadrature rule) The quadrature Q D is defined on C(M) (see (2.16)) by ∀g ∈ C(M) , ∀K ∈ M , (Q D g) |K = i∈I K g |K (x i )1 Ui∩K .
(2.17) (9) (Mesh regularity) There exists ρ > 0 such that: Any K ∈ M is star-shaped with respect to all points in a ball of radius ρh K , For all i ∈ I, ρh Ui ≤ h.

A few remarks are of order.

Remark 2.20 (Local polynomial space). The space P k in Item (3) could be replaced by any of its subspace P K that contains P 1 ; the analysis would not be hindered, and some assumptions could even be weakened (see Remark 2.25). We chose to use P k to simplify the presentation.

Remark 2.21 (Nodes). The same i can belong to several I K , as is the case for conforming finite elements. Conversely, in the case of DG scheme for example, the following may occur (see the numerical example in Section 3.3):

• one can have x i = x j for i = j, • I K does not necessarily contain all the indices i ∈ I such that x i ∈ K, • there can exist i ∈ I ∂Ω \ (∪ K∈M I K ) -but, in that case, U i = ∅.
Remark 2.22 (Quadrature rule). The gradient scheme (2.6) is usually implemented by assembling cell contributions. When the source term f is continuous on each cell, and since Π D v is constant on each U i , it is customary to use the simple -apparently low-order -quadrature rule defined by (2.17). Since u ∈ C(M) by Remark 2.18, the formula (2.17) can thus be used to compute Q D u. The definition of Q D outside C(M) is irrelevant to the analysis that follows, and could for example be taken as the identity outside this space.

In the rest of this section, we write a b as a shorthand for "a ≤ Cb with C not depending on M or U , but possibly depending on ρ, k and C ∇ ". Theorem 2.23 (High-order error estimate). Under Assumption 2.19, let ≥ 0 be an integer and suppose that we have the following exactness of local quadrature rules (where k is the degree in Item (3) of Assumption 2.19):

∀K ∈ M , ∀q ∈ P k+ , K q = i∈I K |U i ∩ K|q(x i ).
(2.18)

Let s ∈ {0, . . . , + 2} be such that Assumption 2.16 holds. Then, the solution u to (2.6) satisfies

∇ D [I D ζ(u) -ζ(u)] L 2 W D * (Λ∇ζ(u)) + ∇ D I D ζ(u) -∇ζ(u) L 2 + h s (1 + C D * ) β(u) -f W s,∞ (M) . (2.19)
Let us first make a few remarks before proving the theorem.

Remark 2.24 (Quadrature rule). The quadrature rule (2.18) bears similarities with the conditions on quadrature rules highlighted for Finite Elements in [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF][START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF]. However, in the proof below, the exactness condition (2.18) responds to a different need than the ones encountered in the analysis of mass-lumped Finite Elements for linear equations.

Remark (i) Term T D (u, u). For all K ∈ M, all i ∈ I K , and all x ∈ U i ∩ K, by Definition 2.3 of Π D and Item [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] in Assumption 2.19,

(h min(s,k) ) estimate on ∇ D [I D ζ(u) -ζ(u)] L 2 ,
Π D I D ζ(u)(x) = (I D ζ(u)) i = ζ(u)(x i ) = (ζ(u)) |K (x i ) = ζ(u |K (x i )) = ζ(Q D u(x)). Hence, Π D I D ζ(u) = ζ(Q D u) and T D (u, u) = 0. (ii) Term R D ,D * (u, f ).
For the sake of brevity, set g = β(u) -f . By definition (2.17

) of Q D , we have Q D g = β(Q D u) -Q D f
and thus, to bound R D ,D * (u, f ) above by the last term in (2.19), we have to establish that, for all v ∈ X D ,0 ,

Ω (Q D g Π D v -g Π D * v) h s (1 + C D * ) g W s,∞ (M) ∇ D v L 2 . ( 2

.20)

Let A D ,D * (g, v) be the integral in the left-hand side of (2.20). We have

A D ,D * (g, v) := K∈M i∈I K |U i ∩ K|g |K (x i )v i - K g Π D * v = K∈M i∈I K |U i ∩ K|g |K (x i )(Π D * v) |K (x i ) - K g Π D * v = K∈M E K (gΠ D * v), (2.21)
where, in the second line, we have used (Π [START_REF] Jund | Arbitrary high-order finite element schemes and high-order mass lumping[END_REF]), and we have defined the error in the local quadrature rule on K by

D * v) |K (x i ) = v i (see Item (5) in Assumption 2.
∀w ∈ C(K) , E K (w) := i∈I K |U i ∩ K|w |K (x i ) - K w.
By (2.18) and a straightforward estimate, ∀q ∈ P k+ , E K (q) = 0, and (2.22)

∀w ∈ C(K) , |E K (w)| ≤ 2|K| w L ∞ (K) . (2.23)
For a polynomial degree r ≥ 0, let Pr r K : L 2 (K) → P r (K) denote the L 2 (K)-orthogonal projector on P r and notice that, since (Π D * v) |K ∈ P k (Item (3) in Assumption 2.19) and k ≥ 1,

q := (Pr K g)(Π D * v) |K + (Pr 0 K (Π D * v) |K )(Pr +1 K g -Pr K g) (2.24)
belongs to P +k + P 0+ +1 ⊂ P k+ . Using (2.22) with this q yields

E K (gΠ D * v) = E K gΠ D * v -(Pr K g)(Π D * v) |K -(Pr 0 K (Π D * v) |K )(Pr +1 K g -Pr K g) = E K [g -Pr K g][(Π D * v) |K -Pr 0 K (Π D * v) |K ] + (Pr 0 K (Π D * v) |K )[g -Pr +1 K g] .
Invoking then the bound (2.23) and the straightforward estimate Pr 0

K (Π D * v) |K L ∞ (K) ≤ Π D * v L ∞ (K) , we infer |E K (gΠ D * v)| ≤ 2 g -Pr K g L ∞ (K) |K| (Π D * v) |K -Pr 0 K (Π D * v) |K L ∞ (K) + 2|K| Π D * v L ∞ (K) g -Pr +1 K g L ∞ (K) .
(2.25)

Under item 9 of Assumption 2.19, [10, Lemma 3.4] shows that, for any natural numbers a ≥ 0 and b ∈ {0, . . . , a + 1}, and any w ∈ W b,∞ (K),

w -Pr a K w L ∞ (K) h b K w W b,∞ (K)
. Applying this estimate with (a, b, w) = ( , min(s, + 1), g), (a, b, w) = (0, 1, (Π D * v) |K ) and (a, b, w) = ( + 1, s, g), (2.25) leads to

|E K (gΠ D * v)| h min(s, +1) K g W min(s, +1),∞ (K) |K|h K ∇(Π D * v) |K L ∞ (K) + |K| Π D * v L ∞ (K) h s K g W s,∞ (K) . The discrete inverse Lebesgue embedding of [10, Lemma 5.1] gives, if q ∈ P k (K), |K| q L ∞ (K) |K| 1 2 q L 2 (K)
. Applied to q = (Π D * v) |K and q = components of ∇(Π D * v) |K , and since min(s, + 1) + 1 = min(s + 1, + 2) ≥ s, we obtain

|E K (gΠ D * v)| h s K g W s,∞ (K) |K| 1 2 ∇(Π D * v) |K L 2 (K) + Π D * v L 2 (K) .
Plugging this estimate into (2.21), using a discrete Cauchy-Schwarz inequality on the sums, and recalling Item (4) in Assumption 2.19, we obtain

|A D (g, v)| h s g W s,∞ (M) ∇ D v L 2 + Π D * v L 2 .
The estimate (2.20) follows recalling the definition (2.3) of C D * .

Numerical illustrations

The tests are performed within the framework of Assumption 2.19, using mass-lumped versions of either conforming P k finite elements, or SIPG method. Each mass-lumped gradient discretisation D used in the tests share the same (X D ,0 , ∇ D , I D , Q D ) as the corresponding D * and thus, when describing D we focus on the particular choices of nodes (x i ) i∈I and partition (U i ) i∈I , as these choices dictate the exactness of the local quadrature rule (2.18).

3.1.

Setting for the tests. The convergences are assessed through the following quantities:

E Π β,I D = β(Q D u) -Π D β(u) L 2 (Ω) , E Π ζ,I D = Π D (I D ζ(u) -ζ(u)) L 2 (Ω) , E ∇ ζ,I D = ∇ D (I D ζ(u) -ζ(u)) L 2 (Ω) , E ∇ ζ = ∇ζ(u) -∇ D ζ(u) L 2 (Ω)
. measuring approximation errors on β(Q D u), the interpolation of ζ(u) (for both function and gradient reconstruction), and on ∇ζ(u) using high-order quadrature rules. A first order polynomial fit is done on the logarithms or these errors with respect to -1 d log(Card(I)), which yields an approximation under the form E CCard(I) -α/d . Our outputs give the numerical values of C and α, the latter providing a numerical convergence order with respect to an evaluation of the mesh size (the number of unknowns, Card(I), growing linearly with the number of cells).

All the 1D or 2D tests refer to the following situations: β = Id, and ζ ∈ {Id, ζ p , ζ s }, where the "porous media" ζ p function is defined by ∀s ∈ R, ζ p (s) = max(s, 0) 2 , and the "Stefan" ζ s function is defined by

∀s ∈ R, ζ s (s) =      s if s < 0, 0 if 0 ≤ s ≤ 1, s -1 if 1 < s.
In all the numerical tests, the approximate solution remains numerically bounded. There is therefore no need to re-define ζ p on the negative axis in order to explicitly satisfy the super-linear bound in Assumption (1.3b). Let us now give the complete continuous cases which are approximated below in 1D or in 2D.

Test case R: Regular problem, f = 0. This problem corresponds to ζ = Id (the model is therefore linear, but we still apply the mass-lumping process) and, for x ∈ (0, 1), the source term and solution are given by f (x) = 4xe x and u(x) = x(1 -x)e x .

Test case P-1: Porous media problem, homogeneous Dirichlet BC, f = 0. This test is on the porous medium equation, with ζ = ζ p . The source term and exact solutions are defined as follows: for x ∈ (0, 1), setting y x = max(x -0.2, 0) and z x = max(0.8 -x, 0), we take

f (x) = (y x z x ) 3/2 -6y x z x (z 2
x -3y x z x + y 2 x ) and u(x) = (y x z x ) 3/2 . To describe u, we first let γ ∈ (0, 1 2 ) such that u(x) = f (x) for g(x) ∈ (γ, 1 2 ) and u(x) ≥ 1 for g(x) ∈ (0, γ). The ODE in (1.1) can then be solved on each sub-interval and gives ζ(u(x)) = 0 for g(x) ∈ (γ, 1 2 ) and, for some a, b ∈ R, ζ(u(x)) = ae g(x) + be -g(x) + 3( 1 2 -g(x)) -1 for g(x) ∈ (0, γ). Hence, u(x) = ae g(x) + be -g(x) + 3( 1 2 -g(x)) for g(x) ∈ (0, γ). These values a, b and γ are found by expressing the matching conditions ensuring that ζ(u) ∈ H 2 (0, 1) (since

Test

(ζ(u)) = u -f ∈ L 2 (0, 1)), namely ζ( 1 2 ± γ) = 0 and ζ ( 1 2 ± γ) = 0;
the symmetry of the problem also imposes ζ ( 12 ) = 0. This leads to the following equations: 3 1 2 -γ -1 + ae γ + be -γ = 0, -3 + ae γ -be -γ = 0, and -3 + a -b = 0.

Numerically solving this nonlinear system of equations gives γ 0.33036, a 1.2545 and b -1.7455.

Test case S-2: Stefan problem, non-homogeneous Dirichlet BC, f = 0. As in the previous test, ζ = ζ s . The source term is fixed at f = 0 and, for any γ ∈ [0, 1], a solution is given by:

u(x) = 1 2 (e x-γ + e -(x-γ) ) = cosh(x -γ) ∀x ∈ (γ, 1), 0
∀x ∈ (0, γ). Non-homogeneous Dirichlet conditions are imposed at x = 1 to match the value of u there, and the tests are run with γ = 1 3 . Test case S-3: Stefan problem, homogeneous Dirichlet BC, f = 0 and F = 0. We let ζ = ζ s and

f (x) = 0 F (x) = 4 sinh(1/4) cosh(1/4) u(x) = 0 ∀x ∈ (0, 1 4 ), f (x) = 5 F (x) = 0 u(x) = 5 -4 cosh(x-1/2) cosh(1/4) ∀x ∈ ( 1 4 , 3 4 ), f (x) = 0 F (x) = -4 sinh(1/4) cosh(1/4) u(x) = 0 ∀x ∈ ( 3 4 , 1).
(3.1) 3.2. Mass-lumped finite elements. For a conforming simplicial mesh M and using the notations in Assumption 2.19, the gradient discretisations D * = D k,fe * , for k ∈ {1, 2, 3}, corresponding to conforming P k finite elements is defined by the following elements. • The points (x i ) i∈I are:

If k = 1: the mesh vertices (in dimension 1 or 2), If k = 2: the mesh vertices and one point in each cell (in dimension 1), or the mesh vertices and the edge midpoints (in dimension 2), If k = 3: the mesh vertices and two points in each cell (in dimension 1), or the mesh vertices, two additional points on each edge, and one point in each cell. I ∂Ω is the set of indices i ∈ I such that x i ∈ ∂Ω and, for K ∈ M, I K := {i ∈ I :

x i ∈ K}. • For each simplex K ∈ M and v = (v i ) i∈I ∈ X D ,0 , (Π D * v) |K is the unique polynomial in P k that takes the values v i at x i for all i ∈ I K . • The gradient reconstruction is given by (∇ D v) |K = ∇(Π D * v) |K for all v ∈ X D ,0 and K ∈ M.

3.2.1.

Numerical tests for mass-lumped finite elements in dimension 1. We consider two families of meshes of Ω = (0, 1) with N cells each, for N ∈ {16, 32, 64, 512, 1024, 2048}. The first one is the uniform mesh M u N with mesh step h = 1/N . The second one is a random mesh M r N such that each cell has size h i = H i / j H j , where H i = (3 + ρ i ) with ρ i following a random uniform law on (0, 1). As mentioned above, all the gradient discretisations are mass-lumped versions of the corresponding P k GD. We describe hereafter the remaining elements to fully define the GD, that is: the degree k, the nodes (x i ) i∈I , and the partition U = (U i ) i∈I . In each case, the elements of this partition are intervals and satisfy Item (6) in Assumption 2.19.

• For k = 1: D 1,fe u (uniform mesh) and D 1,fe r (random mesh) are defined setting I = {0, . . . , N }, (x i ) i∈I are the vertices of the mesh, so that (x i-1 , x i ) i=2,...,N are the cells of the mesh. For each i ∈ I, if x i is an endpoint of K then |U i ∩ K| = 1 2 |K| (providing an exact quadrature on K for polynomials of degree 1).

• For k = 2: D 2,fe u (uniform mesh) and D 2,fe r (random mesh) are defined setting I = {0, . . . , 2N } and the cells of the uniform or random mesh are (x 2j-2 , x 2j ) for j = 1, . . . , N . The remaining nodes are located at the centers of the cells, that is, x 2j-1 = 1 2 (x 2j-2 + x 2j ) for all j = 1, . . . , N . For all i ∈ I, if x i is an endpoint of K then |U i ∩ K| = 1 6 |K| and, if x i is the midpoint of K, |U i ∩ K| = 2 3 |K| (providing an exact quadrature on K for polynomials of degree 3). • For k = 3: D 3,fe u,a (uniform mesh) and D 3,fe r,a (random mesh) are defined setting I = {0, . . . , 3N } and the cells of the uniform or random mesh are (x 3i-3 , x 3i ) for i = 1, . . . , N . The remaining nodes are located at 1/3 and 2/3 inside each cell, that is,

x 3i-2 = 2 3 x 3i-3 + 1 3 x 3i and x 3i-1 = 1 3 x 3i-3 + 2 3 x 3i for i = 1, . . . , N . For all i ∈ I, if x i is an endpoint of K then |U i ∩ K| = 1 6 |K|, and if x i is inside K then |U i ∩ K| = 1
3 |K| (providing an exact quadrature on K only for polynomials of degree 1). • For k = 3: D 3,fe u,b (uniform mesh) and D 3,fe r,b (random mesh) are defined setting I and (x i ) i∈I as for D 3,fe u,a and D 3,fe r,a , For all i ∈ I,

if x i is an enpoint of K then |U i ∩ K| = 1 8 |K|, and if x i is inside K then |U i ∩ K| = 3
8 |K| (providing an exact quadrature on K for polynomials of degree 3). • For k = 3: D 3,fe u,c (uniform mesh) and D 3,fe r,c (random mesh) are defined using the Gauss-Lobatto method, setting I = {0, . . . , 3N } and the cells of the uniform or random mesh are (x 3i-3 , x 3i ) for i = 1, . . . , N . The remaining nodes are located at barycentric coordinates ( 5±

√ 5 10 , 5∓ √ 5 10 ) inside each cell, that is, x 3i-2 = 5+ √ 5 10 x 3i-3 + 5- √ 5 10 x 3i and x 3i-1 = 5- √ 5 10 x 3i-3 + 5+ √ 5 10 x 3i for i = 1, . . . , N . For all i ∈ I, if x i is an endpoint of K then |U i ∩ K| = 1
12 |K|, and if x i is inside K then |U i ∩ K| = 5 12 |K| (providing an exact quadrature on K for polynomials of degree 5). Remark 3.1 (Local quadratures). The exactness of the local quadratures mentioned in the presentations above show that (2.18) is satisfied with ≥ 0 by all the GDs above, except D 3,fe u,a and D 3,fe r,a . For α = u, r, the local quadratures for D 2,fe α (resp. D 3,fe α,c ) even satisfy (2.18) with = 1 (resp. = 2). as expected since this is the optimal rate, when using piecewise P k polynomials, to approximate a smooth non-polynomial function. These rates for the approximation of the gradient are actually above those predicted by our analysis: by Remark 3.1, for all these methods we can take = 0 in Theorem 2.23 and thus, following Remark 2.26, the expected decay of E ∇ ζ,I D is h min(2,k) = h k . Regarding k = 3, the schemes based on the D 3,fe α,c variant appear to have a worse rate than D 3,fe α,a or D 3,fe α,b , but focusing on the constant C we notice that they are actually much better. They also clearly outperform when considering the gradient reconstruction. Focusing on the latter, the convergence rate O(h) for D 3,fe α,a can be explained recalling that this variant does not even satisfy (2.18); even though T D (u, u) = 0 for this method (see the proof of Theorem 2.23), we do not have any better estimate on R D ,D * (u, f ) than in the proof of Corollary 2.14, which was precisely expected to be O(h) (see Remark 2.15).

E Π β,I D E Π ζ,I D E ∇ ζ,I D E ∇ ζ GD C α C α C α C α D 1,
On the contrary, referring again to Remark 3.1, D 3,fe α,b enables us to take = 0 in Theorem 2.23 and we recover the expected O(h 2 ) estimate on E ∇ ζ,I D mentioned in Remark 2.26. For D 3,fe α,c we can even take = 1 and Table 1 clearly shows that this leads to an improved and optimal O(h 3 ) estimate on the gradient (again, something predicted in Remark 2.26).

These results clearly demonstrate that the key factor in choosing a proper mass-lumped version for a high-order scheme is the exactness property (2.18) -not satisfying this property leads to decreased rates of convergence. They also indicate, at least for k = 3, the sharpness of the error estimate established in Theorem 2.23.

Test Case P-1: Porous media problem, homogeneous Dirichlet BC, f = 0 The results are presented in Table 2. The functions f and u are only piecewise smooth, and the discontinuity of their derivatives is not necessarily aligned with the mesh. As a consequence, Assumption 2.16 does not hold. Compared to the smooth case studied in Test Case R, the convergence is overall degraded. However, the rates mostly remain not far from the linear case (especially for gradient approximations), and the main features discussed for the smooth case can also be found here: there is a super-convergence of the function approximations (when using quadrature or interpolation of the exact solution), and the rates of convergence drop drastically if the local quadrature rule (2.18) do not hold with a high enough .

Test Case P-2: Porous media problem, non-homogeneous Dirichlet BC, f = 0 Table 3 details the outcomes of this test. The source term is obviously smooth, but the solution is only piecewise smooth. Despite this, the results show that, except for the very small constants observed for D 3,fe α,c , the schemes behave here in a very similar way as for the completely smooth 2.7e-01 2.40 5.2e-07 2.33 2.7e-04 3.10 9.9e-02 3.00 D 3,fe r,c 1.4e+00 2.76 2.8e-06 2.64 1.4e-03 3.46 1.1e-01 3.00 Table 3. Constants and rates for Test Case P-2.

E Π β,I D E Π ζ,I D E ∇ ζ,I D E ∇ ζ GD C α C α C α C α D 1,
E Π β,I D E Π ζ,I D E ∇ ζ,I D E ∇ ζ GD C α C α C α C α D 1,
situation of Test Case R. Here again we notice the importance of choosing proper local quadrature rules (2.18) when designing mass-lumped schemes from high-order methods.

Test Case S-1: Stefan problem, homogeneous Dirichlet BC, f = 0 The results for this test case are presented in Table 4. This test case is a much more severe one than the porous medium case, since the solution u is discontinuous. This explains the poor convergence of E Π β,I D for all considered methods. On the contrary, ζ(u) is continuous and E Π ζ,I D thus behaves much better, with an order 2 decay for all schemes. The order of decay of E ∇ ζ,I D is also similar for all methods (around 1.6), except for the GDs D 3,fe α,a for which it drops to 1; this reduction can be explained, as in the previous case, by recalling that these GDs do not satisfy the local quadrature rules (2.18) even for = 0.

Based on our previous discussion, we could expect the schemes corresponding to D 3,fe α,c to have a higher rate of convergence than the other methods, but it should be noted that ζ(u) only belongs to H 2 , not H 3 since (ζ(u)) = u -f is discontinuous. This limits the application of Theorem (2.23) to s = 2 (despite = 2 being a valid choice in this case), and leads to the h 2 rate for the approximation of the gradient.

We however notice that E ∇ ζ has a quite poor convergence (or does not seem to converge) on random meshes. Given that the difference between this error and E ∇ ζ,I D solely lies in the interpolation error

∇ D I D ζ(u) -∇ζ(u) L 2 ,
this apparently indicates that this interpolation error does not converge on random meshes. It is actually not the case, but for these meshes the regularity factor and maximum size oscillate a lot from one mesh to the other; combined with the low regularity of the solution (which implies a slow expected rate of convergence), this explains that the regression performed on

E Π β,I D E Π ζ,I D E ∇ ζ,I D E ∇ ζ GD C α C α C α C α D 1,fe u
1.8e+01 0.41 1.2e+01 1.97 1.2e+01 1.87 2.8e+00 1.00 D 1,fe r 3.7e+01 0.54 1.6e+01 2.15 3.2e+00 1.66 2.6e-01 -0.17 5.1e+01 0.61 2.9e+00 2.41 3.9e-01 1.59 4.7e-01 0.37 Table 4. Constants and rates for Test Case S-1. 1.3e-01 1.79 1.3e-01 1.79 7.9e-02 1.71 7.7e-06 -0.82 Table 5. Constants and rates for Test Case S-1, with excluding the discontinuities. the interpolation errors struggles to capture the correct convergence when considering a finite family of meshes.

E Π β,I D E Π ζ,I D E ∇ ζ,I D E ∇ ζ GD C α C α C α C α D 3 u,a
Figure 1 presents a loglog plot of the errors E ∇ ζ,I D vs. h, for k = 3 and when the errors are calculated excluding two intervals of length 2/10 centred at the discontinuity points 1 2 ± γ of u. The corresponding first order regression are given in Table 5 (which also includes E Π β,I D and E Π ζ,I D ). This figure shows why the regression fails to capture proper orders, as the behaviour of the error with respect to the mesh size is not monotonic enough. These results however show that selecting a highorder method with a proper local quadrature rule (e.g. D 3,fe u,c or D 3,fe r,c ) give much more precise results outside singularities than low-order methods or ill-chosen local quadrature rules.

Test Case S-2: Stefan problem, non-homogeneous Dirichlet BC, f = 0 The results presented in Table 6, are comparable to the results obtained with a non-zero source term in Test Case S-1 (reduced convergence for D 3,fe α,a , limitation of the convergence for D 3,fe α,c dues to the limited regularity of ζ(u)). In this case, however, the gradient ∇ D I D ζ(u) of the interpolant seems to enjoy better convergence property even on random meshes, which preserve a reasonable convergence of E ∇ ζ .

Test Case S-3: Stefan problem, homogeneous Dirichlet conditions, f = 0 and F = 0 The term -Ω F • ∇ D v in the gradient scheme (2.6) is exactly computed, without numerical quadrature, using the relation

- Ω F • ∇ D v = - 1 0 F (s)(Π D * v) (s)ds = -4 sinh(1/4) cosh(1/4) Π D * v 1 4 + Π D * v 3 4 .
The outcome of the test can be seen in Table 7. We note that these data do not satisfy the assumptions of Theorem 2.23, and no high-order rate can therefore be expected. Actually, the solution displays a very low regularity since ζ(u) only belongs to H 1 , not even H 2 . This is represented in the results by the fact that, for each given error, the rates of convergence for all schemes are in the same range. We 3.1e+00 0.71 6.0e-02 2.16 7.8e-02 1.78 6.1e-02 1.52 Table 6. Constants and rates for Test Case S-2.

E Π β,I D E Π ζ,I D E ∇ ζ,I D E ∇ ζ GD C α C α C α C α D 1,fe u 2.0e+00 0.
notice also that, across the board, the schemes perform better on regular grids rather than random grids.

3.2.2.

Numerical tests for mass-lumped finite elements in dimension 2. In the following 2D cases, we consider the domain Ω = (0, 1) × (0, 1), the polynomial degrees k = 1, 2, and the following meshes (see Figure 3):

• Triangular meshes which are as equilateral as possible, with edge length 1/N for N ∈ {25, 50, 100}.

The gradient discretisations on these meshes will have the subscript "e", e.g., D k e . • Rectangular triangular meshes obtained by splitting N 2 squares in 2, for N ∈ {25, 50, 100}. We use the subscript "s" for these GDs, e.g., D k s . • Random meshes based on the three meshes mesh1 3, mesh1 4 and mesh1 5 from the FVCA5 benchmark [START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF]. The randomness is obtained moving the internal nodes by a uniform random factor, and we use the subscript "r" for these GDs, such as in D k r . Based on these meshes, we define the mass-lumped version of D k * (k = 1, 2), for each α ∈ {e, s, r}: • D 1,fe α : k = 1, (x i ) i∈I are the vertices of the mesh and, for all i ∈ I, U i is Donald dual cell around x i . • D 2,fe α : k = 2, (x i ) i∈I are the vertices and midpoints of the edges of the mesh, and, for all i ∈ I, if x i is a vertex of the mesh, U i = ∅, if x i is the midpoint of and edge of the mesh, U i is the union of one or two triangles obtained joining the edge midpoints with each of the centers of mass of the cells around the edge (see Figure 2, left) • D 1,fe α,1/4 : we divide each mesh triangle in four (see Figure 2, right), and apply the mass-lumped P 1 method on this sub-mesh. This gradient discretisation has the same unknowns as D 2,fe α , corresponding to k = 2 on the original mesh. matrix is therefore singular, which is of course an issue when considering explicit discretisations of time-dependent (even linear) problems; solving this issue requires the usage of enriched P 2 elements [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF]. However, in the context of implicit time stepping, or equivalently of stationary problems, this is not an issue since the stiffness matrix is always non-singular.

E Π β,I D E Π ζ,I D E ∇ ζ,I D E ∇ ζ GD C α C α C α C α D 1,
The case of stationary nonlinear degenerate equations such as (1.1) requires nonetheless an implementation trick. Since the diffusion term acts on ζ(u), in the nonlinear iterations the stiffness matrix is multiplied by ζ (u) which can vanish, and the diffusion term does not yield in itself a control of all the unknowns u i . It does however enable a control of the unknowns ζ(u) i = ζ(u i ). Even though these unknowns, especially for the Stefan problem, do not determine u entirely, this gives a way to implement the scheme in a non-singular way. Instead of writing an equation on (u i ) i∈I , we write an equation on ((u i ) i∈Ie , (ζ(u) i ) i∈Iv ), where I e is the set of indices corresponding to edge midpoints, and I v the set of indices corresponding to the vertices. The unknowns (u i ) i∈Ie are controlled by the mass-lumped reaction term, and the unknowns (ζ(u) i ) i∈Iv by the diffusion term. This implementation does not entirely determine a solution u to the scheme, only its values at the edge midpoints and the values of ζ(u) at the vertices, but this is expected given Lemma 2.6 and Remark 2.7.

Remark 3.3 (The case k = 3). If k = 3, it is possible to satisfy the local quadrature rules (2.18) with l = 0 (i.e. to have rules exact for third degree polynomials). This is done by fixing α ∈ (0, ( 3 44 ) 1/2 ) and making as nodes x i and proportion |U i ∩ K|/|K| of the weights the following choices:

• the vertices of the mesh, each one associated with proportion 3-44α 2 60(1-4α 2 ) ; • two points on each edge located at the barycentric coordinates ( 1 2 ±α, 1 2 ∓α) on the edge, associated with proportion 1 15(1-4α 2 ) ; • the centers of mass of the triangles, associated with proportion 9/20. To have local quadratures of degree four (that is, (2.18) with l = 1), one must set α 2 = 1/12, which leads to the negative weight proportion -1/60 at the vertices of the triangle, a situation which is incompatible with the mass-lumping setting. To properly mass-lump the P 3 finite elements while preserving their high-order, an enriched version of these elements must be considered [START_REF] Cohen | Higher order triangular finite elements with mass lumping for the wave equation[END_REF].

The data we consider in the following test case are the same as for the 1D case, using the diagonal as 1D coordinate. For example, if g is a solution or source term for a 1D test case, the solution or source term for the corresponding 2D case is computed by setting g(x, y) = f ((x + y)/ √ 2). All these 2D test cases therefore have non-homogeneous Dirichlet boundary conditions. All the considered gradient discretisations satisfy the local quadrature rules (2.18) with = 0. Accordingly, if the solution and source were smooth, rates of convergence for E ∇ ζ,I D should be O(h k ) for k = 1, 2 (see Remark 2.26). The results show that, for the porous medium case, we are above these rates for all meshes, except for the random mesh -probably more representative of genuine situations-where we are at these rates (or slightly above). As in the 1D case, the Stefan problem is more challenging and, probably due to the loss of regularity of the solution, the rates are a little bit worse. They do however remain at or above O(h) for k = 1, and only drop to around O(h 1.5 ) for k = 2.

Test with D To properly assess the interest of using a 2nd order scheme over a 1st order, we now look, on the same triangular mesh, to the outputs of D if k = 2, we compute the ratios r = E 2 /E 1 for all the tests on the three random meshes based on mesh1 3, mesh1 4 and mesh1 5. Assuming that the each error E k is of the form C k ( h h0 ) α k , where h 0 is the size of the reference mesh mesh1 3, the ratio between the two errors should be given by

C α C α C α C α C α C α E Π β,ID
r = E 2 E 1 = C 2 C 1 h h 0 α2-α1
.

Table 12 performs a C(h/h 0 ) α regression of the ratio r. Hence, the C values in this table can be considered as approximations of C2 C1 , and the α values as approximations of α 2 -α 1 . The results show a clear advantage (smaller C k , larger α k ) of the second order method over the first order method, and also that this advantage still holds, albeit reduced, for irregular (Stefan) test cases. for the SIPG method of order k is defined as follows.

• For each cell K ∈ M, points (x i ) i∈I K are chosen such that for each choice of real numbers (w i ) i∈I K there is a unique q ∈ P k such that q(x i ) = w i for all i ∈ I K . Then I = (∪ K∈M I K ) ∪ I ∂Ω is the family that gathers the indices of all these points for all the cells, and of all the boundary points where a jump is accounted for in the expression of ∇ D . • The gradient reconstruction is given by (

• For each K ∈ M and v = (v i ) i∈I ∈ X D ,0 , (Π D * v) |K is
∇ D v) |K = ∇(Π D * v) |K + S K (v) for all v ∈ X D ,0
and K ∈ M, where S K (v) is an appropriate stabilisation term accounting for the jumps appearing in the DG scheme (see [START_REF] Droniou | The gradient discretisation method[END_REF]Definition 11.1] for details).

We take k = 3 and consider the same families of uniform and random meshes of Ω = (0, 1) with N cells each as in Section 3.2.1. We describe hereafter the remaining elements to fully define the GD, that is: the nodes (x i ) i∈I and the partition U = (U i ) i∈I . In each case, the elements of this partition are intervals and satisfy Item (6) in Assumption 2.19. These elements follow closely the choices made for the 1D finite element meshes in Section 3.2.1.

• D 3,dg u,a (uniform mesh) and D 3,dg r,a (random mesh) are defined setting I = {0, . . . , 4N + 1} with

x 0 = x 1 = 0 < x 2 < x 3 < x 4 = x 5 < . . . x 4i = x 4i+1 < x 4i+2 < x 4i+3 < . . . x 4N = x 4N +1 = 1.
The cells of the uniform or random mesh are (x 4i+1 , x 4i+4 ) for i = 0, . . . , N -1. The nodes x 4i+2 and x 4i+3 are located at 1/3 and 2/3 inside each cell. Each cell i = 0, . . . , N -1 is partitioned into four intervals U 4i+1 , . . . U 4i+4 , associated to the points x 4i+1 , . . . x 4i+4 , with respective lengths 1/6, 1/3, 1/3, 1/6 of that of the cell (providing an exact quadrature (2.18) only for polynomials of degree 1), and U 0 and U 4N +1 are empty. • D 3,dg u,b (uniform mesh) and D 3,dg r,b (random mesh) are similar to the previous one. The only difference is that the respective lengths of the intervals are 1/8, 3/8, 3/8, 1/8 of that of the cell (providing an exact quadrature (2.18) for polynomials of degree 3). • D 3,dg u,c (uniform mesh) and D 3,dg r,c (random mesh) are defined using the Gauss-Lobatto method, the differences with D 3,dg α,a are that the nodes x 4i+2 and x 4i+3 are located at 5- √ 5 10 and 5+ √ 5 10 inside each cell, and that the respective lengths are 1/12, 5/12, 5/12, 1/12 of that of the cell (providing an exact quadrature (2.18) for polynomials of degree 5). In all the tests, the penalisation parameter (denoted by β in [START_REF] Droniou | The gradient discretisation method[END_REF]Chapter 11]) is fixed at 0.6. Remark 3.4 (Usage of the GDM). SIPG schemes are usually implemented by assembling consistent and stabilisation terms. The design of a proper stabilisation term for a nonlinear problem as considered in this work is not straightforward, and requires ad-hoc choices. When embedding the SIPG method into the GDM, the stabilisation term is already accounted for in the design of the gradient reconstruction ∇ D . No specific treatment of a separate stabilisation of the nonlinear second order term is thus required.

The results in Table 13 are in line with what we already observed for Finite Element methods. The better local quadrature rule of D 3,dg α,c enable these schemes to outperform D 3,dg α,b and D 3,dg α,a (this latter behind badly hindered by its very low-order local quadrature rule), and preserves the expected order 3 convergence for smooth data and solutions. This optimal convergence is even noticed in the fully nonlinear test case P-2. As before, the Stefan problem is much more challenging due to its reduced regularity, but even for this one we notice an interest in selecting a method with high enough local quadrature rules.

Figure 6 shows the solutions u obtained with FE and DG schemes, for k = 1, 3, on Test Case S-1 and with a relatively coarse mesh (N = 16). As expected, the solutions obtained with k = 3 are much more accurate. They however present oscillations (more severe for DG than for FE) in the viscinity of the discontinuity of u. The solutions for Test Case S-2, corresponding to f = 0, do not present such oscillations.

Conclusion

We presented a generic analysis framework, covering a range of methods, for the numerical approximation of nonlinear degenerate elliptic equations, stationary version of the Stefan or porous medium problems. We identified a particular structure of the method, the piecewise constant function reconstruction, which appears to be necessary to establish the robustness of the schemes, and to obtain error estimates. We showed how to design mass-lumping versions of high-order numerical methods in order to preserve, despite the usage of piecewise constant approximations in the scheme, high-order approximations of the solution to this severely nonlinear model. Our numerical tests on mass-lumped Finite Element and Discontinuous Galerking schemes corroborated the theoretical findings, showing that even for non-smooth solutions an elevated rate of convergence is only obtained if the mass-lumping is designed to satisfy proper local quadrature rules. We finally consider the regularity properties of ζ(u). This function is a weak solution of ζ(u) ∈ H 1 0 (Ω) and -div(Λ∇ζ(u) + F ) = f -β(u) ∈ L 2 (Ω). This readily shows that Λ∇ζ(u) + F ∈ H div (Ω). If d ≤ 3, then L 2 ⊂ W -1,q (Ω) for some q > d and thus, assuming that F ∈ L p (Ω) d for p > d, ζ(u) is a solution in H 1 0 (Ω) of -div(Λ∇ζ(u)) = f + div(F ) -β(u) ∈ W -1,min(q,p) (Ω); the results of [START_REF] Stampacchia | Le problème de Dirichlet pour les équations elliptiques du second ordre à coefficients discontinus[END_REF] then show that ζ(u) has the Hölder-regularity Remark B.2 (On condition (B.3)). We observe that (B.3) holds in the case where w ∈ H 1 0 (Ω) since it can then be taken as a test function in (B.1). But it may also hold in some less regular situations.

Proof. We only sketch the proof. Assuming the existence of a solution w m to the scheme, we let v = v m in (B. Letting ϕ ∈ H 1 0 (Ω) and taking v m := argmin v∈Vm ϕ -v H 1 0 (Ω) in (B.2), the above convergences enable us to take the limit as m → ∞ to see that w is the solution to (B.1). The uniqueness of w shows that the convergence property holds for the whole sequence.

Assuming that (B.

( 2 . 2 )

 22 to write Π D ζ(u) = ζ(Π D u) in the first integral term. By monotonicity of β, ζ and β(0) = ζ(0) = 0, we have β(s)ζ(s) ≥ 0 and the equation above thus gives, by definition of C D and Assumption (1.3e),

  and if Π D is a piecewise constant reconstruction in the sense of Definition 2.3. Remark 2.10. Of course, if D * already has a piecewise constant reconstructions as in Definition 2.3, one can take D = D * .

. 12 )

 12 Remark 2.12 (Choice of I D ζ(u)). The element I D ζ(u) can be any vector in X D ,0 . However, the estimate (2.10) is obviously useful only if ∇ D I D ζ(u) is close to ∇ζ(u). This is usually achieved selecting for I D ζ(u) a suitable interpolant of ζ(u), which is why we used this notation. Remark 2.13 (Approximation of ζ(u)). Introducing ±∇ D (I D ζ(u)) and using a triangle inequality, we have

  as in Theorem 2.11 therefore also yields an estimate on ∇ D ζ(u) -∇ζ(u) and Π D * ζ(u) -ζ(u), modulo the additional interpolation errors ∇ D I D ζ(u) -∇ζ(u) and Π D * I D ζ(u) -ζ(u). If D * has function and gradient reconstructions that are piecewise polynomial of high-order, these interpolation errors can be expected to have a high rate of convergence with respect to the mesh size. The same argument also gives an error estimate on Π D ζ(u) -ζ(u), but the corresponding interpolation error Π D I D ζ(u) -ζ(u) is limited to a first-order convergence since Π D is a piecewise constant reconstruction.

Assumption 2 . 19 (

 219 Structure of D * , D and I D ζ(u)).

  which is a drastic improvement over(2.15) (see Remark 2.15). The same O(h min(s,k) ) bound on ∇ D [I D ζ(u) -ζ(u)] L 2 holdsfor the gradient discretisation corresponding to DG schemes of degree k, provided that Λ∇ζ(u) ∈ H min(s,k) (Ω) d (see[START_REF] Droniou | The gradient discretisation method[END_REF] Lemmas 11.14 and 11.15]). Proof of Theorem 2.23. The inequality (2.19) follows from Theorem 2.11, estimating in the present context the terms T D (u, u) and R D ,D * (u, f ).

case P- 2 :

 2 Porous media problem, non-homogeneous Dirichlet BC, f = 0. Still taking for ζ the porous medium function ζ = ζ p , this test takes u(x) = max(x -1 5 , 0) 2 /12, which corresponds to the source term f = 0, and non-homogeneous Dirichlet boundary conditions are imposed on ζ(u). Test case S-1: Stefan problem, homogeneous Dirichlet BC, f = 0. In this test, the nonlinearity is given by the Stefan-like function ζ = ζ s . The source term is given by f (x) = 3( 1 2 -g(x)), where g(x) = | 1 2 -x|.

D 2,fe u 6 .

 6 0e+01 0.76 1.1e+00 2.04 6.2e-01 1.54 2.5e+00 1.61 D 2,fe r 3.2e+01 0.50 1.2e+00 2.04 1.0e+00 1.65 3.3e-03 -0.16 D 3,fe u,a 7.9e+01 0.84 1.2e+00 2.03 3.7e-01 1.03 4.4e-01 1.06 D 3,fe r,a 1.0e+02 0.83 1.2e+00 2.15 3.8e-01 1.05 3.6e-02 0.28 D 3,fe u,b 8.6e+01 0.84 3.8e-01 1.95 7.2e-01 1.61 8.9e-01 1.53 D 3,fe r,b 5.1e+01 0.64 3.4e-01 1.84 2.8e-01 1.53 1.4e+03 2.77 D 3,fe u,c 5.4e+01 0.67 4.6e-01 2.08 3.6e-01 1.58 8.5e-01 1.56 D 3,fe r,c

  1.2e+00 2.03 1.2e+00 2.03 3.0e-01 1.04 3.0e-01 1.04 D 3,fe r,a 5.9e-01 2.03 5.9e-01 2.03 4.7e-02 0.71 1.6e+02 2.08 D 3,fe u,b 3.3e-01 1.96 3.3e-01 1.96 1.5e+00 2.00 1.2e+00 1.94 D 3,fe r,b 1.9e-01 1.71 1.9e-01 1.71 1.5e+00 2.03 3.7e-01 0.66 D 3,fe u,c 3.9e-01 2.09 3.9e-01 2.09 1.1e-02 2.10 7.1e-05 0.49 D 3,fe r,c

Figure 1 .

 1 Figure 1. Test Case S-1: log(E ∇ ζ,I D ) vs. log(h), excluding discontinuities.

  50 2.6e-01 1.98 1.5e-01 1.48 7.7e-01 1.00 D 1,fe r 4.7e+00 0.67 5.2e-02 1.78 3.8e-02 1.32 7.6e-01 1.00 D 2,fe u 2.3e+00 0.49 1.2e-01 2.02 8.6e-02 1.50 2.0e-01 1.50 D 2,fe r 2.4e+00 0.53 1.6e-01 2.13 1.0e-01 1.61 2.1e-01 1.51 D 3,fe u,a 3.4e+00 0.50 9.3e-02 2.00 8.9e-02 1.01 9.2e-02 1.01 D 3,fe r,a 2.9e-02 -0.21 6.8e-02 1.94 8.5e-02 1.00 8.8e-02 1.00 D 3,fe u,b 4.1e+00 0.53 5.6e-02 2.03 8.0e-02 1.50 1.1e-01 1.50 D 3,fe r,b 1.7e+01 1.10 1.5e-01 2.28 1.8e-01 1.75 9.3e-02 1.51 D 3,fe u,c 3.1e+00 0.50 4.9e-02 2.01 5.3e-02 1.49 9.3e-02 1.50 D 3,fe r,c

Figure 2 .

 2 Figure 2. Left: mass-lumping region for P 2 finite elements. Right: division of a mesh triangle in four to construct D 1,fe α,1/4 .

Figure 3 . 1

 31 Figure3. The three types of 2D meshes: "e" mesh (left), "s" mesh (middle), "r" mesh (right).

Tests with D k α , k = 1, 2

 2  present the results for the 2D versions of the Test Cases P-1, P-2, S-1 and S-2, that is: porous medium with f = 0, porous medium with f = 0, Stefan problem with f = 0, and Stefan problem with f = 0. Plots of solutions for the Stefan problems are given in Figure4(2D version of Test Case S-1, f = 0) and Figure5(2D version of Test Case S-2, f = 0).

Figure 4 ., D 1

 41 Figure 4. Approximate functions (top: u, bottom: ζ(u)) for the 2D version of Test Case S-1. From left to right: D 1,fe e , D 1,fe s , D 1,fe r .

3. 3 .

 3 Numerical tests for mass-lumped DG schemes. The mesh M being a general polytopal mesh as in [12, Definition 7.2], still using the notations in Assumption 2.19 the gradient discretisation D * = D k,dg *

Figure 6 .

 6 Figure 6. Comparison of approximate u on Test Case S-1: k = 1 (left), k = 3 (right), Finite Element (blue), dG (green), exact (red), N = 16, uniform mesh.

Appendix A .

 . Existence and uniqueness of the weak solutionTheorem A.1 (Existence and uniqueness of the weak solution). Under Assumption (1.3), there is a unique solution u to(1.2). This solution has the following regularity properties:• Λ∇ζ(u) + F ∈ H div (Ω); • if d ≤ 3 and F ∈ L p (Ω) d for some p > d, then ζ(u) ∈ C θ (Ω)for some θ ∈ (0, 1) depending only on Ω, Λ and p;• if F = 0, Ω is convex and Λ is Lipschitz-continuous, then ζ(u) ∈ H 2 (Ω).Proof. The existence of a solution is a consequence of Theorem 2.8, together with Lemma A.2 that establishes the existence of a proper sequence of gradient discretisations. To prove the uniqueness of this solution, consider u 1 and u 2 two solutions to (1.2), subtract their respective equations, and takev = ζ(u 1 ) -ζ(u 2 ) ∈ H 1 0 (Ω) as a test function to get Ω (β(u 1 ) -β(u 2 ))(ζ(u 1 ) -ζ(u 2 )) + Ω Λ∇(ζ(u 1 ) -ζ(u 2 )) • ∇(ζ(u 1 ) -ζ(u 2 )) = 0.The first term is non-negative since β and ζ are non-decreasing, and thus ∇(ζ(u 1 ) -ζ(u 2 )) = 0. This shows that that ζ(u 1 ) = ζ(u 2 ). The weak formulation (1.2) also shows that β(u 1 ) -∆ζ(u 1 ) = f + div(F ) = β(u 2 ) -∆ζ(u 2 ) in the sense of distributions on Ω; since ζ(u 1 ) = ζ(u 2 ), this yields β(u 1 ) = β(u 2 ). Hence, β(u 1 ) + ζ(u 1 ) = β(u 2 ) + ζ(u 2 ) and Hypothesis (1.3d) shows that u 1 = u 2 .

  2), use the monotonicity of µ and ρ, the relation µ (w m )|∇w| 2 = |∇ν(w m )| 2 , the coercivity of Λ and the Poincaré inequality, we write (with a b meaning a ≤ Cb with C independent of m):λ ∇ν(w m ) 2 L 2 ≤ f L 2 w m L 2 f L 2 (1 + µ(w m ) L 2 ) f L 2 ∇µ(w m ) L 2 ). (B.4) We have |∇µ(w m )| 2 = µ (w m )|∇ν(w m )| 2 ≤ |∇ν(w m )| 2and the estimate above therefore gives a bound on ν(w m ) in H 1 0 (Ω), and thus also on µ(w m ). Using a coercivity property of µ similar to that of ζ we infer bounds in L 2 (Ω) on w m and ρ(w m ). A topological degree argument, similar to the one developed in the proof of Lemma 2.6 below, then ensures the existence of at least one solution w m to (B.2).These bounds give v ∈ H 1 0 (Ω) and w ∈ L 2 (Ω) such that, up to a subsequence, µ(w m ) → v strongly in L 2 (Ω), ∇µ(w m ) → ∇v weakly in L 2 (Ω) d and w m → w weakly in L 2 (Ω). By weak/strong convergence we infer thatlim m→∞ Ω w m µ(w m ) = Ω w vand a Minty argument[START_REF] Droniou | The gradient discretisation method[END_REF] Lemma D.10] yields v = µ(w), and thus ρ(w m ) → w -µ(w) = ρ(w) weakly in L 2 (Ω). We have (ν(a) -ν(b)) 2 ≤ (b -a)(µ(b) -µ(a)) and the strong convergence of µ(w m ) in L 2 therefore shows that ν(w m ) → ν(w) in L 2 (Ω). Since (ν(w m )) m∈N is bounded in H 1 0 (Ω), this convergence also holds weakly in this space.

5 )Ω w 2 ,

 52 3) holds, we apply (B.2) with v = w m to get lim m→∞ Ω ρ(w m )w m + Ω Λ∇ν(w m ) • ∇ν(w m ) = The weak convergence of ν(w m ) in H 1 0 (Ω) ensures that lim inf m→∞ Ω Λ∇ν(w m ) • ∇ν(w m ) ≥ Ω Λ∇ν(w) • ∇ν(w). (B.6) Developing the relation Ω (ρ(w m ) -ρ(w))(w m -w) ≥ 0 and using the weak convergences w m → w and ρ(w m ) → ρ(w) in L 2 (Ω) we have lim inf m→∞ Ω ρ(w m )w m ≥ Ω ρ(w)w. (B.7) Using (B.6) and (B.7) together with (B.5) yields Ω Λ∇ν(w m ) • ∇ν(w m ) → Ω Λ∇ν(w) • ∇ν(w) and Ω ρ(w m )w m → Ω ρ(w)w. The first relation classically shows that ∇ν(w m ) → ∇ν(w) strongly in L 2 (Ω). Using the second relation and a weak/strong convergence argument on µ(w m )w m we infer that m )w m + µ(w m )w m → Ω ρ(w)w + µ(w)w = which gives the strong convergence in L 2 (Ω) of w. Remark B.3 (About the assumption F = 0). If F = 0, then an additional term Ω F •∇w m appears in the sequence (B.4), which cannot be estimated since no a priori bound is expected on w m in H 1 0 (Ω).

  [START_REF] Eymard | Error estimate for approximate solutions of a nonlinear convectiondiffusion problem[END_REF] (Rate of convergence). For all classical mass-lumping of schemes based on a mesh of size h, we have α D ,D * = O(h) (see, e.g.,[START_REF] Droniou | The gradient discretisation method[END_REF] Eqs. (8.18) and (9.46)]). Likewise, any reasonable quadrature rule is locally exact on piecewise constant functions and thus, if β, ζ are globally Lipschitzcontinuous and u, f are locally H 1 , we expect O(h) estimates on β

  .16) Assumption 2.16 only imposes a local smoothness of f and β(u), which can in particular be discontinuous across cell interfaces. It is also worthwhile noticing that, since F = 0, ζ(u) is continuous (see Theorem A.1). Hence, the values of u at one of its discontinuities must belong to a plateau of ζ; in particular, if ζ does not have any plateau, then u is globally continuous.

Remark 2.18 (u is in C(M)). Theorem A.1 and Assumption 2.16 show that ζ(u) ∈ C(Ω) and

  2.25 (Local polynomial space). Following Remark 2.20, if P k is replaced by P K in Item (3) of Assumption 2.19, then an inspection of the proof below (see in particular the polynomial (2.24)) shows that (2.18) only has to be assumed on the smaller space P K × P l . This is similar to what was noticed in [17], in a the context of mass-lumped P k Finite Elements for linear equations. Remark 2.26 (Rates of convergence). If D * is the gradient discretisation corresponding to conforming P k finite elements, we have W D L 2 h k ; see [12, Proposition 8.11 and Remark 8.12]. In this case, (2.19) yields an O

* ≡ 0 and, if ζ(u) ∈ H k+1 (Ω), ∇ D I D ζ(u) -∇ζ(u)

Table 1 .

 1 Constants and rates for Test Case R. Test case R: Regular problem, f = 0 The results provided in Table 1 show a super-convergence, for k = 1, 2, of the function and gradient reconstruction when quadrature or interpolation are accounted for in the measure of the error: the errors E Π β,I D and E Π ζ,I D appear to be O(h k+2 ), and E ∇ ζ,I D decays as h k+1 . The rate for E ∇ ζ falls to h k

	u	fe	4.6e-01 2.00 4.6e-01 2.00 4.4e-01 2.00 1.3e+00 1.00
	D 1,fe r	2.5e-01 1.89 2.5e-01 1.89 3.1e-01 1.90 1.2e+00 0.99
	D 2,fe u	8.8e-02 3.83 8.8e-02 3.83 1.4e-01 3.00 4.4e-01 2.00
	D 2,fe r	7.4e-02 3.77 7.4e-02 3.77 1.3e-01 2.98 4.2e-01 1.98
	D 3,fe u,a	1.8e-01 2.00 1.8e-01 2.00 1.5e-01 1.00 1.5e-01 1.00
	D 3,fe r,a	2.0e-01 2.01 2.0e-01 2.01 1.5e-01 1.00 1.5e-01 1.00
	D 3,fe u,b	9.4e-02 3.00 9.4e-02 3.00 2.0e-01 2.00 2.0e-01 2.00
	D 3,fe r,b	9.6e-02 2.99 9.6e-02 2.99 2.0e-01 1.99 2.0e-01 1.99
	D 3,fe u,c	6.4e-08 1.73 6.4e-08 1.73 2.0e-04 2.95 7.2e-02 3.00
	D 3,fe r,c	9.0e-08 1.80 9.0e-08 1.80 2.4e-04 2.97 7.6e-02 3.00

Table 2 .

 2 Constants and rates for Test Case P-1.

	u	fe	2.3e+02 1.68 5.6e+00 2.01 1.2e+01 2.00 3.2e+00 1.00
	D 1,fe r	3.4e+02 1.71 5.3e+00 2.00 1.2e+01 1.98 3.4e+00 1.01
	D 2,fe u	1.9e+02 1.71 1.3e+00 2.69 4.3e+00 2.45 6.9e+00 2.01
	D 2,fe r	9.5e+01 1.50 1.3e+01 3.16 1.6e+01 2.69 6.4e+00 1.99
	D 3,fe u,a	8.0e+01 1.82 4.4e-01 2.01 4.1e-01 1.03 4.0e-01 1.02
	D 3,fe r,a	9.2e+01 1.74 5.2e-01 2.03 4.2e-01 1.03 4.2e-01 1.03
	D 3,fe u,b	8.6e+01 1.74 2.8e+00 2.90 2.7e+00 1.99 2.7e+00 1.99
	D 3,fe r,b	3.0e+01 1.46 3.9e+00 3.01 2.3e+00 1.95 2.4e+00 1.96
	D 3,fe u,c	1.7e+01 1.41 1.0e+00 2.92 1.2e+00 2.42 2.7e+00 2.41
	D 3,fe r,c	3.5e+01 1.52 7.0e-02 2.17 2.2e-01 1.98 1.7e+00 2.28

Table 7 .

 7 Constants and rates for Test Case S-3.

	u	fe	3.8e+01 0.50 3.5e+01 2.01 7.7e+00 1.49 1.2e+00 0.71
	D 1,fe r	2.3e+01 0.42 4.0e-01 0.81	5.8e-01 0.82 5.7e-01 0.34
	D 2,fe u	2.2e+01 0.50 3.6e+00 2.00 1.6e+00 1.50 3.7e-01 0.51
	D 2,fe r	1.2e+01 0.42 3.1e-03 -0.29 7.2e-02 0.26 6.6e-01 0.44
	D 3,fe u,a	2.2e+01 0.50 3.3e+00 2.01	6.5e-01 1.18 3.6e-01 0.51
	D 3,fe r,a	2.9e+00 0.17 4.7e-03 -0.03 6.3e-02 0.14 3.6e-01 0.35
	D 3,fe u,b	1.8e+01 0.50 2.3e+00 2.00 1.0e+00 1.50 3.6e-01 0.50
	D 3,fe r,b	5.4e+00 0.22 4.3e-01 0.95	2.2e-01 0.46 7.0e-01 0.50
	D 3,fe u,c	1.5e+01 0.50 8.8e-01 2.00	5.7e-01 1.50 3.5e-01 0.50
	D 3,fe r,c	5.5e+00 0.26 9.1e-03 -0.05 8.3e-02 0.33 6.4e-01 0.48
			x i
			U i

Table 8 .

 8 Constants and rates for the 2D version of Test Case P-1.

		9.1e-03 2.06 2.9e-02 2.95 5.0e-03 2.03 3.0e-03 2.59 2.0e-02 2.02 8.1e-03 2.50
	E Π ζ,ID	3.3e-04 2.07 2.8e-03 3.53 1.5e-04 2.04 2.0e-03 4.04 1.5e-03 2.03 1.2e-03 3.01
	E ∇ ζ,ID	1.4e-03 1.51 7.7e-03 2.52 7.5e-04 2.04 6.9e-03 3.02 2.0e-03 1.02 2.8e-03 2.02

Table 9 .

 9 Constants and rates for the 2D version of Test Case P-2.

	1,fe α,1/4 and D 2,fe α : comparison between degree 1 and degree 2

Table 10 .

 10 Constants and rates for the 2D version of Test Case S-1.

	1.2	1.2	1.2
	1.0	1.0	1.0
	0.8	0.8	0.8
	0.6	0.6	0.6
	0.4	0.4	0.4
	0.2	0.2	0.2
	0.0	0.0	0.0
	0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0	0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0	0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
	0.008	0.008	0.008
	0.006	0.006	0.006
	0.004	0.004	0.004
	0.002	0.002	0.002
	0.000	0.000	0.000
	0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0	0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0	0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Table 11 .

 11 Constants and rates for the 2D version of Test Case S-2.

		2.5e-01 0.48 4.0e-01 0.58 8.1e-02 0.35 5.0e-01 0.68 1.4e-01 0.37 7.7e-01 0.72
	E Π ζ,ID	2.4e-02 2.10 1.6e-02 2.10 2.4e-02 2.24 3.2e-02 2.23 3.8e-02 2.06 1.5e-02 1.86
	E ∇ ζ,ID	7.5e-02 1.57 9.5e-02 1.51 7.2e-02 1.71 9.5e-02 1.52 5.0e-02 1.02 9.8e-02 1.49

Table 12 .

 12 Constants and rates for the comparison first/second order with the same number of degrees of freedom.

	2,fe r	and D 1,fe r,1/4 . This makes for a fair comparison since

these two schemes have the same number of unknowns. For each errors E = E Π β,I D , E = E Π ζ,I D and E = E ∇ ζ,I D , letting E k be the error corresponding to D 1,fe r,1/4 if k = 1, or to D 2,fe r

Table 13 .

 13 the unique polynomial in P k that takes the values v i at x i for all i ∈ I K . Constants and rates for E ∇ ζ,I D with DG applied to some 1D test cases.

		Test Case R	Test Case P-1 Test Case P-2 Test Case S-2
	GD	C	α	C	α	C	α	C	α
	D 3,dg u,a	1.5e-01 1.01 4.2e-01 1.03 1.5e-01 1.01 9.0e-02 1.01
	D 3,dg r,a	1.7e-01 1.02 4.2e-01 1.03 1.5e-01 1.01 8.5e-02 1.00
	D 3,dg u,b	2.2e-01 2.00 2.9e+00 1.98 2.7e-01 2.00 8.2e-02 1.50
	D 3,dg r,b	2.2e-01 1.99 2.9e+00 1.97 2.8e-01 2.00 7.4e-02 1.57
	D 3,dg u,c	2.3e-02 3.25 1.4e+00 2.39 3.9e-02 3.42 5.4e-02 1.49
	D 3,dg r,c	1.1e-02 3.01 1.0e+00 2.32 1.9e-02 3.08 5.8e-02 1.58

stated in the theorem. Finally, the H 2 regularity property is a straightforward consequence of the optimal elliptic regularity on convex domains for Lipschitz-continuous diffusion tensor.

Lemma A.2 (Existence of suitable sequences of GDs). Under Assumption (1.3a), there exists a sequence (D m ) m∈N = (X D m ,0 , Π D m , ∇ D m , Q D m ) m∈N of gradient discretisations, with piecewise constant reconstructions, that satisfy the coercivity, consistency, limit-conformity and compactness properties stated in Theorem 2.8.

Proof. Let ( M m ) be a sequence of conformal simplicial meshes of R d (see, e.g., [START_REF] Droniou | The gradient discretisation method[END_REF]Definition 7.4]), such that max T ∈ Mm diam(T ) → 0 and (M m ) m∈N is regular in the sense that the ratio of the diameter of T ∈ M m over the largest ball inside T is bounded uniformly with respect to T and m. We let M m = {T ∈ M m : T ⊂ Ω} and define the polyhedral set Ω m ⊂ Ω as the interior of ∪ T ∈Mm T m .

The gradient discretiation

) is defined as the mass-lumped conforming P 1 gradient discretisation on the mesh M m of Ω m [START_REF] Droniou | The gradient discretisation method[END_REF]Section 8.4], with extensions to Ω by 0 outside Ω m , and no quadrature rule. Letting V m be the set of vertices of M m , we therefore set

and

is on Ω m the gradient of the conforming P 1 reconstruction from the vertex values (v i ) i∈Vm , and

Since the functions and gradient reconstructions are extended by 0 outside Ω m , C D m and W D m can be computed using norms and integrals over Ω m . The properties of mass-lumped P 1 GDs on Ω m (see [START_REF] Droniou | The gradient discretisation method[END_REF]Theorem 8.17]) then show that (D m ) m∈N is coercive, limit-conforming and compact. It remains to analyse the consistency of (D m ) m∈N .

As seen in [12, Lemma 2.16], the consistency follows if we prove that S D m (ϕ) → 0 when ϕ ∈ C ∞ c (Ω). In that case, for m large enough, ϕ ∈ C ∞ c (Ω m ) and the norms in S D m (ϕ) can be restricted to Ω m . The estimate in [START_REF] Droniou | The gradient discretisation method[END_REF]Remark 8.18] then show that S D m (ϕ) ≤ C ϕ max T ∈Mm diam(T ) with C ϕ not depending on m. This shows that S D m (ϕ) → 0 as m → ∞, as required.

Appendix B. Conforming scheme

Throughout this section, we assume that F = 0. Using Assumptions (1.3b), (1.3c) and (1.3d), we see that β + ζ : R → R is bijective and we can therefore set µ(t) = ζ((β + ζ) -1 (t)) and ρ(t) := t -µ(t) = β((β + ζ) -1 (t). These functions are non-decreasing and 1-Lipschitz continuous and, setting w = (β + ζ)(u), we see that (1.2) is equivalent to: find w ∈ L 2 (Ω) such that µ(w) ∈ H 1 0 (Ω) and

Given a family (V m ) m∈N of finite dimensional subspaces of H 1 0 (Ω), conforming schemes for (B.1) are written: find

Introducing the function ν : R → R defined by ν(s) = s 0 µ (r)dr, we can then state the following convergence theorem.