Discontinuous Skeletal Gradient Discretisation Methods on polytopal meshes
DDiscontinuous Skeletal Gradient Discretisation Methods on polytopalmeshes (cid:73)
Daniele A. Di Pietro a , J´erˆome Droniou b , Gianmarco Manzini c a Institut Montpelli´erain Alexander Grothendieck, CNRS, Univ. Montpellier (France) b School of Mathematical Sciences, Monash University, Melbourne (Australia) c T-5 Applied Mathematics and Plasma Physics Group, Los Alamos National Laboratory, Los Alamos, New Mexico (USA)
Abstract
In this work we develop arbitrary-order Discontinuous Skeletal Gradient Discretisations (DSGD) on gen-eral polytopal meshes. Discontinuous Skeletal refers to the fact that the globally coupled unknowns arebroken polynomials on the mesh skeleton. The key ingredient is a high-order gradient reconstructioncomposed of two terms: (i) a consistent contribution obtained mimicking an integration by parts for-mula inside each element and (ii) a stabilising term for which sufficient design conditions are provided.An example of stabilisation that satisfies the design conditions is proposed based on a local lifting ofhigh-order residuals on a Raviart–Thomas–N´ed´elec subspace. We prove that the novel DSGDs satisfycoercivity, consistency, limit-conformity, and compactness requirements that ensure convergence for a va-riety of elliptic and parabolic problems. Links with Hybrid High-Order, non-conforming Mimetic FiniteDifference and non-conforming Virtual Element methods are also studied. Numerical examples completethe exposition.
Keywords:
Gradient discretisation methods, Gradient Schemes, high-order Mimetic Finite Differencemethods, Hybrid High-Order methods, Virtual Element methods, non-linear problems
1. Introduction
The numerical resolution of (linear or non-linear) partial differential equations (PDEs) is nowadaysubiquitous in the engineering practice. In this context, the design of convergent numerical schemesis a very active research topic. The Gradient Discretisation Method (GDM) is a recently introducedframework which identifies key design properties to obtain convergent schemes for a variety of linear andnon-linear elliptic and parabolic problems. Several models of current use in fluid mechanics fall into thelatter categories including, e.g., porous media flows governed by Darcy’s law, phase change problemsgoverned by the Stefan problem [3636], as well as simplified models of the viscous terms in power-lawfluids corresponding the Leray–Lions elliptic operators. The latter also appear in the modelling of glaciermotion [3838], of incompressible turbulent flows in porous media [2626], and in airfoil design [3737].A Gradient Discretisation (GD) is defined by a finite-dimensional space encoding the discrete un-knowns, as well as two linear operators acting on the latter, and corresponding to reconstructions ofscalar functions and of their gradient. For a given PDE problem, convergent GDs are characterisedby four properties, which can also serve as guidelines for the design of new schemes: coercivity , which (cid:73)
The first author was supported by the ANR grant HHOMM (ANR-15-CE40-0005). The second author was supportedby the ARC Discovery Projects funding scheme (project number DP170100605). The third author was funded by the Lab-oratory Directed Research and Development program, under the auspices of the National Nuclear Security Administrationof the U.S. Department of Energy by Los Alamos National Laboratory, operated by Los Alamos National Security LLCunder contract DE-AC52-06NA25396. a r X i v : . [ m a t h . NA ] D ec orresponds to a discrete Poincar´e inequality; GD-consistency , which expresses the ability of the scalarand gradient reconstructions to approximate functions in the space where the continuous problem is set; limit-conformity , linking the two reconstructions through an approximate integration by parts formula; compactness , corresponding to a discrete counterpart of the Rellich theorem.In the recent monograph [2828], several classical discretisation methods have been interpreted in theGDM framework. These include: arbitrary-order conforming, nonconforming, and mixed Finite Elements(FE) on standard meshes; arbitrary-order discontinuous Galerkin (DG) schemes in their SIPG form [11](see, in particular, [3535] on this point); various lowest-order Finite Volume methods on specific grids;lowest-order methods belonging to the Hybrid Mixed Mimetic family (see the unified presentation in [2929]of the methods originally proposed in [88, 2727, 3434]) as well as nodal Mimetic Finite Differences (MFD) [99]on arbitrary polyhedral meshes; see also [44].In this paper we present an important addition to the GDM framework: arbitrary-order DiscontinuousSkeletal (DS) methods [1818], characterised by globally coupled unknowns that are broken polynomials onthe mesh skeleton. Specifically, the primary source of inspiration are the recently introduced Hybrid High-Order (HHO) methods for linear [2222, 2020] and non-linear [1616, 1717] diffusion problems, and the high-ordernon-conforming MFD (ncMFD) method of [4141]; see also [22] for an interpretation in the Virtual Elementframework and [33] for an introduction to the latter. We also cite here the Hybridizable DiscontinuousGalerkin methods of [1414], whose link with the former methods has been studied in [1313]; see also [66]for a unified formulation. Like DG methods, DS methods support arbitrary approximation orders ongeneral polytopal meshes. DS methods are, in addition, amenable to static condensation for linear(ised)problems, which can significantly reduce the number of unknowns in some configurations. They alsohave better data locality, which can ease parallel implementations. Moreover, lowest-order versions areoften available that can be easily fitted into traditional Finite Volume simulators. Finally, unlike DGmethods, DS methods admit a Fortin operator in general meshes, a crucial property in the context ofincompressible or quasi-incompressible problems in solid- and fluid-mechanics; see, e.g., [2020, 2323].Let a polynomial degree k ě k on the mesh skeleton, as well as locally coupled element-basedunknowns that correspond to broken polynomials of total degree up to l P t k ´ , k, k ` u on the meshitself. The reconstruction of scalar functions is defined in a straightforward manner through the latterif l ě
0, or by a suitable combination of face-based unknowns if l “ ´
1. The gradient reconstruction,on the other hand, requires a more careful design. The seminal ideas to devise high-order gradientreconstructions on general meshes are already present, among others, in HHO methods (see, e.g., [2222,Eq. (13)] and [1616, Eq. (4.3)]) as well as in ncMFD methods (see [4141, Eq. (21)]). These gradientreconstructions, however, are not suitable to define a convergent DSGD because they fail to satisfy thecoercivity requirement. In addition, when considering non-linear problems, the codomain of the gradientreconstruction has to be carefully selected in order for the GD-consistency requirement to be satisfiedwith optimal scaling in the meshsize for k ě L -orthogonality to vector-valued polynomials of degree up to k . Whenconsidering problems posed in a non-Hilbertian setting, an additional condition is added stipulating thatthe stabilisation is built on a piecewise polynomial space. An example of stabilisation term that meetsall of the above requirement is proposed based on a Raviart–Thomas–N´ed´elec space on a submesh.The rest of the paper is organised as follows. In Section 22 we recall the basics of the GDM and give afew examples of linear and non-linear problems for which GDs are convergent under the coercivity, GD-2onsistency, limit-conformity, and compactness properties discussed above. The construction of arbitrary-order DSGD is presented in Section 33, the main results are stated in Section 3.53.5, and numerical examplesare collected in Section 3.73.7. The links with HHO, ncMFD, and ncVEM schemes are studied in detailin Section 44. Appendix AAppendix A contains the proofs of the main results. The material is organised so thatmultiple levels of reading are possible: readers mainly interested in the numerical recipe and results canprimarily focus on Sections 22–33; readers also interested in the relations with other polytopal methods canconsult Section 44.
2. The Gradient Discretisation Method
We give here a brief presentation of the Gradient Discretisation Method (GDM) in the context ofhomogeneous Dirichlet boundary conditions, and we refer to the monograph [2828] for more details andother boundary conditions.
Let Ω be a bounded polytopal domain in R d , where d ě W ,p p Ω q , where p P p , `8q denotes a Sobolevexponent which we assume fixed in what follows.A Gradient Discretisation (GD) is a triplet D “ p X D , , Π D , ∇ D q where:(i) X D , is a finite dimensional vector space on R encoding the discrete unknowns , and accounting forthe homogeneous Dirichlet boundary condition;(ii) Π D : X D , Ñ L p p Ω q is a linear mapping that reconstructs scalar functions in L p p Ω q from thediscrete unknowns in X D , ;(iii) ∇ D : X D , Ñ L p p Ω q d is a linear mapping that reconstructs the “gradient” of scalar functions in L p p Ω q d from the unknowns in X D , . This reconstruction must be defined such that } ∇ D ¨} L p p Ω q d isa norm on X D , .In a nutshell, the GDM consists in selecting a GD and in replacing, in the weak formulation ofthe PDE, the continuous space and operators by the discrete ones provided by the GD. The schemethus obtained is called a Gradient Scheme (GS). To illustrate this procedure, consider the simple linearproblem: Find u : Ω Ñ R such that ´ ∇ ¨p Λ ∇ u q “ f in Ω ,u “ B Ω , (1)with diffusion tensor Λ bounded and uniformly coercive, and source term f P L p Ω q . The weak formu-lation of (11) is Find u P H p Ω q such that, for all v P H p Ω q , ż Ω Λ ∇ u ¨ ∇ v “ ż Ω f v. (2)Given a gradient discretisation D , the gradient scheme for (22) is thenFind u D P X D , such that, for all v D P X D , , ż Ω Λ ∇ D u D ¨ ∇ D v D “ ż Ω f Π D v D . (3)The same procedure applies to non-linear problems. Consider, e.g., the following generalisation of (11)that corresponds to Leray–Lions operators: Find u : Ω Ñ R such that ´ ∇ ¨ σ p x , u, ∇ u q “ f in Ω ,u “ B Ω , (4)3here the flux function σ : Ω ˆ R ˆ R d Ñ R d satisfies the requirements detailed in [2828, Eq. (2.85)]. Aparadigmatic example of this class of problems is the p -Laplace equation which, for a fixed p P p , `8q ,corresponds to the flux function σ p x , u, ∇ u q “ | ∇ u | p ´ ∇ u. (5)Assuming f P L p p Ω q with p – pp ´ , Problem (44) admits the following weak formulation:Find u P W ,p p Ω q such that, for all v P W ,p p Ω q , ż Ω σ p x , u, ∇ u q¨ ∇ v “ ż Ω f v . (6)Given a gradient discretisation D , the gradient scheme for (66) is thenFind u D P X D , such that, for all v D P X D , , ż Ω σ p x , Π D u, ∇ D u D q¨ ∇ D v D “ ż Ω f Π D v D . (7) The accuracy and convergence of GSs, for linear and non-linear problems, can be assessed by a fewproperties of the underlying GDs. In many situations, and in all cases considered in this paper, GDsare obtained starting from a mesh of the domain. We consider here polytopal meshes corresponding tocouples M h – p T h , F h q , where T h is a finite collection of polytopal elements T of maximum diameterequal to h ą
0, while F h is a finite collection of hyperplanar faces F . It is assumed henceforth that themesh M h matches the weak geometrical requirements detailed in [2828, Definition 7.2]; see also [2525, Section2]. Our focus is on the so-called h -convergence analysis, where we consider a sequence of refined meshes p M h q h P H whose sizes are collected in a countable set H Ă R `˚ having 0 as its unique accumulation point.We further assume that the polytopal mesh sequences that we deal with are regular in the sense of [2525,Definition 3], and we denote by (cid:37) ą p D h q h P H “ p X D h , , Π D h , ∇ D h q h P H of GDsthat lead to gradient schemes that converge, for both linear and non-linear problems: (GD1) Coercivity.
Consider, for all h P H , the norm of the linear mapping Π D h defined by: C D h – max v P X D h, zt u } Π D h v } L p p Ω q } ∇ D h v } L p p Ω q d . Then, there exists a real number C P ą C D h ď C P for all h P H . (GD2) GD-Consistency.
For all h P H , let S D h : W ,p p Ω q Ñ r , `8q be such that S D h p φ q – min v P X D h, ` } Π D h v ´ φ } L p p Ω q ` } ∇ D h v ´ ∇ φ } L p p Ω q d ˘ @ φ P W ,p p Ω q . Then, it holds that lim h Ñ S D h p φ q “ @ φ P W ,p p Ω q . (8) (GD3) Limit-conformity.
Let p – pp ´ denote the dual exponent of p , and set W p p div; Ω q – t ψ P L p p Ω q d : ∇ ¨ ψ P L p p Ω qu . For all h P H , let W D h : W p p div; Ω q Ñ r , `8q be such that, for all ψ P W p p div; Ω q , W D h p ψ q – sup v P X D h, zt u } ∇ D h v } L p p Ω q d ˇˇˇˇż Ω ´ ∇ D h v p x q¨ ψ p x q ` Π D h v p x q ∇ ¨ ψ p x q ¯ d x ˇˇˇˇ . Then, it holds that lim h Ñ W D h p ψ q “ @ ψ P W p p div; Ω q . (9)4 GD4)
Compactness.
For any v h P X D h , such that p} ∇ D h v h } L p p Ω q d q h P H is bounded, the sequence p Π D h v h q h P H is relatively compact in L p p Ω q .A few comments are of order. Property (GD1)(GD1) is linked to the stability of the method, and stipulatesthat the L p -norm of the reconstruction of scalar functions is uniformly controlled by the L p -norm ofthe reconstruction of their gradient. It readily implies the uniform Poincar´e inequality } Π D h v h } L p p Ω q ď C P } ∇ D h v h } L p p Ω q d valid for all h P H and all v h P X D h , .Properties (GD2)(GD2) and (GD3)(GD3) are linked to the consistency of the method. More specifically, property (GD2)(GD2) states that the reconstructions Π D h of scalar functions and ∇ D h of their gradients are able toapproximate functions that lie in the space W ,p p Ω q where the continuous problem is set. In the context ofthe FE convergence analysis, this property is an attribute of the underlying discrete space, and is usuallycalled approximability ; see, e.g., [3333, Definition 2.14]. Property (GD3)(GD3) , on the other hand, establishesa link between Π D h and ∇ D h in the form of a discrete integration by parts formula. Its counterpart inthe context of the FE convergence analysis for linear problems is asymptotic consistency ; see, e.g., [3333,Definition 2.15]. Notice, however, that the formulation in (GD3)(GD3) is in a sense more general, as it isnot linked to a specific underlying problem and is in particular readily applicable to non-linear problems(whereas [3333, Definition 2.15] is restricted to linear problems).Finally, property (GD4)(GD4) is a discrete Rellich compactness theorem, and can be regarded as the keyingredient to obtain strong convergence results by compactness techniques. Remark . Either one of (GD3)(GD3) or (GD4)(GD4) imply (GD1)(GD1) , see [2828, Lemmas 2.7 and 2.11]. The coercivity is however kept as a separate property to highlightits importance.The above properties are sufficient to carry out a convergence analysis, either by error estimates(when the model is amenable to these) or by compactness, for a variety of linear and non-linear elliptic orparabolic models. An example of such convergence results for gradient discretisations of the Leray–Lionsproblem (66) is provided by Theorems 22 and 33 below; see [2828] for a comprehensive collection of convergenceresults for various linear and non-linear elliptic and parabolic problems. Theorem 2 (Convergence) . We assume that σ satisfies the classical properties of Leray–Lions operators(see [2828, Eqs. (2.85) and (2.87)]). Let p D h q h P H denote a sequence of GDs satisfying (GD1)(GD1) – (GD4)(GD4) .Then, for all h P H , there exists at least one u D h P X D h , solution to (77) and, along a subsequence as h Ñ , (i) Π D h u D h converges strongly in L p p Ω q to a solution u of (66) ; (ii) ∇ D h u D h converges stronglyin L p p Ω q d to ∇ u .Proof. This is a special case of [2828, Theorem 2.45].
Theorem 3 (Error estimates) . Let D be a gradient discretisation, and let σ be the Leray–Lions operatorcorresponding to the p -Laplace equation (see (55) ). Then there exists a unique u D solution to (77) and, if u is the solution to (66) , then there exists C depending only on p , f and an upper bound of C D such that • If ă p ď , } u ´ Π D u D } L p p Ω q ` } ∇ u ´ ∇ D u D } L p p Ω q ď C “ S D p u q ` S D p u q p ´ ` W D p σ p ∇ u qq ‰ . • If ď p , } u ´ Π D u D } L p p Ω q ` } ∇ u ´ ∇ D u D } L p p Ω q ď C ” S D p u q ` S D p u q p ´ ` W D p σ p ∇ u qq p ´ ı . Proof.
These error estimates are simplified forms of the ones in [2828, Theorem 2.39].
3. Discontinuous Skeletal Gradient Discretisations
In this section, we construct a family of Discontinuous Skeletal Gradient Discretisations (DSGD). Thenotation is closely inspired by HHO methods; see, e.g., [2525].5 .1. Local polynomial spaces and projectors
Local polynomial spaces on mesh elements or faces and projectors thereon play a crucial role in thedesign and analysis of DSGD methods.For any X Ă Ω, we denote by p¨ , ¨q X the standard L p X q - or L p X q d -products. This notation isused in place of integrals when dealing with quantities that are inherently L -based. Let now X be amesh element or face. For an integer (cid:96) ě P (cid:96) p X q denotes the space spanned by the restriction to X ofscalar-valued, d -variate (if X is a mesh element) or p d ´ q -variate (if X is a face) polynomials of totaldegree (cid:96) or less, and conventionally set P ´ p X q – t u .Let again X denote a mesh element or face. The L -orthogonal projector π ,(cid:96)X : L p X q Ñ P (cid:96) p X q isdefined as follows: For all v P L p X q , π ,(cid:96)X is the unique polynomial in P (cid:96) p X q such that p π ,(cid:96)X v ´ v, w q X “ @ w P P (cid:96) p X q . (10)In the vector case, the L -projector is defined component-wise and denoted by π ,(cid:96)X .For any mesh element T P T h , we also define the elliptic projector π ,(cid:96)T : W , p T q Ñ P (cid:96) p T q as follows:For all v P W , p T q , π ,(cid:96)T v is the unique polynomial in P (cid:96) p T q that satisfies p ∇ p π ,(cid:96)T v ´ v q , ∇ w q T “ w P P (cid:96) p T q and p π ,(cid:96)T v ´ v, q T “ π ,(cid:96)T and π ,(cid:96)T have optimal approximation properties in P (cid:96) p T q (see Theorem 1.1, Theorem 1.2, and Lemma 3.1 in [1717]): For any α P t , u and s P t α, . . . , (cid:96) ` u ,there exists a real number C ą h , but possibly depending only on d , p , (cid:37) , (cid:96) , α , and s ,such that, for all T P T h , and all v P W s,p p T q , | v ´ π α,(cid:96)T v | W r,p p T q ď Ch s ´ rT | v | W s,p p T q @ r P t , . . . , s u , (11a)and, if s ě h p T | v ´ π α,(cid:96)T v | W r,p p F T q ď Ch s ´ rT | v | W s,p p T q @ r P t , . . . , s ´ u , (11b)where W r,p p F T q – (cid:32) v P L p pB T q : v | F P W r,p p F q for all F P F T ( and h T denotes the diameter of theelement T . We continue our discussion with a crucial remark concerning the computation of the L -orthogonalprojection of the gradient from L -orthogonal projections of a scalar function and its traces. This remarkwill inspire the choice of the discrete unknowns as well as the definition of the gradient reconstructionin DSGD methods. In what follows, we work on a fixed mesh element T P T h , denote by F T the set ofmesh faces that lie on the boundary of T and, for all F P F T , by n T F the normal vector to F pointingout of T .Consider a function v P W , p T q . We note the following integration by parts formula, valid for all φ P C p T q d : p ∇ v, φ q T “ ´p v, ∇ ¨ φ q T ` ÿ F P F T p v, φ ¨ n T F q F . (12)Let now an integer k ě φ P P k p T q d , we obtain p π ,kT ∇ v, φ q T “ ´p π ,k ´ T v, ∇ ¨ φ q T ` ÿ F P F T p π ,kF v, φ ¨ n T F q F , (13)where we have used (1010) to insert π ,kT into the left-hand side, and π ,k ´ T and π ,kF into the right-handside after observing that ∇ ¨ φ P P k ´ p T q and, since we are considering planar faces, φ | F ¨ n T F P P k p F q forall F P F T . The relation (1313) shows that computing the L -orthogonal projection of ∇ v on P k p T q d doesnot require a full knowledge of the function v . All that is required is6i) π ,k ´ T v , the L -projection of v on P k ´ p T q . Other possible choices are π ,kT v or π ,k ` T v (in fact,any polynomial degree larger than or equal to k ´ F P F T , π ,kF v , the L -projection on P k p F q of the trace of v on F . Inspired by the previous remark, for two given integers k ě l P t k ´ , k, k ` u we consider thefollowing set of discrete unknowns: U k,lh – ˜ ą T P T h P l p T q ¸ ˆ ˜ ą F P F h P k p F q ¸ . The choice l “ k ´ l “ k to the Hybrid High-Ordermethod of [2222], and the choice l “ k ` k “ l “ k ´
1, element-based unknowns are not present.For a generic element of U k,lh , we use the standard HHO underlined notation v h “ pp v T q T P T h , p v F q F P F h q , and we define the interpolator I k,lh : W , p Ω q Ñ U k,lh such that, for all v P W , p Ω q , I k,lh v – ´ p π ,lT v q T P T h , p π ,kF v | F q F P F h ¯ . To account for Dirichlet boundary conditions strongly, we introduce the subspace U k,lh, – ! v h P U k,lh : v F ” F P F b h ) , where F b h is the set collecting the mesh faces that lie on the boundary of Ω.The restrictions of U k,lh , I k,lh and v h P U k,lh to a generic mesh element T P T h are denoted by U k,lT , I k,lT ,and v T , respectively. That is, U k,lT – t v T “ p v T , p v F q F P F T q : v T P P l p T q , v F P P k p F q @ F P F T u and, for all v P W , p T q , I k,lT v – p π ,lT v, p π ,kF v | F q F P F T q . Moreover, we adopt the convention that, for all T P T h , v T – ÿ F P F T ω T F v F if l ă , (14)where, following [4141, Appendix A], the weights t ω T F u F P F T are defined in such a way that ř F P F T ω T F p q, q F “p q, q T for all q P P p T q (this condition is required in the above reference to obtain L -superconvergence,not treated in this work). For all v h P U k,lh , we also define the broken polynomial field v h such that v h | T – v T @ T P T h . (15)The space of discrete unknowns and the reconstruction of the scalar variable for a DSDG are givenby, respectively, X D h , – U k,lh, and Π D h v h – v h for all v h P U k,lh, . (16) To complete the definition of a DSGD, it remains to identify a reconstruction of the gradient, whichmakes the object of this section. 7 .4.1. A consistent and limit-conforming high-order gradient
Let a mesh element T P T h be fixed. Taking inspiration from the integration by parts formula (1313),we define the gradient reconstruction G kT : U k,lT Ñ P k p T q d such that, for any v T “ p v T , p v F q F P F T q P U k,lT , G kT v T satisfies, for all φ P P k p T q d , p G kT v T , φ q T “ ´p v T , ∇ ¨ φ q T ` ÿ F P F T p v F , φ ¨ n T F q F . (17)By construction, it holds for all v P W , p T q , G kT I k,lT v “ π ,kT p ∇ v q . (18)Recalling the estimates (1111) on π ,kT , this implies that G kT I k,lT v optimally approximates ∇ v in P k p T q d when v is smooth enough.A reconstruction of the gradient that meets the consistency requirement expressed by (GD2)(GD2) can beobtained at this point letting ∇ D h be such that, for all v h P U k,lh, , ∇ D h v h “ G kh v h , (19)where G kh : U k,lh Ñ P k p T h q d is the global consistent gradient reconstruction operator obtained patchingthe local reconstructions: For all v h P U k,lh , p G kh v h q | T – G kT v T @ T P T h . (20)However, for general element shapes, the L p -norm of this gradient reconstruction is not a norm onthe space X D h , “ U k,lh, , hence the coercivity requirement expressed by (GD1)(GD1) cannot be met. Thisinitial choice of reconstructed gradient therefore has to be stabilised by accounting for jumps betweenelement and face unknowns. These jumps can be controlled in turn via a discrete counterpart of the W ,p -seminorm, which gives us an emulated Sobolev structure on U k,lh . Remark P finite elements) . If T is a d -simplex (i.e., a triangle if d “
2, a tetrahedronif d “
3, etc.) and we take k “ l “ ´
1, the gradient reconstruction G kT v T defined by (1717) coincideswith the gradient of the non-conforming P function ϕ such that | F | ´ ş F ϕ “ v F for all F P F T . Inthis case, the L p -norm of the global gradient given by (1919) defines a norm on the space of discreteunknowns, and stabilisation is not needed. To recover the usual non-conforming P scheme (often calledthe Crouzeix–Raviart scheme, although historically this name refers to the usage of non-conforming P – P discretisations for the velocity–pressure unknowns in Stokes and Navier–Stokes equations [1515]), (1616)has to be modified setting Π D h v h | T – r T v T for all T P T h and all v T P U , ´ T , where r T is the high-orderreconstruction of scalar function defined in the following section. W ,p -seminorm Let v T P U k,lT . Recalling the convention (1414), v T defines a reconstruction of scalar functions inside T of degree max p , l q . However, taking again inspiration from the integration by parts formula (1313), thistime specialised to φ “ ∇ w with w P P k ` p T q , one can define a higher-order reconstruction r k ` T : U k,lT Ñ P k ` p T q such that, for all v T P U k,lT , r k ` T v T satisfies, for all w P P k ` p T q , p ∇ r k ` T v T , ∇ w q T “ ´p v T , (cid:52) w q T ` ÿ F P F T p v F , ∇ w ¨ n T F q F . (21a)Equation (21a21a) defines r k ` T v T up to an additive constant, which we fix by imposing p r k ` T v T ´ v T , q T “ . (21b)8 emark k ` T ˝ I k,lT ) . When l ě
0, following the reasoning of [2222,Lemma 3], it can be proved that r k ` T ˝ I k,lT “ π ,k ` T , and optimal approximation properties in P k ` p T q follow from (1111) with α “ (cid:96) “ k `
1. The case l ă
0, on the other hand, can only occur when k “
0. Owing to the specific choice for the reconstruction of a (constant) element value in (1414), optimalapproximation properties analogous to (1111) with α “ (cid:96) “ U k,lh , for all T P T h we introduce the difference operators δ lT : U k,lT Ñ P l p T q and, for all F P F T , δ kT F : U k,lT Ñ P k p F q such that, for all v T P U k,lT , δ lT v T – π ,lT p r k ` T v T ´ v T q , δ kT F v T – π ,kF p r k ` T v T ´ v F q @ F P F T . (22)The role of these difference operators in the context of HHO methods has been highlighted in [2525,Section 3.1.4]. We also note here the following relation: p δ lT v T , p δ kT F v T q F P F T q “ I k,lT r k ` T v T ´ v T , (23)which will be exploited in Section 4.54.5 and Lemma 2121 below.The discrete W ,p -seminorm is defined setting } v h } p ,p,h – ÿ T P T h } v T } p ,p,T , where, for all T P T h , the local seminorm is such that, denoting by h F the diameter of the face F , } v T } p ,p,T – } G kT v T } pL p p T q d ` | v T | pp, B T , | v T | pp, B T – ÿ F P F T h ´ pF }p δ kT F ´ δ lT q v T } pL p p F q . (24)As a result of Proposition 2020 below, }¨} ,p,h is a norm on the subspace U k,lh, . We can now describe the general form of the gradient ∇ D h : U k,lh Ñ L p Ω q d , built inside each meshelement from the consistent and limit-conforming part G kT and a stabilising contribution: p ∇ D h v h q | T “ G T v T – G kT v T ` S T v T @ v h P U k,lh , @ T P T h , (25)where S T : U k,lT Ñ L p T q d satisfies the following design conditions: (S1) L -stability and boundedness. For all T P T h and all v T P U k,lT , it holds that } S T v T } L p T q d » | v T | , B T , (26)where a » b means Ca ď b ď C ´ a with real number C ą h and of T , but possiblydepending on d and on discretisation parameters including (cid:37) , k , and l . (S2) Orthogonality.
For all v T P U k,lT and all φ P P k p T q d , it holds p S T v T , φ q T “ . (27) (S3) Image. If p ‰
2, there exists k S P N independent of h and of T such that the image of S T iscontained in P k S p P T q d , the space of vector-valued broken polynomials of total degree up to k S on aregular polytopal partition P T of T . Here, regular means that, for all P P P T , denoting by r P and h P the inradius and diameter of P , respectively, it holds that (cid:37)h P ď r P , (cid:37)h T ď h P . (28)9 emark L -based stabilising contribution) . Property (S2)(S2) , which is crucial to ensure the stabilisingproperties of S T , requires to work with an inner product space. In our case, a natural choice is L p T q d .The role of orthogonality properties analogous to (S2)(S2) has been previously recognised in the context ofspecific stabilised method, see for example [77, Proposition 7] for the lowest-order Compatible Discretisa-tion Operator methods, [2828, Theorems 13.7 and 14.5] for the Hybrid Mimetic Mixed methods and thenodal Mimetic Mixed Methods, and [3232, Section 4.2] for numerical methods for elasticity models. Remark L p -stability of S T ) . Property (S3)(S3) is required to extend the stability properties expressed by (S1)(S1) to L p ; see the proof of point (i) in Proposition 1919 for further details.The above construction of a DSGD is summarised in the following Definition 8 (Discontinuous Skeletal Gradient Discretisation) . Given a polytopal mesh M h , a Discon-tinuous Skeletal Gradient Discretisation (DSGD) is given by D h “ p X D h , , Π D h , ∇ D h q where X D h , andΠ D h are defined by (1616), and ∇ D h is given by (2525) with a family of stabilisations t S T : T P M h u satisfying properties (S1)(S1) – (S3)(S3) . The construction detailed above yields a GD that meets properties (GD1)(GD1) – (GD4)(GD4) identified inSection 22, as summarised in the following Theorem 9 (Properties of DSGD) . If p M h q h P H is a regular sequence of polytopal meshes, then thesequence of the corresponding DSGDs p D h q h P H given by Definition 88 satisfies properties (GD1)(GD1) – (GD4)(GD4) .Proof. See Appendix A.3.1Appendix A.3.1.Since we are dealing with arbitrary-order methods, given the error estimates in Theorem 33, a relevantpoint consists in estimating the convergence rates of the quantities S D h p φ q (see (GD2)(GD2) ) and W D p ψ q (see (GD3)(GD3) ) when their arguments exhibit further regularity. This makes the object of the following Proposition 10 (Estimates on S D and W D ) . Let M h be a polytopal mesh and D h be a DSGD as inDefinition 88. Then, denoting by a À b the inequality a ď Cb with real number C ą not depending on h , but possibly depending on d , p , (cid:37) , k , l , and k S , it holds with l ` – max p l, q , @ φ P W ,p p Ω q X W l ` ` ,p p T h q , } Π D h I k,lh φ ´ φ } L p p Ω q À h l ` ` | φ | W l `` ,p p T h q , (29a) @ φ P W ,p p Ω q X W k ` ,p p T h q , } ∇ D h I k,lh φ ´ ∇ φ } L p p Ω q d À h k ` | φ | W k ` ,p p T h q . (29b) As a consequence, @ φ P W ,p p Ω q X W min p k,l ` q` ,p p T h q , S D h p φ q À h min p k,l ` q` } φ } W min p k,l `q` ,p p T h q . (30) Moreover, @ ψ P W p p div; Ω q X W k ` ,p p T h q d , W D h p ψ q À h k ` } ψ } W k ` ,p p T h q d . (31) Here, for an integer s ě and a real number q P r , `8q , W s,q p T h q – t v P L q p Ω q : v | T P W s,q p T q @ T P T h u is the broken space on T h constructed on W s,q and endowed with the norm } v } W s,q p T h q – ˜ ÿ T P T h } v | T } qW s,q p T q ¸ q . Proof.
See Appendix A.3.2Appendix A.3.2.
Remark
11 (Order of S D h ) . In the case l ě k , (3030) yields the optimal order O p h k ` q for interpolations ofsmooth enough functions. If l “ k ´
1, one order is lost and (3030) gives an O p h k q estimate (but, as shownby (2929), this loss is only perceptible on the approximations of the functions, not of their gradients).10 .6. Local stabilising contribution based on a Raviart–Thomas–N´ed´elec subspace We construct in this section a stabilising contribution that fulfils the requirements expressed by (S1)(S1) – (S3)(S3) . Let, for the sake of brevity, δ k ∇ ,T – ∇ r k ` T ´ G kT . We start by observing that it holds, for all v T P U k,lT and all φ P P k p T q d , p δ k ∇ ,T v T , φ q T “ p ∇ r k ` T v T , φ q T ´ p G kT v T , φ q T “ p v T ´ r k ` T v T , ∇ ¨ φ q T ` ÿ F P F T p r k ` T v T ´ v F , φ ¨ n T F q F “ p π ,lT p v T ´ r k ` T v T q , ∇ ¨ φ q T ` ÿ F P F T p π ,kF p r k ` T v T ´ v F q , φ ¨ n T F q F “ ´p δ lT v T , ∇ ¨ φ q T ` ÿ F P F T p δ kT F v T , φ ¨ n T F q F “ p ∇ δ lT v T , φ q T ` ÿ F P F T pp δ kT F ´ δ lT q v T , φ ¨ n T F q F , where we have used the definition of δ k ∇ ,T in the first line, an integration by parts together with thedefinition (1717) of G kT in the second line, (1010) together with the fact that ∇ ¨ φ P P k ´ p T q Ă P l p T q since l ě k ´ φ | F ¨ n T F P P k p F q for all F P F T to introduce the projectors in the third line, thedefinition (2222) of δ lT and δ kT F in the fourth line, and an integration by parts to conclude. Rearrangingthe terms, we arrive at pp δ k ∇ ,T ´ ∇ δ lT q v T , φ q T “ ÿ F P F T pp δ kT F ´ δ lT q v T , φ ¨ n T F q F . (32)A few remarks are of order to illustrate the consequences of the above relation. Remark
12 (Control of the element-based difference through face-based differences) . A first notableconsequence is that the element-based difference p δ k ∇ ,T ´ ∇ δ lT q v T can be controlled in terms of the face-based differences tp δ kT F ´ δ lT q v T : F P F T u : For all T P T h and all v T P U k,lT , it holds }p δ k ∇ ,T ´ ∇ δ lT q v T } L p T q d À | v T | , B T , (33)where a À b means a ď Cb with real number C ą h and of T , but possibly dependingon d , (cid:37) , k , and l . To prove (3333), it suffices to observe that }p δ k ∇ ,T ´ ∇ δ lT q v T } L p T q d “ sup φ P P k p T q d , } φ } L p T q d “ pp δ k ∇ ,T ´ ∇ δ lT q v T , φ q T “ sup φ P P k p T q d , } φ } L p T q d “ ÿ F P F T pp δ kT F ´ δ lT q v T , φ ¨ n T F q F ď sup φ P P k p T q d , } φ } L p T q d “ | v T | , B T h T } φ ¨ n T } L pB T q À | v T | , B T , where we have used the fact that p δ k ∇ ,T ´ ∇ δ lT q v T P P k p T q d in the first line, (3232) in the second line, theCauchy–Schwarz inequality in the third line, and the discrete trace inequality (A.6A.6) below with p “ h T } φ ¨ n T } L pB T q À } φ } L p T q d and conclude. 11 x T P T F T Figure 1: Illustration of P TF . Remark
13 (Stabilisation based on a lifting of face-based differences) . Let now v T P U k,lT be fixed.Relation (3232) no longer holds true in general if we replace φ by a function η belonging to a space S T larger than P k p T q d . It then makes sense to define the nontrivial residual linear form R T p v T ; ¨q : S T Ñ R such that R T p v T ; η q – ´pp δ k ∇ ,T ´ ∇ δ lT q v T , η q T ` ÿ F P F T pp δ kT F ´ δ lT q v T , η ¨ n T F q F . Assume now S T large enough for the L p T q d -norm of the Riesz representation L T v T P S T of R T p v T ; ¨q to control | v T | , B T (hence also }p δ k ∇ ,T ´ ∇ δ lT q v T } L p T q d by (3333)). Property (S1)(S1) is then fulfilled lettingthe stabilising contribution in (2525) be such that S T v T “ L T v T for all v T . This choice also satisfies (S2)(S2) by construction since P k p T q d Ă S T and R T p v T ; ¨q vanishes on P k p T q d for any v T P U k,lT owingto (3232). Finally, property (S3)(S3) is satisfied provided that S T is a piecewise polynomial space on a regularpolytopal partition of T .The above procedure can be interpreted as a lifting on S T of the face-based differences tp δ kT F ´ δ lT q v T : F P F T u realised by means of the operator L T . This interpretation justifies the terminology employed inSection 3.6.33.6.3 below. In this section we define a good candidate to play the role of the space S T in Remark 1313. From thispoint on, we work on a fixed mesh element T P T h and assume, for the sake of simplicity, that (i) thefaces of T are p d ´ q -simplices and that (ii) T is star-shaped with respect to a point x T whose ortogonaldistance d T F from each face F P F T satisfies d T F ě (cid:37)h T , (34)where, as in Section 2.22.2, (cid:37) denotes the mesh regularity parameter. These assumptions can be relaxedusing a simplicial submesh of T and at the price of a heavier notation. For all F P F T , we denote by P T F the d -simplex of base F and apex x T , and by F T F the set of p d ´ q -simplicial faces of P T F , see Figure11. In what follows, we work on the face-based simplicial partition P T – t P T F : F P F T u .For an integer m ě F P F T , we let RT m p P T F q – P m p P T F q d ` x P m p P T F q denote theRaviart–Thomas–N´ed´elec space [4444, 4343] of degree m on the simplex P T F . Each function η P RT m p P T F q is uniquely identified by the following degrees of freedom (see, e.g., [55, Proposition 2.3.4]): p η ¨ n σ , q q σ @ σ P F T F , @ q P P m p σ q , p η , χ q P TF @ χ P P m ´ p P T F q d , σ P F T F , the normal n σ points out of P T F . Additionally, we note the following relation,valid for all η P RT m p P T F q : } η } L p P TF q d » } π ,m ´ P TF η } L p P TF q d ` ÿ σ P F TF h F } η ¨ n σ } L p σ q , (35)where » means Ca ď b ď C ´ a with real number C ą h and of P T F , but possiblydepending on (cid:37) and m .The candidate to play the role of the space S T in Remark 1313 is RT k ` p P T q , the broken Raviart–Thomas–N´ed´elec space of degree p k ` q on the submesh P T . We are now ready to construct the lifting of face-based differences. Owing to the specific choice of S T , we can proceed face by face. Specifically, for all F P F T , we define the lifting operator L k ` T F : U k,lT Ñ RT k ` p P T F q such that, for all v T P U k,lT , L k ` T F v T satisfies for all η P RT k ` p P T F qp L k ` T F v T , η q P TF “ ´pp δ k ∇ ,T ´ ∇ δ lT q v T , η q P TF ` pp δ kT F ´ δ lT q v T , η ¨ n T F q F . (36)In what follows, we extend L k ` T F v T by zero outside P T F . Proposition 14 (Stabilisation based on a Raviart–Thomas–N´ed´elec subspace) . The following stabilisingcontribution satisfies properties (S1)(S1) – (S3)(S3) : S T – ÿ F P F T L k ` T F . (37) Proof. (i)
Proof of (S1)(S1) . We abridge by a À b the inequality a ď Cb with real number C independent ofboth h and T , but possibly depending on d , (cid:37) , k , and l . We start by proving that | v T | , B T À } S T v T } L p T q d . (38)Let η P RT k ` p P T F q be such that p η ¨ n T F , q q F “ h ´ F pp δ kT F ´ δ lT q v T , q q F @ q P P k ` p F q , p η ¨ n σ , q q σ “ @ σ P F T F zt F u , @ q P P k ` p σ q , p η , χ q P TF “ @ χ P P k p P T F q d . (39)Plugging this definition into (3636) and using the Cauchy–Schwarz inequality, we infer that h ´ F }p δ kT F ´ δ lT q v T } L p F q ď } L k ` T F v T } L p P TF q d } η } L p P TF q d À } L k ` T F v T } L p P TF q d h ´ F }p δ kT F ´ δ lT q v T } L p F q , where we have used (3535) and (3939) with q “ η ¨ n T F P P k ` p F q or q “ η ¨ n σ P P k ` p σ q to estimate the L -norm of η . Dividing by h ´ F }p δ kT F ´ δ lT q v T } L p F q , squaring, summing over F P F T , and taking thesquare root of the resulting inequality proves (3838).Let us now prove that } S T v T } L p T q d À | v T | , B T . (40)13etting in (3636) η “ L k ` T F v T , summing over F P F T , and using multiple times the Cauchy–Schwarzinequality, we obtain } S T v T } L p T q d “ ÿ F P F T } L k ` T F v T } L p P TF q d ď ˆ }p δ k ∇ ,T ´ ∇ δ lT q v T } L p T q d ` ÿ F P F T h ´ F }p δ kT F ´ δ lT q v T } L p F q ˙ ˆ ˆ } S T v T } L p T q d ` h T } S T v T ¨ n T } L pB T q ˙ À | v T | , B T } S T v T } L p T q d , where we have used (3333) and the discrete trace inequality (A.6A.6) below with p “ (S1)(S1) .(ii) Proof of (S2)(S2) . Let φ P P k p T q d , set η “ φ | P TF P P k p P T F q d Ă RT k ` p P T F q in (3636), and sum over F P F T to obtain p S T v T , φ q T “ ÿ F P F T p L k ` T F v T , φ q P TF “ ´pp δ k ∇ ,T ´ ∇ δ lT q v T , φ q T ` ÿ F P F T pp δ kT F ´ δ lT q v T , φ ¨ n T F q F “ , where we have used (3232) to conclude.(iii) Proof of (S3)(S3) . The regularity of the face-based partition P T follows from (3434) in view of [2424,Lemma 3]. Property (S3)(S3) is then satisfied with k S – k ` Remark
15 (The lowest-order case) . The lifting operators L T F and stabilisation S T described above aresome of the possible choices that satisfy (S1)(S1) – (S3)(S3) . In certain cases, simpler liftings can be designed.Consider for example k “ l P t´ , u . In this case, (i) δ k ∇ ,T v T “ since P p T q d “ ∇ P p T q , so that G T “ ∇ r T ; (ii) δ ´ T “ δ T “ p L T F v T , η q P TF “ p δ kT F v T , η ¨ n T F q F . (41)An appropriate lifting can then be constructed in P p P T F q d instead of RT p P T F q as described hereafter,and the assumption that the faces of T are simplices can be removed (for any F P F T , P T F is then thepyramid with base F and apex x T ). Define L HMM
T F : U ,lT Ñ P p P T F q d such that, for all v T P U ,lT , L HMM
T F v T “ | F || P T F | δ T F v T n T F “ dd T F δ T F v T n T F . Notice that this lifting is designed to satisfy (4141) for all η P P p P T F q d . Defining, for all T P T h , S T by(3737), the discrete elements of the corresponding DSGD D are then X D , “ t v h “ pp v T q T P T h , p v F q F P F h q : v T P R , v F P R , v F “ @ F P F b h u , and, for all v h P X D , and all T P T h , expressing (1717), (2121) and (2222) for k “ p Π D v h q | T “ v T , p ∇ D v h q | P TF “ G T v h ` dd T F ` v T ` G T v T ¨ p x F ´ x T q ´ v F ˘ n T F @ F P F T where x F is the centre of mass of F and G T v T – | T | ř F P F T | F | v F n T F .14 igure 2: Triangular, Cartesian, hexagonal and locally refined meshes for the numerical examples of Section 3.73.7
This gradient discretisation corresponds, for k “ l “
0, to the Hybrid-Mixed-Mimetic methods of [2929](with the isomorphism A T “ ´ ? d Id in [2929, Section 5.3.1]). The coercivity, GD-consistency, limit–conformity and compactness of families of such GDs are established in [2828, Chapter 13], and can also beproved by checking that S T satisfies (S1)(S1) – (S3)(S3) (use η “ h ´ F δ T F v T n T F in (4141) to establish Á in (2626),and [2828, Eq. (13.10)] to see that ş T S T v T “ and thus that (S2)(S2) holds). In this section we provide numerical evidence to support our theoretical results. p ě “ p , q the p -Laplace problem (66) correspondingto the exact solution u p x q “ sin p πx q sin p πx q , (42)with p P t , , u and source term inferred from u . We consider the matching triangular, Cartesian,(predominantly) hexagonal, and locally refined mesh families depicted in Figure 22 and polynomial degreesranging from 0 to 4. The first, second, and fourth mesh families originate from the FVCA5 benchmark [3939],whereas the third from [2424]. The local refinement in the third mesh family has no specific meaning here:its purpose is to demonstrate the seamless treatment of nonconforming interfaces.We report in Figure 33 the error } G kh p I k,lh u ´ u h q} L p p Ω q d versus the meshsize h , where G kh is theconsistent (but in general not stable) global gradient reconstruction defined by (2020). The referenceslopes correspond to the convergence rates resulting from Theorem 33 together with Proposition 1010 (moreprecisely, the order corresponds to the dominating term). The theoretical orders of convergence areperfectly matched for p “
2. Similar considerations hold for p “ k ă k ě ψ ÞÑ | ψ | p ´ ψ , which impacts onthe regularity of | ∇ u | p ´ ∇ u . Finally, for p “ k P t , , u , while faster convergence than predicted by the error estimates is observed for k P t , u .This phenomenon will be further investigated in the future. For a comparison with the HHO method of[1616], see Figure 55 below. p ă ă p ă
2. For this reason, we consider instead theexponential solution u p x q “ exp p x ` πx q , and solve the p -Laplace problem with p “ , Dirichlet boundary conditions on B Ω, and right-hand side f inferred from the expression of u . The error } G kh p I k,lh u ´ u h q} L p p Ω q d versus the meshsize h is plotted inFigure 44 for the mesh families illustrated in Figure 22 and polynomial degrees ranging from 0 to 4. Forthe sake of completeness, a comparison with the HHO method (4949) is also included. We observe thatthe DSGD and HHO methods give similar results in terms of the selected error measure (which accounts15 “ k “ k “ k “ k “ k “ k “ k “ k “ k “ ´ . ´ ´ . ´ ´ ´ ´ ´ (a) p “
2, triangular ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (b) p “
3, triangular ´ . ´ ´ . ´ ´ ´ ´ ´ ´ ´ (c) p “
4, triangular ´ . ´ ´ . ´ ´ ´ ´ ´ (d) p “
2, Cartesian ´ . ´ ´ . ´ ´ ´ ´ ´ (e) p “
3, Cartesian ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (f) p “
4, Cartesian ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ ´ (g) p “
2, hexagonal ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ ´ (h) p “
3, hexagonal ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ ´ ´ (i) p “
4, hexagonal ´ . ´ ´ . ´ ´ ´ ´ ´ (j) p “
2, locally refined ´ . ´ ´ . ´ ´ ´ ´ ´ (k) p “
3, locally refined ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (l) p “
4, locally refined
Figure 3: } G kh p I k,lh u ´ u h q} L p p Ω q d v. h . Trigonometric test case, p P t , , u , DSGD. “ k “ k “ k “ k “ k “ k “ k “ k “ k “ ´ . ´ ´ . ´ ´ ´ ´ ´ (a) Triangular, DSDG ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (b) Triangular, HHO ´ . ´ ´ . ´ ´ ´ ´ ´ (c) Cartesian, DSGD ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (d) Cartesian, HHO ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ ´ (e) Hexagonal, DSGD ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ ´ (f) Hexagonal, HHO ´ . ´ ´ . ´ ´ ´ ´ ´ (g) Locally refined, DSGD ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (h) Locally refined, HHO Figure 4: Error } G hh p I k,lh u ´ u h q} L p p Ω q d v. h , exponential test case, p “ .
4. Alternative gradient and links with other methods
In this section we discuss an alternative to the gradient reconstruction G T defined by (2525), and linkswith other methods. r k ` T An alternative to using G kT in (2525) would be to use ∇ r k ` T . This would lead to a local gradient definedby p r ∇ D h v h q | T “ r G T v T : “ ∇ r k ` T v T ` r S T v T @ v h P U k,lh , @ T P T h . (43)The stabilisation term r S T is still required to satisfy the L -stability and boundedness property (S1)(S1) and the image property (S3)(S3) . As shown in (A.13A.13) below, the property (S2)(S2) is required on S T v T toensure that it is orthogonal to the consistent part G kT v T of G T v T . For r G T v T , the consistent part is ∇ r k ` T v T P ∇ P k ` p T q , and the orthogonality property on r S T v T can therefore be relaxed into ( Ă S2)
For all v T P U k,lT and all φ P ∇ P k ` p T q , p r S T v T , φ q T “ (S2)(S2) , the space for φ is ∇ P k ` p T q instead of P k p T q d .Performing integration-by-parts and introducing projection operators in (21a21a) leads to ´ p ∇ δ lT v T , ∇ w q T “ ÿ F P F T pp δ kT F ´ δ lT q v T , ∇ w ¨ n T F q T , @ w P P k ` p T q . (44)Comparing with (3232), we see that the volumetric term δ k ∇ ,T v T has disappeared, and the only remainingvolumetric term ∇ δ lT v T belongs to ∇ P l p T q . This is a gain with respect to the situation in Section 3.63.6,in which the degree of the volumic term was constrained by δ k ∇ ,T v T to be k , whatever the choice of l .As a consequence, the construction of r S T can be done in a (possibly) smaller space than RT k ` p P T F q ,as detailed in what follows.First, considering w “ ´ δ lT v T in (4444) and using the trace inequality (A.6A.6) with p “ ∇ δ lT v T yields (compare with Remark 1212) } ∇ δ lT v T } L p T q d À | v T | , B T . (45)Then, following the ideas of Section 3.6.33.6.3, we construct the lifting r L kT F : U k,lT Ñ RT max p l,k q p P T F q suchthat, for all η P RT max p l,k q p P T F q , p r L kT F v T , η q P TF “ p ∇ δ lT v T , η q P TF ` pp δ kT F ´ δ lT q v T , η ¨ n T F q F . (46)Since ∇ δ lT v T P ∇ P l p T q Ă P l ´ p T q d and p δ kT F ´ δ lT q v T P P max p l,k q p F q , following the proof of (S1)(S1) inProposition 1414 shows that, with this choice of r L kT F , we only need (3939) to hold for q P P max p l,k q p F q and χ P P l ´ p T q d . The space RT max p l,k q p P T F q enables these choices of q and χ , and the L -stability of r S T : “ ř F P F T r L kT F therefore follows. When l P t k ´ , k u , this stabilisation term is constructed on apiecewise RT k space, instead of a piecewise RT k ` space for S T in Section 3.6.33.6.3.Although the choice (4343) leads to coercive, consistent, limit-conforming, and compact families ofgradient discretisations, it can turn out to be far from optimal for general problems. More precisely, itslimit-conformity properties are much worse than those of (2525). Indeed, the full orthogonality property (S2)(S2) is essential to establish the O p h k ` q estimate (3131) on W D h (see Remark 2424). In general, we cannot18stablish more than an O p h q estimate on W D h , irrespective of k , if the gradient is reconstructed via(4343). For anisotropic linear problems, a modification of r k ` T embedding a dependence on the diffusioncoefficient can be constructed to recover optimal rates of convergence [2121]; for fully non-linear models,though, the only option to recover a truly high-order method seems to be using the gradient G kT v T inthe full polynomial space P k p T q d .To illustrate numerically the loss of convergence experienced when using the gradient reconstruc-tion (4343) in the context of fully non-linear problems, we solve the problem described in Section 3.7.13.7.1using two numerical methods: the HHO method of [1616] (see (4949) below), and the method obtained fromthe latter replacing G kT by ∇ r k ` T . We report in Figure 55 the error } I k,lh u ´ u h } ,p,h versus the meshsize h , with reference slopes corresponding to the estimates of convergence rates derived in [1717, Theorem 3.2].The leftmost column, corresponding to the Poisson problem with p “
2, shows that both G kT and ∇ r k ` T can be used when W D h is applied to ψ “ ∇ u . For p P t , u , on the other hand, a significant loss in theconvergence rate is observed for k ą ∇ r k ` T departs from the solidline corresponding to G kT ). Notice that, for k ě p “
3, the order of convergence is again limited bythe regularity of the function ψ ÞÑ | ψ | p ´ ψ . The results presented here replace the ones of [1717, Figure3], which were affected by a bug in one of the libraries used in our code. From these new tests, the errorestimates of [1717, Theorem 3.2] appear to be sharp also for p ą The HHO method proposed in [2222] for problem (22) with Λ “ I d readsFind u h P U k,lh, such that, for all v h P U k,lh, , a hho h p u h , v h q “ p f, v h q , (47)where the broken polynomial function v h is defined by (1515), and the bilinear form a hho h : U k,lh ˆ U k,lh Ñ R is assembled from the elementary contributions a hho T p u T , v T q – p ∇ r k ` T u T , ∇ r k ` T v T q T ` ÿ F P F T h ´ F pp δ kT F ´ δ lT q u T , p δ kT F ´ δ lT q v T q F . (48)Here, the consistent gradient is ∇ r k ` T , as in (4343), and the stabilisation is not incorporated in the gradientreconstruction, but rather added as a separate term in the bilinear form.For the non-linear problem (66), on the other hand, the HHO method considered in [1616, 1717] readsFind u h P U k,lh, such that, for all v h P U k,lh, , A hho h p u h , v h q “ p f, v h q ,with function A hho h : U k,lh ˆ U k,lh Ñ R assembled from the elementary contributions A hho T p u T , v T q – ż T σ p G kT u T q¨ G kT v T ` ÿ F P F T h p ´ F ż F |p δ kT F ´ δ lT q u T | p ´ p δ kT F ´ δ lT q u T p δ kT F ´ δ lT q v T . (49)Also in this case, stability is achieved by a separate term, and only the consistent (but not stable)gradient reconstruction appears in the consistency term. Notice that, unlike the linear case, the gradientreconstructed in the full polynomial space P k p T q d is present here; see comments at the end of Section4.14.1. The ncMFD method of [4141] hinges on degrees of freedom (DOFs) that are the polynomial momentsof degree up to p k ´ q inside the mesh elements and the polynomial moments of degree up to k on themesh faces. X h is the space of vectors gathering such DOFs. Given a function v P H p Ω q , we denote by19 “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) k “ G kT ) k “ ∇ r k ` T ) ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (a) p “
2, triangular ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (b) p “
3, triangular ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (c) p “
4, triangular ´ . ´ ´ . ´ ´ ´ ´ ´ (d) p “
2, Cartesian ´ . ´ ´ . ´ ´ ´ ´ (e) p “
3, Cartesian ´ . ´ ´ . ´ ´ ´ ´ ´ ´ (f) p “
4, Cartesian ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ ´ (g) p “
2, hexagonal ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ (h) p “
3, hexagonal ´ . ´ . ´ ´ . ´ . ´ . ´ . ´ ´ ´ ´ (i) p “
4, hexagonal ´ . ´ ´ . ´ ´ ´ ´ ´ (j) p “
2, locally refined ´ . ´ ´ . ´ ´ ´ ´ (k) p “
3, locally refined ´ . ´ ´ . ´ ´ ´ ´ ´ (l) p “
4, locally refined
Figure 5: } G kh p I k,lh u ´ u h q} L p p Ω q d v. h . Trigonometric test case, p P t , , u , HHO. I its interpolation in X h , that is, the vector collecting its moments in the elements and on the faces. Ifv P X h , we denote by v T P P k ´ p T q , T P T h , and v F P P k p F q , F P F h , the polynomials reconstructedfrom the moments represented by v. This defines an isomorphismv P X h ÞÑ v h “ pp v T q T P T h , p v F q F P F h q P U k,k ´ h . (50)Denoting by v | T the sub-vector made of the DOFs of v P X h in the mesh element T P T h and on themesh faces in F T , the ncMFD method for problem (22) with Λ “ I d readsFind u P X h such that, for all v P X h , ÿ T P T h u t | T M T v | T “ ÿ T P T h L f p v | T q , (51)where L f p v T q is a discretisation of ş T f v , and the matrix M T is positive semi-definite with suitableconsistency and stability properties. Setting N d,k “ dim p P k ´ p T qq and selecting a basis p q i q i “ ,...,N d,k ´ of P k ´ p T q with q “
1, the required consistency and stability properties on M T lead to the followingdecomposition (see [4141, Eq. (35)]): M T “ M T ` M T “ p R T p p R tT p N q ´ p R tT ` M T where A t is the transpose of A , ker p M T q “ tp v I q | T : v P P k ´ p T qu , p N has columns p N i “ p q i q I for i “ , . . . , N d,k ´
1, and p R T is the matrix with columns p p R T,i q i “ ,...,N d,k defined by @ v P X h , p p R T,i q t v | T “ ´p v T , ∆ q i q ` ÿ F P F T p v F , ∇ q i ¨ n T F q F . (52)The stabilising matrix M T does not have any impact on the consistency of the method; its sole role isto stabilise the matrix M so that its kernel is R I (which is expected: a matrix M T representing thebilinear form ş T ∇ u ¨ ∇ v should vanish on interpolants of constant functions). The matrix M T , on theother hand, contains all the consistency properties of the method. The analysis of the stabilisation partis made in Section 4.54.5, alongside the analysis of the stabilisation in the HHO and ncVEM methods.Let us analyse here the consistent part M T of M T . Take v P X h , and let v h P U k,k ´ h correspondingto v through the isomorphism (5050). Comparing (21a21a) and (5252) shows that p ∇ r k ` T v T , ∇ q i q T “ p p R T,i q t v | T “ p p R tT v | T q i , @ i “ , . . . , N d,k ´ . Hence, with obvious notations,v t | T M T w | T “ p p R tT v | T q t p p R tT p N q ´ p p R tT w | T q“ rp ∇ r k ` T v T , ∇ q i q T s ti “ ,...,N d,k ´ p p R tT p N q ´ rp ∇ r k ` T w T , ∇ q i q T s i “ ,...,N d,k ´ . (53)Applying (5252) to v | T “ p q j q I and integrating by parts shows that p p R tT p N q ij “ p ∇ q i , ∇ q j q T , that is, p R tT p N is the Gram matrix of p ∇ q i q i “ ,...,N d,k ´ in L p T q d . Equation (5353) can therefore be re-written asv t | T M T w | T “ p ∇ r k ` T v T , ∇ r k ` T w T q T . Thus, recalling the isomorphism (5050) between X h and U k,k ´ h , the ncMFD method (5151) is equivalent toFind u h P U k,k ´ h such that, for all v h P U k,k ´ h , ÿ T P T h p ∇ r k ` T u T , ∇ r k ` T v T q T ` s mfd T p u , v q “ ÿ T P T h L f p v | T q , where s mfd T is the stabilising bilinear form assembled from the local matrices p M T q T P T h . Here, the consis-tent part is constructed from ∇ r k ` T , as in (4343), but the stabilisation is external to the gradient. In thencMFD method, the loading term is discretised (for k ě
2) by L f p v | T q “ ż T π ,k ´ T p f q v T “ ż T f v T . (54)A modification of (5454) is necessary for the low orders k “ , .4. Non-conforming Virtual Element Method For a given mesh element T P T h , we define the local non-conforming virtual element space asfollows [22]: V h,nck p T q – (cid:32) v h P H p T q : n T F ¨ ∇ v h | F P P k p F q @ F P F T , ∆ v h P P k ´ p T q ( . (55)This space is finite dimensional and its functions are described by their face moments against the poly-nomials of degree up to k on each face F P F T and the cell moments against polynomials of degree up to k ´ T . The unisolvency of these degrees of freedom has been proved in [22].The global non-conforming virtual element space of degree k is given by: V h,nck, p Ω q – " v h P L p Ω q : ż F rr v h ss q “ @ F P F h , @ q P P k p F q and v h | T P V h,nck p T q @ T P T h * , (56)where rr v h ss denotes the jump operator with the usual definition at interfaces (the sign is not relevant),and extended to boundary faces setting rr v h ss – v h .The polynomials of degree up to p k ` q are a subspace of V h,nck p T q . The jump conditions on theelements of V h,nck, p Ω q ensure that I k,k ´ h : V h,nck, p Ω q Ñ U k,k ´ h, is well defined (there is only one interpolanton each face, and the polynomial moments up to degree k of functions in V h,nck, p Ω q vanish on each F P F b h ). Moreover, the unisolvent property of the set of DOFs shows that I k,k ´ h is an isomorphismbetween V h,nck, p Ω q and U k,k ´ h, .The virtual element discretisation of problem (22) with Λ “ Id reads as:Find u h P V h,nck, p Ω q such that, for all v h P V h,nck, p Ω q , a h p u h , v h q “ p f h , v h q (57)where the bilinear form a h p u h , v h q approximates the left-hand side of (22) and f h | T “ π ,k ´ T f .Mimicking the additivity of integrals, we assume that the virtual element bilinear form is the sum-mation of local elemental terms a h p u h , v h q “ ÿ T P T h a h,T p u h , v h q . (58)Two different formulations of the local bilinear form a h,T can be found in the literature, both includinga consistency and a stability term: – first formulation [22]: for every u h , v h P V h,nck p T q : a h,T p u h , v h q – p ∇ π ,k ` T u h , ∇ π ,k ` T v h q T ` S ` p Id ´ π ,k ` T q u h , p Id ´ π ,k ` T q v h ˘ ; (59) – second formulation [1111, 1010]: for every u h , v h P V h,nck p T q : a h,T p u h , v h q – p π ,kT ∇ u h , π ,kT ∇ v h q T ` S ` p Id ´ π ,k ` T q u h , p Id ´ π ,k ` T q v h ˘ . (60)In these definitions, the first term on the right is the consistency term designed to provide the exactnessof the integration whenever at the least one of the entries u h or v h is a polynomial of degree up to p k ` q .The second term is a stabilisation, and S can be any symmetric and positive definite bilinear form forwhich there exist two positive constants c ˚ and c ˚ such that c ˚ } ∇ v h } L p T q d ď S p v h , v h q ď c ˚ } ∇ v h } L p T q d @ v h P V h,nck p T q such that π ,k ` T v h “ . (61) Remark . The connection between formulation (5959) and the ncMFD method of [4141] reviewed in Sec-tion 4.34.3 has been established in [22]. 22 .5. HHO, ncMFD and ncVEM are gradient discretisation methods
In the previous sections, we showed that the consistent gradient in the HHO, ncMFD and ncVEMmethods (first formulation) is, for the Poisson problem, identical to the consistent gradient in (4343).We show here that the stabilisations used in these methods can actually be represented by well-chosenstabilisation terms S T satisfying (S1)(S1) – (S3)(S3) , and thus that these methods are gradient discretisationmethods.Following the discussion in the previous sections, for the Poisson problem the HHO, ncMFD andncVEM (first formulation) can be written, upon an isomorphism of the space of discrete unknowns andwith the proper choice of l P t k ´ , k, k ` u : Find u h P U k,lh, such that, for all v h P U k,lh, , ÿ T P T h p ∇ r k ` T u T , ∇ r k ` T v T q T ` ÿ T P T h Stab T p u T ´ I k,lT r k ` T u T , v T ´ I k,lT r k ` T v T q “ ÿ T P T h L T,f p v T q , (62)where L T,f : U k,lT Ñ R is a specific linear form and Stab T is a symmetric bilinear form on U k,lT that iscoercive and stable on Im p Id ´ I k,lT r k ` T q , that is, @ v T P U k,lT , Stab T p v T ´ I k,lT r k ` T v T , v T ´ I k,lT r k ` T v T q » | v T | , B T , (63)where a » b means Ca ď b ď C ´ a with real number C ą h and of T . Remark
17 (Stabilisation term) . For the HHO scheme (4848), by (2323) the vector pp δ kT F ´ δ lT q v T q F P F T isthe difference of the face- and cell-unknowns of v T ´ I k,lT r k ` T v T . Hence, the definitionStab T p u T ´ I k,lT r k ` T u T , v T ´ I k,lT r k ` T v T q – ÿ F P F T h ´ F pp δ kT F ´ δ lT q u T , p δ kT F ´ δ lT q v T q F is valid and ensures that the stabilisation terms in (4848) and (6262) coincide. This also shows that Stab T p v T ´ I k,lT r k ` T v T , v T ´ I k,lT r k ` T v T q “ | v T | , B T , and thus that (6363) holds.For the ncVEM version (any formulation), comparing the stabilisations in (6060)–(5959) and in (6262) leadsto defining, for all u h , v h P V h,nck p T q and setting u T “ I k,lT u h | T and v T “ I k,lT v h | T ,Stab T p u T ´ I k,lT r k ` T u T , v T ´ I k,lT r k ` T v T q – S ` p Id ´ π ,k ` T q u h , p Id ´ π ,k ` T q v h ˘ . Estimate (6363) then follows from (6161) and from (A.17A.17) in Lemma 2222.Due to the orthogonality condition (S2)(S2) , the GS (33) for Problem (11) with Λ “ I d , based on thegradient reconstructions (4343) (with r S T “ S T ) and some function reconstruction r Π D h , is given by ÿ T P T h p ∇ r k ` T u T , ∇ r k ` T v T q T ` ÿ T P T h p S T u T , S T v T q T “ ż Ω f r Π D h v h . (64)For all methods except the ncMFD with k “
1, accounting for (1414) when k “ l “ ´
1, we have L f,T p v T q “ p f, v T q T and thus the choice r Π D h “ Π D h defined in (1616) ensures that the right-hand sides of(6262) and (6464) coincide. For the ncMFD with k “
1, a slightly different discretisation of the right-handside has to be considered in order to ensure optimal L -error estimates under elliptic regularity; see [4141,Section 2.7] for further details.To prove that (6262) can be written as (6464), it remains to show that for any stabilisation Stab T asabove, there exists S T satisfying (S1)(S1) – (S3)(S3) and such that @ u T , v T P U k,lT , p S T u T , S T v T q T “ Stab T p u T ´ I k,lT r k ` T u T , v T ´ I k,lT r k ` T v T q . (65)Let us fix an initial S T satisfying the design properties (for example, the stabilisation defined by (3737)). Theproperty (S1)(S1) on S T and Lemma 2121 in Appendix Appendix A.2Appendix A.2 below show that S T : U k,lT Ñ Im p S T q p Id ´ I k,lT r k ` T q : U k,lT Ñ Im p Id ´ I k,lT r k ` T q have the same kernel. As shown by (6363) and Lemma2121, Stab T is an inner product on Im p Id ´ I k,lT r k ` T q . Applying thus [2929, Lemma A.3] produces an innerproduct x¨ , ¨y T on Im p S T q such that @ u T , v T P U k,lT , x S T u T , S T v T y T “ Stab T p u T ´ I k,lT r k ` T u T , v T ´ I k,lT r k ` T v T q . (66)Using then [3030, Lemma 5.2] with the inner products x¨ , ¨y T and p¨ , ¨q T on Im p S T q , we find an isomorphism L T of Im p S T q such that @ u T , v T P U k,lT , x S T u T , S T v T y T “ p L T S T u T , L T S T v T q T . (67)Combining (6666) and (6767) shows that (6565) holds with S T “ L T S T . The proof that S T defined above satisfies the design properties is easy. The L -stability and boundedness (S1)(S1) is a direct consequence of (6565) with v T “ u T and of (6363). The image of S T is, by assumption, L p T q d -orthogonal to P k p T q d and contained in some P k S p P T q . Since L T is an isomorphism of Im p S T q ,we have Im p S T q “ Im p S T q and Properties (S2)(S2) and (S3)(S3) on S T therefore follow. Remark
18 (Formulations based on G kT ) . With minor modifications, it is easy to construct a stabilisation S T that satisfies ( Ă S2) instead of (S2)(S2) . This can be done at a reduced cost, as discussed in Section 4.14.1.The reasoning above can also be easily adapted to the first formulation of ncVEM, provided that theterms ∇ r k ` T in (6262) are replaced with G kT . Appendix A. Proofs of the results on DSGDs
This section contains the proofs of Theorem 99 and Proposition 1010 preceded by some preliminaryresults: the study of the properties of stabilising contributions satisfying (S1)(S1) – (S3)(S3) and uniform equiv-alences of discrete W ,p -seminorms. We also include lemmas used in Section 44 to show that the HHOmethod, ncMFD method and ncVEM are GDMs. Appendix A.1. Properties of the stabilising contribution
Proposition 19 (Properties of S T ) . Let t S T : T P T h u be a family of stabilising contributions satisfyingassumptions (S1)(S1) – (S3)(S3) . Then, the following properties hold:(i) L p -stability and boundedness. For all T P T h and all v T P U k,lT , } S T v T } L p p T q d » | v T | p, B T , (A.1) with hidden constant as in (2626) and additionally depending on p and k S .(ii) Consistency.
For all T P T h and all v P W k ` ,p p T q , } S T I k,lT v } L p p T q d À h k ` T | v | W k ` ,p p T q , (A.2) where a À b means a ď Cb with real number C ą independent of both h and T , but possiblydepending on d , (cid:37) , k , l , p and k S .As a consequence of (A.2A.2) , if v P P k ` p T q , then S T I k,lT v “ .
24n the proof, we will need the following direct and reverse Lebesgue embeddings, proved in [1616,Lemma 5.1]: Let X denote a measurable subset of R d with inradius r X and diameter h X , and let tworeals r, s P r , `8s and an integer (cid:96) P N be fixed. Then, for all q P P (cid:96) p X q , it holds that } q } L r p X q » | X | r ´ s } q } L s p X q , (A.3)where a » b means Ca ď b ď C ´ a with real number C ą d , a lower bound of theratio r X h X , r , s , and (cid:96) .We will also need the following L p -trace inequality (see, e.g., [1616, Eq. (A.10)]): For all T P T h and all v P W ,p p T q , h p T } v } L p pB T q À } v } L p p T q ` h T } ∇ v } L p p T q d , (A.4)where a À b means a ď Cb with real number C ą h and of T , but possibly dependingon d , p , and (cid:37) . When v P P (cid:96) p T q for some integer (cid:96) ě
0, combining (A.4A.4) with the following inverseinequality (see, e.g., [1616, Remark A.2]): } ∇ v } L p p T q d À h ´ T } v } L p p T q , (A.5)yields h p T } v } L p pB T q À } v } L p p T q , (A.6)where the hidden multiplicative constant in (A.5A.5) and (A.6A.6) can additionally depend on (cid:96) .We are now ready to prove Proposition 1919. Proof of Proposition 1919. (i) L p -stability and boundedness. When p “
2, (A.1A.1) coincides with (2626). Letus now consider the case p ‰ n ě q P r , `8q , a i ě ď i ď n ): n ÿ i “ a qi ď ˜ n ÿ i “ a i ¸ q ď n q ´ n ÿ i “ a qi , n ´ qq n ÿ i “ a q i ď ˜ n ÿ i “ a i ¸ q ď n ÿ i “ a q i . (A.7)We also notice that, owing to (2828), it holds for all T P T h that | P | » | T | @ P P P T , card p P T q » . (A.8)To prove | P | » | T | , it suffices to observe that | P | » h dP » h dT » | T | , where we have used, respectively, thefirst and second conditions in (2828) and the mesh regularity to conclude. The bound on card p P T q followsby writing | T | “ ř P P P T | P | » ř P P P T | T | “ card p P T q| T | .Let now T P T h and v T P U k,lT be fixed. Since S T v T P P k S p P q for all P P P T , we have that } S T v T } pL p p T q d “ ÿ P P P T } S T v T } pL p p P q d » ÿ P P P T | P | p p p ´ q} S T v T } pL p P q d Eqs. (2828) and (A.3A.3) » | T | p p p ´ q ÿ P P P T } S T v T } pL p P q d Eq. (A.8A.8) . In the second line, condition (2828) is invoked to use (cid:37) as a lower bound for r P h P when applying (A.3A.3). Usingthe first pair of inequalities in (A.7A.7) with q “ p if p ě
2, the second pair of inequalities in (A.7A.7) with q “ p if p ă a i “ } S T v T } L p P q d and n “ card p P T q »
1, we infer } S T v T } pL p p T q d » | T | p p p ´ q ˜ ÿ P P P T } S T v T } L p P q d ¸ p “ | T | p p p ´ q} S T v T } pL p T q d . p th root of the above relation, we arrive at } S T v T } L p p T q d » | T | p ´ } S T v T } L p T q d . (A.9)Proceeding similarly using (A.3A.3) repeatedly on the faces of T and using | F | h F » | T | , we can provethat | v T | , B T » | T | ´ p | v T | p, B T . (A.10)Combining (A.9A.9) and (A.10A.10) with (2626), (A.1A.1) follows.(ii) Consistency.
Let a mesh element T P T h be fixed and set, for the sake of brevity, p v T – I k,lT v .Using the uniform equivalence (A.1A.1) proved in the first point, and recalling the definition (2424) of the |¨| p, B T -seminorm, it is inferred that } S T p v T } pL p p T q d À | p v T | pp, B T “ ÿ F P F T h ´ pF }p δ kT F ´ δ lT q p v T } pL p p F q . For all F P F T , using the triangle inequality followed by the discrete trace inequality (A.6A.6), we infer that }p δ kT F ´ δ lT q p v T } pL p p F q À } δ kT F p v T } pL p p F q ` } δ lT p v T } pL p p F q À } δ kT F p v T } pL p p F q ` h ´ T } δ lT p v T } pL p p T q . Using the above inequality together with the uniform bound on the number of faces of T (see [1919,Lemma 1.41]) for the second term, and expanding the difference operators according to their defini-tions (2222), we arrive at } S T p v T } pL p p T q d À h ´ pT } π ,lT p r k ` T p v T ´ v q} pL p p T q ` ÿ F P F T h ´ pF } π ,kF p r k ` T p v T ´ v q} pL p p F q . We notice that the projectors in the above bound can be removed invoking the L p p T q -boundedness of π ,lT for the first term and the L p p F q -boundedness of π ,kF for the second (see [1616, Lemma 3.2]). Combiningthis observation with the optimal approximation properties of r k ` T ˝ I k,lT discussed in Remark 55, we thenconclude that } S T p v T } pL p p T q d À h ´ pT } r k ` T p v T ´ v } pL p p T q ` ÿ F P F T h ´ pF } r k ` T p v T ´ v } pL p p F q À h p p k ` q T | v | pW k ` ,p p T q . Appendix A.2. Uniform equivalence of discrete W ,p -seminorms The second preliminary result is the uniform equivalence of various W ,p -seminorms on the globalspace of discrete unknowns U k,lh . Proposition 20 (Uniform equivalence of discrete W ,p -seminorms) . Define the discrete seminorm ~¨~ ,p,h such that, for all v h P U k,lh , ~ v h ~ p ,p,h – ÿ T P T h ~ v T ~ p ,p,T , ~ v T ~ p ,p,T – } ∇ v T } pL p p T q d ` ÿ F P F T h ´ pF } v F ´ v T } pL p p F q @ T P T h . (A.11) Then, ~¨~ ,p,h is a norm on the subspace U k,lh, and, denoting by t S T : T P T h u a family of stabilisingcontributions that satisfy (S1)(S1) – (S3)(S3) and defining ∇ D h by (2525) , it holds for all v h P U k,lh and all T P T h , ~ v T ~ ,p,T » } v T } ,p,T » } ∇ D h v h } L p p T q d , (A.12) where a » b means Ca ď b ď C ´ a with real number C ą independent of T and h , but possiblydepending on d , p , (cid:37) , k , l , and k S . As a consequence, for all v h P U k,lh , ~ v h ~ ,p,h » } v h } ,p,h » } ∇ D h v h } L p p Ω q d . (A.13)26 roof. The fact that ~¨~ ,p,h is a norm on U k,lh, can be proved in a similar manner as for the case p “ l “ k considered in [2020, Proposition 5]. The local seminorm equivalence ~ v h ~ ,p,T » } v h } ,p,T validfor all T P T h , on the other hand, is proved in [1616, Lemma 5.2] (see also references therein) for the case l “ k , and the same reasoning extends to l “ k ´ l “ k ` T P T h be fixed. We have that } ∇ D h v h } L p p T q d “ } G T v T } L p p T q d » | T | p ´ } G T v T } L p T q d “ | T | p ´ ´ } G kT v T } L p T q d ` } S T v T } L p T q d ¯ » | T | p ´ ´ } G kT v T } L p T q d ` | v T | , B T ¯ » ´ } G kT v T } L p p T q d ` | v T | p, B T ¯ » } v T } ,p,T , (A.14)where we have used a reasoning similar to the one leading to (A.9A.9) in the first line (recall that G T v T is piecewise polynomial on T owing to (S3)(S3) ), the orthogonality property (S2)(S2) in the second line, thestability and boundedness property (S1)(S1) in the third line, and (A.10A.10) together with the discrete Lebesgueembeddings (A.3A.3) in the last line. This concludes the proof of (A.12A.12). The global version (A.13A.13) followsby raising (A.12A.12) to the power p and summing over T P T h .The following lemma, which justifies the importance of the seminorm | ¨ | p, B T , was used in Section 4.54.5to prove that HHO, ncMFD and ncVEM are GDM. Lemma 21.
For any T P T h and any v T P U k,lT , it holds that ~ v T ´ I k,lT r k ` T v T ~ ,p,T » | v T | p, B T , (A.15) where a » b means C ´ a ď b ď Ca with real number C ą depending only on d , (cid:37) , p , k , and l . As aconsequence, v T ´ I k,lT r k ` T v T “ if and only if | v T | p, B T “ .Proof. Here, a À b means that a ď Cb for some real number C ą } ∇ δ lT v T } L p p T q d À | v T | p, B T . Hence, using the relation (2323) together with the definitions (A.11A.11) of ~¨~ ,p,T and (2424) of | v T | p, B T , we obtain ~ v T ´ I k,lT r k ` T v T ~ p ,p,T “ } ∇ δ lT v T } pL p p T q d ` ÿ F P F T h ´ pF }p δ kT F ´ δ lT q v T } pL p p F q » | v T | pp, B T , which is (A.15A.15). If v T ´ I k,lT r k ` T v T “
0, the relation above shows that | v T | p, B T “
0. Conversely, if | v T | p, B T “ w T – v T ´ I k,lT r k ` T v T , we have that ~ w T ~ ,p,T “
0, which implies in turn that w T is constant equal to c and that w F “ w T “ c for all F P F T . Then, r k ` T w T “ r k ` T I k,lT c “ π ,k ` T c “ c (see Remark 55 and additionally observe, for k “ l “ ´
1, that w T “ c owing to the choice of theweights ω T F ). We then have that c “ r k ` T w T “ r k ` T p v T ´ I k,lT r k ` T v T q “ r k ` T v T ´ r k ` T I k,lT r k ` T v T “ r k ` T v T ´ r k ` T v T “ , where we have used the definition of w T in the second equality, the linearity of r k ` T in the third, thefact that r k ` T I k,lT r k ` T “ r k ` T in the fourth (since r k ` T I k,lT “ π ,k ` T preserves polynomials up to degree ď k ` c is 0, which shows that v T ´ I k,lT r k ` T v T “ emma 22. For all T P T h and all v h P V h,nck p T q (with V h,nck p T q defined by (5555) ), it holds that } ∇ v h } L p T q d » ~ I k,k ´ T v h ~ , ,T , (A.16) where ~¨~ , ,T is defined in (A.11A.11) and a » b means C ´ a ď b ď Ca with real number C ą independenton h , but possibly depending on d , k , l and (cid:37) . As a consequence,for all v h P V h,nck p T q such that π ,kT v h “ , } ∇ v h } L p T q d » | I k,k ´ T v h | , ,T . (A.17) Proof.
Here, a À b means a ď Cb with C as in the statement. Since ∇ v h ¨ n T F P P k p F q for all F P F T and ∆ v h P P k ´ p T q , integrating by parts and setting v T – I k,k ´ T v h , } ∇ v h } L p T q d “ ´ p v h , ∆ v h q T ` ÿ F P F T p v h , ∇ v h ¨ n T F q F “ p´ v T , ∆ v h q T ` ÿ F P F T p v F , ∇ v h ¨ n T F q F “ p ∇ v T , ∇ v h q T ` ÿ F P F T p v F ´ v T , ∇ v h ¨ n T F q F . The Cauchy–Schwarz inequality and the trace inequality (A.6A.6) applied on ∇ v h ¨ n T F then yield } ∇ v h } L p T q d À~ v T ~ , ,T , which is half of (A.16A.16).To prove the second half for k ě
1, recall first that v T “ π ,k ´ T v h . A triangle inequality and (11a11a)with (cid:96) “ k ´ p “ α “ r “ s “ } ∇ v T } L p T q d À } ∇ v h } L p T q d . (A.18)Then, write } v F ´ v T } L p F q “ } π ,kF p v h ´ π ,k ´ T v h q} L p F q ď } v h ´ π ,k ´ T v h } L p F q ď h T } ∇ v h } L p T q d , where (11b11b) was used with (cid:96) “ k ´ p “ α “ r “ s “
1. Combining this estimate with (A.18A.18)concludes the proof of (A.16A.16). The above argument can be extended along similar lines to the case k “ v T is not an L -orthogonal projection of v h (see (1414)).To prove (A.17A.17), notice that if π ,k ` T v h “ k ` T v T “ r k ` T I k,lT v h “ π ,k ` T v h “
0. Hence, (A.16A.16)yields } ∇ v h } L p T q d » ~ v T ´ I k,lT r k ` T v T ~ , ,T , and the proof is complete by invoking (A.15A.15). Appendix A.3. Main results
We are now ready to prove the main results stated in Section 3.53.5.
Appendix A.3.1. Properties of Discontinuous Skeletal Gradient DiscretisationsProof of Theorem 99.
We use a polytopal toolbox in the spirit of [3131] and [2828, Section 7.2]. Let X T h , – (cid:32) w h “ pp w T q T P T h , p w F q F P F h q : w T P R , w F P R , w F “ @ F P F b h ( “ U , h, and Π T h : X T h , Ñ L p p Ω q , ∇ T h : X T h , Ñ L p p Ω q d be defined by, for all w h P X T h , and all T P T h , p Π T h w h q | T “ w T and p ∇ T h w h q | T “ | T | ÿ F P F T | F | w F n T F .
28y [2828, Corollary 7.12], the coercivity, limit-conformity, consistency, and compactness of p D h q h P H followif we find a mapping Φ : U k,lh, Ñ X T h , (“control” of D h ) such that, recalling the definition (A.11A.11) of ~¨~ ,p,h and setting } Φ } D h , T h – max v h P U k,lh, zt h u ~ Φ p v h q~ ,p,h } ∇ D h v h } L p p Ω q d , (A.19a) ω Π p D h , T h , Φ q – max v h P U k,lh, zt h u } Π D h v h ´ Π T h Φ p v h q} L p p Ω q } ∇ D h v h } L p p Ω q d , (A.19b) ω ∇ p D h , T h , Φ q – max v h P U k,lh, zt h u `ř T P T h | T | ´ p ˇˇş T p ∇ D h v h ´ ∇ T h Φ p v h qq ˇˇ p ˘ p } ∇ D h v h } L p p Ω q d , (A.19c)we have } Φ } D h , T h À h Ñ , ω Π p D h , T h , Φ q Ñ ω ∇ p D h , T h , Φ q Ñ , (A.20)where a À b means a ď Cb with real number C ą h , but possibly depending on d , (cid:37) , k , l , and k S .Let Φ be defined the following way. For all v h “ pp v T q T P T h , p v F q F P F h q P U k,lh , we let Φ p v h q – v h “ ` p v T q T P T h , p v F q F P F h ˘ P X T h , be such that v T “ π , T v T for all T P T h and v F “ π , F v F for all F P F h .Properties (A.20A.20) follow if we establish that, for all v h P U k,lh, , ~ v h ~ ,p,h À } ∇ D h v h } L p p Ω q d , (A.21a) } Π D h v h ´ Π T h v h } L p p T q À h } ∇ D h v h } L p p Ω q d , (A.21b)and, for all T P T h and all v T P U k,lT , p G T v T , η q T “ ÿ F P F T p v F , η ¨ n T F q F @ η P P p T q d . (A.21c)Indeed, (A.21aA.21a) gives a bound on } Φ } D h , T h , (A.21bA.21b) gives an O p h q estimate on ω Π p D h , T h , Φ q , and (A.21cA.21c)shows that ω ∇ p D h , T h , Φ q “ Proof of (A.21aA.21a) . By the definition (A.11A.11) of the ~¨~ ,p,h -seminorm along with that of v h , we havethat ~ v h ~ p ,p,h “ ÿ T P T h ÿ F P F T h ´ pF } v F ´ v T } pL p p F q . (A.22)Let now a mesh element T P T h be fixed and observe that, for all F P F T , } v F ´ v T } L p p F q “ } π , F p v F ´ v T q} L p p F q ď } v F ´ v T } L p p F q ď } v F ´ v T } L p p F q ` } v T ´ v T } L p p F q À } v F ´ v T } L p p F q ` h ´ p T } v T ´ v T } L p p T q À } v F ´ v T } L p p F q ` h ´ p T } ∇ v T } L p p T q d , (A.23)where we have used the L p -boundedness of the L -orthogonal projector (see [1616, Lemma 3.2]) in the firstline, the triangle inequality in the second line, the discrete L p -trace inequality (A.6A.6) in the third line, anda local Poincar´e–Wirtinger inequality which can be inferred from (11a11a) with α “ (cid:96) “ r “ s “ p th power of (A.23A.23), multiplying by h ´ pF » h ´ pT and summing over F P F T and T P T h leads to ~ v h ~ ,p,h À ~ v h ~ ,p,h . Estimate (A.21aA.21a) follows by using (A.13A.13).29ii) Proof of (A.21bA.21b) . Using a local Poincar´e–Wirtinger inequality as above we infer, for all T P T h , } v T ´ v T } L p p T q À h T } ∇ v T } L p p T q d . Taking the p th power of this inequality, summing over T P T h , andusing h T ď h and the uniform norm equivalence (A.13A.13) to bound the right-hand side, (A.21bA.21b) follows.(iii) Proof of (A.21cA.21c) . Let η P P p T q d Ă P k p T q d . Since ∇ ¨ η “
0, using the orthogonality property (S2)(S2) followed by the definition (1717) of G kT we infer that p G T v T , η q T “ ř F P F T p v F , η ¨ n T F q F . Equa-tion (A.21cA.21c) then follows by noticing that, η ¨ n T F being constant on F , p v F , η ¨ n T F q F “ p π , F v F , η ¨ n T F q F “p v F , η ¨ n T F q F .The GD-consistency follows from Proposition 1010 (proved below), and from [2828, Lemma 2.17] whichshows that the consistency holds provided that S D h p φ q Ñ φ in a dense subset of W ,p p Ω q . Remark
23 (Condition (A.21aA.21a)) . In [2828], a slightly different norm is considered in the argument of themaximum in (A.19aA.19a). The original expression is obtained replacing ~ v h ~ p ,p,h by ÿ T P T h ÿ F P F T d ´ pT F } v F ´ v T } pL p p F q , where the only difference with respect to (A.22A.22) is that the role of the local length scale is played by d T F ,the orthogonal distance between a point x T inside T and the face F , instead of h F . If, for all T P T h , wechoose the point x T such that condition (3434) is verified, it can easily be proved that d T F » h F , and thetwo norms are uniformly equivalent. Appendix A.3.2. Estimates on S D and W D Proof of Proposition 1010. (i)
Estimates on the addends in S D h . Take φ P W ,p p Ω q X W l ` ,p p T h q and let v h “ I k,lh φ P U k,lh . For all T P T h , if l ě p Π D h v h q | T “ v T “ π ,lT φ on T so the approximationestimate (11a11a) applied with α “ (cid:96) “ l , s “ (cid:96) ` r “ } Π D h v h ´ φ } L p p T q À h l ` T | φ | W l ` ,p p T q . (A.24)On the other hand, if l “ ´
1, the specific choice (1414) of v T yields } Π D h v h ´ φ } L p p T q À h T | φ | W ,p p T q . (A.25)Combining (A.24A.24) and (A.25A.25), taking the p th power, summing over T P T h , and taking the p th root ofthe resulting inequality gives (29a29a).For a fixed mesh element T P T h , use the definition (2525) of ∇ D h , the commutativity property (1818) of G kT , the approximation property (11a11a) of π ,kT with s “ k ` r “
0, and the consistency (A.2A.2) of S T to obtain } ∇ D h v h ´ ∇ φ } L p p T q d ď } G kT v T ´ ∇ φ } L p p T q d ` } S T v T } L p p T q d À h k ` T | ∇ φ | W k ` ,p p T q d ` h k ` T | φ | W k ` ,p p T q À h k ` T | φ | W k ` ,p p T q . The estimate (29b29b) follows taking the p th power, summing over T P T h , and taking the p th root of theresulting inequality.If l P t k, k ` u or l “ ´ k “ S D h p φ q is an immediateconsequence of (2929). Consider now l “ k ´ k ě
1. An easy modification of the proof above showsthat, for all φ P W ,p p Ω q X W k ` ,p p T h q , } ∇ D h v h ´ ∇ φ } L p p Ω q d À h k } φ } W k ` ,p p T h q . Then (3030) follows fromthis modified version of (29b29b) and from (29a29a).(ii) Estimate on W D h p ψ q . For all v h P X D h , , ż Ω ∇ D h v h p x q¨ ψ p x q d x “ ÿ T P T h p ∇ D h v h , ψ ´ π ,kT ψ q T ` ÿ T P T h p ∇ D h v h , π ,kT ψ q T — ÿ T P T h A T ` ÿ T P T h B T . (A.26)30he approximation property (11a11a) of π ,kT with s “ k ` r “ p instead of p yields ÿ T P T h | A T | À ÿ T P T h h k ` T } ∇ D h v h } L p p T q d } ψ } W k ` ,p p T q d ď h k ` ˜ ÿ T P T h } ∇ D h v h } pL p p T q d ¸ { p ˜ ÿ T P T h } ψ } p W k ` ,p p T q d ¸ { p “ h k ` } ∇ D h v h } L p p Ω q d } ψ } W k ` ,p p T h q d . (A.27)By definitions (2525) and (1717) of ∇ D h and G kT , and by the orthogonality property (S2)(S2) of S T , B T “ p G kT v T , π ,kT ψ q T “ ´ p v T , ∇ ¨ π ,kT ψ q T ` ÿ F P F T p v F , π ,kT ψ ¨ n T F q F “ ´ p v T , ∇ ¨p π ,kT ψ ´ ψ qq T ` ÿ F P F T p v F , p π ,kT ψ ´ ψ q¨ n T F q F ´ p v T , ∇ ¨ ψ q T ` ÿ F P F T p v F , ψ ¨ n T F q F “ p ∇ v T , p π ,kT ψ ´ ψ qq T ` ÿ F P F T p v F ´ v T , p π ,kT ψ ´ ψ q¨ n T F q F ´ p Π D h v h , ∇ ¨ ψ q T ` ÿ F P F T p v F , ψ ¨ n T F q F — B T, ´ p Π D h v h , ∇ ¨ ψ q T ` ÿ F P F T p v F , ψ ¨ n T F q F , (A.28)where we used an integration-by-parts and the definition (1616) of Π D h in the penultimate line. For anyinterface F with T , T as neighbouring mesh elements, since ψ P W p p div; Ω q we have ψ ¨ n T F ` ψ ¨ n T F “ F . Moreover, v F “ F is a boundary face. Hence ÿ T P T h ÿ F P F T p v F , ψ ¨ n T F q F “ ÿ F P F i h p v F , ψ ¨ n T F ` ψ ¨ n T F q F ` ÿ F P F b h p v F , ψ ¨ n T F q F “ . Summing (A.28A.28) over T P T h and using the previous relation leads to ÿ T P T h B T “ ÿ T P T h B T, ´ ż Ω Π D h v h p x q ∇ ¨ ψ p x q d x . (A.29)Recalling the definition (1010) of π ,kT , it is readily inferred that the first term in B T, is zero since ∇ v T P ∇ P l p T q Ă P k p T q d . Moreover, using again the approximation properties (1111) of π ,kT with r “ s “ k ` p instead of p , we can write | B T, | À ÿ F P F T } v F ´ v T } L p p F q } π ,kT ψ ´ ψ } L p p F q d À ÿ F P F T h k ` ´ p T } v F ´ v T } L p p F q | ψ | W k ` ,p p T q d À h k ` | ψ | W k ` ,p p T q d ˜ ÿ F P F T h p ´ F } v F ´ v T } L p p F q ¸ , where we used h F ď h T in the last line. Sum over T P T h and invoke H¨older’s inequality, the property31ard p F T q À
1, and the norm equivalence (A.13A.13) to deduce ÿ T P T h | B T, | À h k ` | ψ | W k ` ,p p T h q d ˜ ÿ T P T h ÿ F P F T h ´ pF } v F ´ v T } pL p p F q ¸ { p À h k ` } ψ } W k ` ,p p T h q d } ∇ D h v h } L p p Ω q d . (A.30)Finally, using (A.27A.27), (A.29A.29) and (A.30A.30) in (A.26A.26), we get ˇˇˇˇż Ω ∇ D h v h p x q¨ ψ p x q d x ` ż Ω Π D h v h p x q ∇ ¨ ψ p x q d x ˇˇˇˇ À h k ` } ψ } W k ` ,p p T h q d } ∇ D h v h } L p p Ω q d . The estimate (3131) follows immediately.
Remark
24 (Choice of the gradient reconstruction) . An inspection of the above proof shows that using (4343)in place of (2525) can lead to significant losses in the order of convergence for W D h p ψ q (while the convergenceexpressed by (99) still holds true). As a matter of fact, with this choice one would have to replacethroughout the proof π ,kT ψ by the L -orthogonal projection of ψ on ∇ P k ` p T q . The latter quantity hasoptimal approximation properties only if either k “ ∇ P p T q “ P p T q d ) or there exists w P L p p Ω q such that ψ | T “ ∇ w | T for all T P T h . Recalling Theorem 33 with p “ σ p x , u, ∇ u q “ ∇ u , we seethat for the Poisson equation, W D h is applied to ψ “ ∇ u . In this case, the gradient reconstruction (4343)leads to optimal convergence rates. However, there is a real loss of estimate for more general problemsfor which error estimates are written in terms of W D h p Λ ∇ u q (for anisotropic linear diffusion, see [2828,Theorem 2.29]) or W D h p| ∇ u | p ´ ∇ u q (for the p -Laplace equation, see Theorem 33). ReferencesReferences [1] Arnold, D. N., 1982. An interior penalty finite element method with discontinuous elements. SIAMJ. Numer. Anal. 19, 742–760.[2] Ayuso de Dios, B., Lipnikov, K., Manzini, G., 2016. The nonconforming virtual element method.ESAIM: Math. Model Numer. Anal. 50 (3), 879–904.[3] Beir˜ao da Veiga, L., Brezzi, F., Cangiani, A., Manzini, G., Marini, L. D., Russo, A., 2013. Basicprinciples of virtual element methods. Math. Models Methods Appl. Sci. 199 (23), 199–214.[4] Beir˜ao da Veiga, L., Lipnikov, K., Manzini, G., 2014. The mimetic finite difference method for ellipticproblems. Vol. 11 of MS&A. Modeling, Simulation and Applications. Springer, Cham.URL http://dx.doi.org/10.1007/978-3-319-02663-3http://dx.doi.org/10.1007/978-3-319-02663-3 [5] Boffi, D., Brezzi, F., Fortin, M., 2013. Mixed finite element methods and applications. Vol. 44 ofSpringer Series in Computational Mathematics. Springer, Heidelberg.URL http://dx.doi.org/10.1007/978-3-642-36519-5http://dx.doi.org/10.1007/978-3-642-36519-5 [6] Boffi, D., Di Pietro, D. A., 2017. Unified formulation and analysis of mixed and primal discontinuousskeletal methods on polytopal meshes. ESAIM: Math. Model Numer. Anal.Published online.URL http://dx.doi.org/10.1051/m2an/2017036http://dx.doi.org/10.1051/m2an/2017036 [7] Bonelle, J., Di Pietro, D. A., Ern, A., 2015. Low-order reconstruction operators on polyhedralmeshes: application to compatible discrete operator schemes. Computer Aided Geometric Design35–36, 27–41.URL http://dx.doi.org/10.1016/j.cagd.2015.03.015http://dx.doi.org/10.1016/j.cagd.2015.03.015 http://dx.doi.org/10.1137/040613950http://dx.doi.org/10.1137/040613950 [9] Brezzi, F., Lipnikov, K., Simoncini, V., 2005. A family of mimetic finite difference methods onpolygonal and polyhedral meshes. Mathematical Models and Methods in Applied Sciences 15 (10).[10] Cangiani, A., Gyrya, V., Manzini, G., 2016. The nonconforming virtual element method for theStokes equations. SIAM Journal on Numerical Analysis 54 (6), 3411–3435.[11] Cangiani, A., Manzini, G., Sutton, O. J., 2017. Conforming and nonconforming virtual elementmethods for elliptic problems. IMA J. Numer. Anal. 37 (3), 1317–1354.[12] Castillo, P., Cockburn, B., Perugia, I., Sch¨otzau, D., 2000. An a priori error analysis of the localdiscontinuous Galerkin method for elliptic problems. SIAM J. Numer. Anal. 38, 1676–1706.[13] Cockburn, B., Di Pietro, D. A., Ern, A., 2016. Bridging the Hybrid High-Order and HybridizableDiscontinuous Galerkin methods. ESAIM: Math. Model. Numer. Anal. 50 (3), 635–650.URL http://dx.doi.org/10.1051/m2an/2015051http://dx.doi.org/10.1051/m2an/2015051 [14] Cockburn, B., Gopalakrishnan, J., Lazarov, R., 2009. Unified hybridization of discontinuousGalerkin, mixed, and continuous Galerkin methods for second order elliptic problems. SIAM J.Numer. Anal. 47 (2), 1319–1365.URL http://dx.doi.org/10.1137/070706616http://dx.doi.org/10.1137/070706616 [15] Crouzeix, M., Raviart, P.-A., 1973. Conforming and nonconforming finite element methods for solv-ing the stationary Stokes equations. I. Rev. Fran¸caise Automat. Informat. Recherche Op´erationnelleS´er. Rouge 7 (R-3), 33–75.[16] Di Pietro, D. A., Droniou, J., 2017. A Hybrid High-Order method for Leray–Lions elliptic equationson general meshes. Math. Comp. 86 (307), 2159–2191.URL http://dx.doi.org/10.1090/mcom/3180http://dx.doi.org/10.1090/mcom/3180 [17] Di Pietro, D. A., Droniou, J., 2017. W s,p -approximation properties of elliptic projectors on polyno-mial spaces, with application to the error analysis of a Hybrid High-Order discretisation of Leray–Lions problems. Math. Models Methods Appl. Sci. 27 (5), 879–908.URL http://dx.doi.org/10.1142/S0218202517500191http://dx.doi.org/10.1142/S0218202517500191 [18] Di Pietro, D. A., Droniou, J., Ern, A., 2015. A discontinuous-skeletal method for advection-diffusion-reaction on general meshes. SIAM J. Numer. Anal. 53 (5), 2135–2157.URL https://dx.doi.org/10.1137/140993971https://dx.doi.org/10.1137/140993971 [19] Di Pietro, D. A., Ern, A., 2012. Mathematical aspects of discontinuous Galerkin methods. Vol. 69of Math´ematiques & Applications. Springer-Verlag, Berlin.[20] Di Pietro, D. A., Ern, A., 2015. A hybrid high-order locking-free method for linear elasticity ongeneral meshes. Comput. Meth. Appl. Mech. Engrg. 283, 1–21.URL http://dx.doi.org/10.1016/j.cma.2014.09.009http://dx.doi.org/10.1016/j.cma.2014.09.009 [21] Di Pietro, D. A., Ern, A., 2015. Hybrid high-order methods for variable-diffusion problems on generalmeshes. C. R. Acad. Sci. Paris, Ser. I 353, 31–34.[22] Di Pietro, D. A., Ern, A., Lemaire, S., 2014. An arbitrary-order and compact-stencil discretizationof diffusion on general meshes based on local reconstruction operators. Comput. Meth. Appl. Math.14 (4), 461–472.URL http://dx.doi.org/10.1515/cmam-2014-0018http://dx.doi.org/10.1515/cmam-2014-0018 http://dx.doi.org/10.1007/s10915-017-0512-xhttp://dx.doi.org/10.1007/s10915-017-0512-x [24] Di Pietro, D. A., Lemaire, S., 2015. An extension of the Crouzeix–Raviart space to general mesheswith application to quasi-incompressible linear elasticity and Stokes flow. Math. Comp. 84 (291),1–31.URL http://dx.doi.org/10.1090/S0025-5718-2014-02861-5http://dx.doi.org/10.1090/S0025-5718-2014-02861-5 [25] Di Pietro, D. A., Tittarelli, R., 2017. Numerical methods for PDEs. Lectures from the fall 2016thematic quarter at Institut Henri Poincar´e. SEMA SIMAI series. Springer, Ch. An introduction toHybrid High-Order methods, accepted for publication. Preprint arXiv: 1703.051361703.05136 [math.NA].[26] Diaz, J. I., de Thelin, F., 1994. On a nonlinear parabolic problem arising in some models related toturbulent flows. SIAM J. Math. Anal. 25 (4), 1085–1111.[27] Droniou, J., Eymard, R., 2006. A mixed finite volume scheme for anisotropic diffusion problems onany grid. Numer. Math. 105 (1), 35–71.URL http://dx.doi.org/10.1007/s00211-006-0034-1http://dx.doi.org/10.1007/s00211-006-0034-1 [28] Droniou, J., Eymard, R., Gallou¨et, T., Guichard, C., Herbin, R., 2017. The gradient discretisationmethod: A framework for the discretisation and numerical analysis of linear and nonlinear ellipticand parabolic problems. Maths & Applications. Springer, to appear. Preprint hal-01382358hal-01382358, version4.URL https://hal.archives-ouvertes.fr/hal-01382358https://hal.archives-ouvertes.fr/hal-01382358 [29] Droniou, J., Eymard, R., Gallou¨et, T., Herbin, R., 2010. A unified approach to mimetic finitedifference, hybrid finite volume and mixed finite volume methods. Math. Models Methods Appl. Sci.20 (2), 1–31.[30] Droniou, J., Eymard, R., Gallou¨et, T., Herbin, R., 2013. Gradient schemes: a generic framework forthe discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations. Math. ModelsMethods Appl. Sci. 23 (13), 2395–2432.URL http://dx.doi.org/10.1142/S0218202513500358http://dx.doi.org/10.1142/S0218202513500358 [31] Droniou, J., Eymard, R., Herbin, R., 2016. Gradient schemes: Generic tools for the numericalanalysis of diffusion equations. ESAIM Mathematical Modelling and Numerical Analysis 50 (3),749–781.[32] Droniou, J., Lamichhane, B. P., 2015. Gradient schemes for linear and non-linear elasticity equations.Numer. Math. 129 (2), 251–277.URL http://dx.doi.org/10.1007/s00211-014-0636-yhttp://dx.doi.org/10.1007/s00211-014-0636-y [33] Ern, A., Guermond, J.-L., 2004. Theory and Practice of Finite Elements. Vol. 159 of Applied Math-ematical Sciences. Springer-Verlag, New York, NY.[34] Eymard, R., Gallou¨et, T., Herbin, R., 2010. Discretization of heterogeneous and anisotropic diffu-sion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybridinterfaces. IMA J. Numer. Anal. 30 (4), 1009–1043.URL http://dx.doi.org/10.1093/imanum/drn084http://dx.doi.org/10.1093/imanum/drn084 [35] Eymard, R., Guichard, C., 2017. The discontinuous Galerkin gradient discretisation. Preprinthal-01535147hal-01535147.URL https://hal.archives-ouvertes.fr/hal-01535147https://hal.archives-ouvertes.fr/hal-01535147 http://dx.doi.org/10.1016/j.jcp.2014.04.021http://dx.doi.org/10.1016/j.jcp.2014.04.021 [42] Lorentz, J., Neilan, M., Smears, I., 2015. Stable Discontinuous Galerkin FEM without penaltyparameter. In: Numerical Mathematics and Advanced Applications ENUMATH 2015. Vol. 112 ofLecture Notes in Computational Science and Engineering. Springer.[43] N´ed´elec, J.-C., 1980. Mixed finite elements in R . Numer. Math. 35 (3), 315–341.URL http://dx.doi.org/10.1007/BF01396415http://dx.doi.org/10.1007/BF01396415http://dx.doi.org/10.1007/BF01396415http://dx.doi.org/10.1007/BF01396415