Skip to content
Publicly Available Published by De Gruyter June 17, 2017

An hp-Hybrid High-Order Method for Variable Diffusion on General Meshes

Joubine Aghili, Daniele A. Di Pietro and Berardo Ruffini


In this work, we introduce and analyze an hp-hybrid high-order (hp-HHO) method for a variable diffusion problem. The proposed method is valid in arbitrary space dimension and for fairly general polytopal meshes. Variable approximation degrees are also supported. We prove hp-convergence estimates for both the energy- and L2-norms of the error, which are the first of this kind for Hybrid High-Order methods. These results hinge on a novel hp-approximation lemma valid for general polytopal elements in arbitrary space dimension. The estimates are additionally fully robust with respect to the heterogeneity of the diffusion coefficient, and show only a mild dependence on the square root of the local anisotropy, improving previous results for HHO methods. The expected exponential convergence behavior is numerically demonstrated on a variety of meshes for both isotropic and strongly anisotropic diffusion problems.

MSC 2010: 65N30; 65N08; 76R50

1 Introduction

Over the last few years, discretization technologies have appeared that support arbitrary approximation orders on general polytopal meshes. In this work, we focus on a particular instance of such technologies, the so-called Hybrid High-Order (HHO) methods originally introduced in [18, 17]; see also [21] for an introduction. So far, the literature on HHO methods has focused on the h-version of the method with uniform polynomial degree. Our goal is to provide a first example of variable-degree hp-HHO method and carry out a full hp-convergence analysis valid for fairly general meshes and arbitrary space dimension. Let Ωd, d1, denote a bounded connected polytopal domain with Lipschitz boundary. We consider the variable diffusion model problem: Find u:Ω such that

(1.1){-(𝜿u)=fin Ω,u=0on Ω,

where 𝜿 is a uniformly positive, symmetric, tensor-valued field on Ω, while fL2(Ω) denotes a volumetric source. For the sake of simplicity, we assume that 𝜿 is piecewise constant on a partition PΩ of Ω into polytopes. The weak formulation of problem (1.1) reads: Find uU:=H01(Ω) such that

(1.2)(𝜿u,v)=(f,v)for all vU,

where we have used the notation (,) for the usual inner products of both L2(Ω) and L2(Ω)d. Here, the scalar-valued field u represents a potential, and the vector-valued field 𝜿u the corresponding flux.

For a given polytopal mesh 𝒯h={T} of Ω, the hp-HHO discretization of problem (1.2) is based on two sets of degrees of freedom (DOFs):

  1. skeletal DOFs, consisting in (d-1)-variate polynomials of total degree pF0 on each mesh face F, and

  2. elemental DOFs, consisting in d-variate polynomials of degree pT on each mesh element T, where pT denotes the lowest degree of skeletal DOFs on the boundary of T.

Skeletal DOFs are globally coupled, and can be alternatively interpreted as traces of the potential on the mesh faces or as Lagrange multipliers enforcing the continuity of the normal flux at the discrete level; cf. [1, 12] for further insight. Elemental DOFs, on the other hand, are auxiliary DOFs that can be locally eliminated by static condensation, as detailed in [12, Section 2.4] for the case where pF=p for all mesh faces F.

Two key ingredients are devised locally from skeletal and elemental DOFs attached to each mesh element T:

  1. a reconstruction of the potential of degree (pT+1) (i.e., one degree higher than elemental DOFs in T) obtained solving a small Neumann problem, and

  2. a stabilization term penalizing face-based residuals and polynomially consistent up to degree (pT+1).

The local contributions obtained from these two ingredients are then assembled following a standard, finite element-like procedure. The resulting discretization has several appealing features, the most prominent of which are summarized hereafter:

  1. It is valid for fairly general polytopal meshes.

  2. The construction is dimension-independent, which can significantly ease the practical implementation.

  3. It enables the local adaptation of the approximation order, a highly desirable feature when combined with a regularity estimator (whose development will be addressed in a separate work).

  4. It exhibits only a moderate dependence on the diffusion coefficient 𝜿.

  5. It has a moderate computational cost thanks to the possibility of eliminating elemental DOFs locally via static condensation.

  6. Parallel implementations can be simplified by the fact that processes communicate via skeletal unknowns only.

The seminal works on the p- and hp-conforming finite element method on standard meshes date back to the early 80s; cf. [7, 6, 5]. Starting from the late 90s, nonconforming methods on standard meshes supporting arbitrary-order have received a fair amount of attention; a (by far) nonexhaustive list of contributions focusing on scalar diffusive problems similar to the one considered here includes [29, 10, 28, 24, 31]. The possibility of refining both in h and in p on general meshes is, on the other hand, a much more recent research topic. We cite, in particular, hp-composite [3, 25] and polyhedral [9] discontinuous Galerkin methods, and the two-dimensional virtual element method proposed in [8].

The main novelty of the paper are hp-energy- and L2-estimates for HHO methods summarized in Section 3.2. These are the first results of this kind for HHO methods, and among the first for discontinuous skeletal methods in general (a prominent example of discontinuous skeletal methods are the hybridizable discontinuous Galerkin methods of [13]; cf. [12] for a precise study of their relation with HHO methods).

The cornerstone of the analysis is a novel extension of the classical Babuška–Suri hp-approximation results to regular mesh sequences in the sense of [16, Chapter 1] and arbitrary space dimension d1; see Lemma 2.3 below. Similar results had been derived in [8] for d=2 and, under assumptions on the mesh different from the ones proposed here, in [9] for d{2,3}. A relevant difference with respect to classical approximation results based on the Bramble–Hilbert lemma (see, e.g., [23, Theorem 1.103]) is that we cannot rely here on the notion of reference element, since general polyhedral shapes are allowed. As a result, it turns out that a key point consists in showing that the regularity assumptions adapted to polyhedral meshes imply uniform bounds for the Lipschitz constant of mesh elements.

The resulting energy-norm estimate confirms the characteristic h-superconvergence behavior of HHO methods, whereas we have a more standard scaling as (pT+1)-pT with respect to the polynomial degree pT of elemental DOFs. This scaling is analogous to the best available results for discontinuous Galerkin (dG) methods on rectangular meshes based on polynomials of degree pT, cf. [24] (on more general meshes, the scaling for the symmetric interior penalty dG method is pT-(pT-1/2), half a power less than for the hp-HHO method studied here). Classically, when elliptic regularity holds, the h-convergence order can be increased by 1 for the L2-norm. In our error estimates, the dependence on the diffusion coefficient is carefully tracked, showing full robustness with respect to its heterogeneity and only a moderate dependence with a power of 12 on its local anisotropy when the error in the energy-norm is considered. These results improve the ones in [19] where, considering a different norm for the error, a power of 1 was found. Numerical experiments confirm the expected exponentially convergent behavior for both isotropic and strongly anisotropic diffusion coefficients on a variety of two-dimensional meshes.

The rest of the paper is organized as follows. In Section 2 we introduce the notation and state the basic results required in the analysis including, in particular, Lemma 2.3 (whose proof is detailed in Appendix A). In Section 3 we formulate the hp-HHO method, state our main results, and provide some numerical examples. The proofs of the main results, preceded by the required preparatory material, are collected in Section 4.

2 Setting

In this section we introduce the notation and state the basic results required in the analysis.

2.1 Mesh and Notation

Let *+ denote a countable set of meshsizes having 0 as its unique accumulation point. We consider mesh sequences (𝒯h)h where, for all h, 𝒯h={T} is a finite collection of nonempty disjoint open polytopal elements such that Ω¯=T𝒯hT¯ and h=maxT𝒯hhT (hT stands for the diameter of T). A hyperplanar closed connected subset F of Ω¯ is called a face if it has positive (d-1)-dimensional measure and

  1. either there exist distinct T1,T2𝒯h such that F=T1T2 (and F is an interface)

  2. or there exists T𝒯h such that F=TΩ (and F is a boundary face).

Interfaces are collected in hi, boundary faces in hb, and we let h:=hihb. It is assumed that the set of faces h partitions the mesh skeleton in the sense that distinct faces have disjoint interiors and that T𝒯hT=FhF. For all T𝒯h, the set T:={FhFT} collects the faces lying on the boundary of T and, for all FT, we denote by 𝒏TF the normal vector to F pointing out of T.

The following assumptions on the mesh will be kept throughout the exposition.

Assumption 2.1 (Admissible Mesh Sequence).

We assume that (𝒯h)h is admissible in the sense of [16, Chapter 1], i.e., for all h, 𝒯h admits a matching simplicial submesh 𝔗h and there exists a real number ϱ>0 (the mesh regularity parameter) independent of h such that the following conditions hold:

  1. ϱhSrS for all h and all simplices S𝔗h of diameter hS and inradius rS.

  2. ϱhThS for all h, all T𝒯h, and all S𝔗h such that ST.

  3. Every mesh element T𝒯h is star-shaped with respect to every point of a ball of radius ϱhT.

Assumption 2.2 (Compliant Mesh Sequence).

We assume that the mesh sequence is compliant with the partition PΩ on which the diffusion tensor 𝜿 is piecewise constant, so that jumps only occur at interfaces and, for all T𝒯h, we set


In what follows, for all T𝒯h, κ¯T and κ¯T denote the largest and smallest eigenvalue of 𝜿T, respectively, and λ𝜿,T:=κ¯T/κ¯T the local anisotropy ratio.

2.2 Basic Results

For a subset X of Ω and an index q, Hq(X) denotes the Hilbert space of functions which are in L2(X) together with their weak derivatives of order q, equipped with the usual inner product (,)q,X and associated norm q,X. When q=0, we recover the Lebesgue space L2(X), and the subscript 0 is omitted from both the inner product and the norm. The subscript X is also omitted when X=Ω. For a given integer l0, we denote by l(X) the space spanned by the restriction to X of d-variate polynomials of degree l. For further use, we also introduce the L2-projector πXl:L1(X)l(X) such that, for all wL1(X),

(2.1)(πXlw-w,v)X=0for all vl(X).

We recall hereafter a few known results on admissible mesh sequences and refer to [16, Chapter 1] and [14, 15] for a more comprehensive collection. By [16, Lemma 1.41], there exists an integer N(d+1) (possibly depending on d and ϱ) such that the maximum number of faces of one mesh element is bounded,


The following multiplicative trace inequality, valid for all h, all T𝒯h, and all vH1(T), is proved in [16, Lemma 1.49]:


where C only depends on d and ϱ. We also note the following local Poincaré’s inequality valid for all T𝒯h and all vH1(T) such that (v,1)T=0:


where CP=π-1 when T is convex, while it can be estimated in terms of ϱ for nonconvex elements (cf., e.g., [34]).

The following functional analysis results lie at the heart of the hp-analysis carried out in Section 4.

Lemma 2.3 (Approximation).

There is a real number C>0 (possibly depending on d and ϱ) such that, for all hH, all TTh, all integers l1, all s0, and all vHs+1(T), there exists a polynomial ΠTlvPl(T) satisfying, for all 0qs+1,



See Appendix A. ∎

Lemma 2.4 (Discrete Trace Inequality).

There is a real number C>0 (possibly depending on d and ϱ) such that, for all hH, all TTh, all integers l1, and all vPl(T), it holds



When all meshes in the sequence (𝒯h)h are simplicial and conforming, the proof of (2.6) can be found in [30, Theorem 4.76] for d=2; for d2 the proof is analogous. The extension to admissible mesh sequences in the sense of Assumption 2.1 can be done following the reasoning in [16, Lemma 1.46]. ∎

3 Discretization

In this section, we formulate the hp-HHO method, state our main results, and provide some numerical examples.

3.1 The hp-HHO Method

We present in this section an extension of the classical HHO method of [18] accounting for variable polynomial degrees. Let a vector p¯h=(pF)Fhh of skeletal polynomial degrees be given. For all T𝒯h, we denote by p¯T=(pF)FT the restriction of p¯h to T, and we introduce the following local space of DOFs:


The standard HHO notation v¯T=(vT,(vF)FT) is used for a generic element of U¯Tp¯T. We define the local potential reconstruction operator


such that, for all v¯TU¯Tp¯T and all wpT+1(T),




Note that computing rTpT+1v¯T according to (3.2) requires to invert the 𝜿T-weighted stiffness matrix of k+1(T), which can be efficiently accomplished by a Cholesky solver.

We define on U¯Tp¯T×U¯Tp¯T the local bilinear form aT such that


where, for all FT, we let κF:=𝜿T𝒏TF𝒏TF and the face-based residual operator


is such that, for all v¯TU¯Tp¯T,


The first contribution in aT is in charge of consistency, whereas the second ensures stability by a least-square penalty of the face-based residuals δTFp¯T. This subtle form for δTFp¯T ensures that the residual vanishes when its argument is the interpolate of a function in pT+1(T), and is required for high-order h-convergence (a detailed motivation is provided in [18, Remark 6]).

The global space of DOFs and its subspace with strongly enforced boundary conditions are defined, respectively, as

U¯hp¯h:=(×T𝒯hpT(T))×(×FhpF(F)),U¯h,0p¯h:={v¯hU¯hp¯hvF0 for all Fhb}.

Note that interface DOFs in U¯hp¯h are single-valued, meaning that they match from one element to the adjacent one. We use the notation


for a generic DOF vector in U¯hp¯h and, for all T𝒯h, we denote by v¯TU¯Tp¯T its restriction to T. For further use, we also introduce the global interpolator


such that, for all vH1(Ω),


and denote by I¯Tp¯T its restriction to T𝒯h.

The hp-HHO discretization of problem (1.2) consists in seeking u¯hU¯h,0p¯h such that

(3.7)ah(u¯h,v¯h)=lh(v¯h)for all v¯hU¯h,0p¯h,

where the global bilinear form ah on U¯hp¯h×U¯hp¯h and the linear form lh on U¯hp¯h are assembled element-wise setting


Remark 3.1 (Static Condensation).

Using a standard static condensation procedure, it is possible to eliminate element-based DOFs locally and solve (3.7) by inverting a system in the skeletal unknowns only. For the sake of conciseness, we do not repeat the details here and refer instead to [12, Section 2.4]. Accounting for the strong enforcement of boundary conditions, the size of the system after static condensation is


Remark 3.2 (Finite Element Interpretation).

A finite element interpretation of the scheme (3.7) is possible following the extension proposed in [12, Remark 3] of the ideas originally developed in [4] in the context of nonconforming virtual element methods. For all Fhi, we denote by []F the usual jump operator (the sign is irrelevant), which we extend to boundary faces Fhb setting [φ]F:=φ. Let

𝔘h,0p¯h:={𝔳hL2(Ω)𝔳h|T𝔘Tp¯T for all T𝒯h and πFpF([𝔳T]F)=0 for all Fh},

where, for all T𝒯h, we have introduced the local space

𝔘Tp¯T:={𝔳TH1(T)(𝜿T𝔳T)pT(T) and 𝜿T𝔳T|F𝒏TFpF(F) for all FT}.

It can be proved that, for all T𝒯h,


is an isomorphism. Thus, the triplet (T,𝔘Tp¯T,I¯Tp¯T) defines a finite element in the sense of Ciarlet [11]. Additionally, problem (3.7) can be reformulated as the nonconforming finite element method: Find 𝔲h𝔘h,0p¯h such that

𝔞h(𝔲h,𝔳h)=lh(𝔳h)for all 𝔳h𝔘h,0p¯h,



and it can be proved that 𝔲h is the unique element of 𝔘h,0p¯h such that u¯h=I¯hp¯h𝔲h with u¯h unique solution to (3.7).

3.2 Main Results

We next state our main results. The proofs are postponed to Section 4. For all T𝒯h, we denote by a,T and ||s,T the seminorms defined on U¯Tp¯T by the bilinear forms aT and sT, respectively, and by a,h the seminorm defined by the bilinear form ah on U¯hp¯h. We also introduce the penalty seminorm ||s,h such that, for all v¯hU¯hp¯h,


Note that a,h is a norm on the subspace U¯h,0p¯h with strongly enforced boundary conditions (the arguments are analogous to those of [17, Proposition 5]). We will also need the global reconstruction operator


such that, for all v¯hU¯hp¯h,

(rhp¯hv¯h)=|TrTpT+1v¯Tfor all T𝒯h.

Finally, for the sake of conciseness, throughout the rest of the paper we denote by ab the inequality aCb with real number C>0 independent of h, p¯h, 𝜿 and, for local inequalities on a mesh element T𝒯h, also of T.

Our first estimate concerns the error measured in energy-like norms.

Theorem 3.3 (Energy Error Estimate).

Let uU and u¯hU¯h,0p¯h denote the unique solutions of problems (1.2) and (3.7), respectively, and set


Assuming that u|THpT+2(T) for all TTh, it holds


Consequently, denoting by h the broken gradient on Th, we have that



See Section 4.3. ∎

In (3.10) and (3.11), we observe the characteristic improved h-convergence of HHO methods (cf. [18]), whereas, in terms of p-convergence, we have a more standard scaling as (pT+1)-pT (i.e., half a power more than discontinuous Galerkin methods based on polynomials of degree pT, cf., e.g., [28]). In (3.11), we observe that the left-hand side has the same convergence rate (both in h and in p) as the interpolation error


as can be verified combining (4.2) and (4.4) below. Note that, in this case, the p-convergence is limited by the second term, which measures the discontinuity of the potential reconstruction at interfaces. An inspection of formulas (3.10) and (3.11) also shows that the method is fully robust with respect to the heterogeneity of the diffusion coefficient, while only a moderate dependence (with a power of 12) is observed with respect to its local anisotropy ratio.

For the sake of completeness, we also provide an estimate of the L2-error between the piecewise polynomial fields uh and u^h such that

uh:=|TuTandu^h:=|Tu^T=πTpTufor all T𝒯h.

To this end, we need elliptic regularity in the following form: For all gL2(Ω), the unique element zU such that

(3.12)(𝜿z,v)=(g,v)for all vU,

satisfies the a priori estimate


The following result is proved in Section 4.4.

Figure 1

Convergence with p-refinement (regular test case).

(a) ∥u¯h-u¯^h∥a,h\lVert\underline{u}_{h}-\widehat{\underline{u}}_{h}\rVert_{\mathrm{a},h} vs. NdofN_{\mathrm{dof}}


u¯h-u¯^ha,h vs. Ndof

(b) ∥uh-u^h∥\lVert u_{h}-\widehat{u}_{h}\rVert vs. NdofN_{\mathrm{dof}}


uh-u^h vs. Ndof

(c) ∥u¯h-u¯^h∥a,h\lVert\underline{u}_{h}-\widehat{\underline{u}}_{h}\rVert_{\mathrm{a},h} vs. p


u¯h-u¯^ha,h vs. p

(d) ∥uh-u^h∥\lVert u_{h}-\widehat{u}_{h}\rVert vs. p


uh-u^h vs. p

Figure 2

Convergence with p-refinement (Le Potier’s test case).

(a) ∥u¯h-u¯^h∥a,h\lVert\underline{u}_{h}-\widehat{\underline{u}}_{h}\rVert_{\mathrm{a},h} vs. NdofN_{\mathrm{dof}}


u¯h-u¯^ha,h vs. Ndof

(b) ∥uh-u^h∥\lVert u_{h}-\widehat{u}_{h}\rVert vs. NdofN_{\mathrm{dof}}


uh-u^h vs. Ndof

(c) ∥u¯h-u¯^h∥a,h\lVert\underline{u}_{h}-\widehat{\underline{u}}_{h}\rVert_{\mathrm{a},h} vs. p


u¯h-u¯^ha,h vs. p

(d) ∥uh-u^h∥\lVert u_{h}-\widehat{u}_{h}\rVert vs. p


uh-u^h vs. p

Theorem 3.4 (L2-Error Estimate).

Under the assumptions of Theorem 3.3, and further assuming elliptic regularity (3.13) and that fHpT+ΔT(T) for all TTh with ΔT=1 if pT=0 while ΔT=0 otherwise, it holds


with λ𝛋:=maxTThλ𝛋,T and κ¯:=maxTThκ¯T.

3.3 Numerical Examples

We close this section with some numerical examples. The h-convergence properties of the method (3.7) have been numerically investigated in [18, Section 4]. To illustrate its p-convergence properties, we solve on the unit square domain Ω=(0,1)2 the homogeneous Dirichlet problem with exact solution u=sin(πx1)sin(πx2) and right-hand side f chosen accordingly. We consider two values for the diffusion coefficients:


where 𝑰2 denotes the identity matrix of dimension 2, 𝒙¯:=-(0.1,0.1), and ϵ=110-2. The choice 𝜿=𝜿1 (“regular” test case) yields a homogeneous isotropic problem, while the choice 𝜿=𝜿2 (“Le Potier’s” test case [27]) corresponds to a highly anisotropic problem where the principal axes of the diffusion tensor vary at each point of the domain. Figures 1 and 2 depict the energy- and L2-errors as a function of the number of skeletal DOFs Ndof (cf. (3.8)) when pF=p for all Fh and p{0,,9} for the proposed choices for 𝜿 on the meshes of Figure 3. In all the cases, the expected exponentially convergent behavior is observed. Interestingly, the best performance in terms of error vs. Ndof is obtained for the Cartesian and Voronoi meshes. A comparison of the results for the two values of the diffusion coefficients allows to appreciate the robustness of the method with respect to anisotropy.

Figure 3

Meshes considered in the p-convergence test of Section 3.3. The triangular, Cartesian, refined, and staggered meshes originate from the FVCA5 benchmark [26]; the hexagonal mesh was originally introduced in [20]; the Voronoi mesh was obtained using the PolyMesher algorithm of [33].

(a) Triangular



(b) Cartesian



(c) Refined



(d) Staggered



(e) Hexagonal



(f) Voronoi



4 Convergence Analysis

In this section we prove the results stated in Section 3.2.

4.1 Consistency of the Potential Reconstruction

Preliminary to the convergence analysis is the study of the approximation properties of the potential reconstruction rTpT+1 defined by (3.2) when its argument is the interpolate of a regular function. Let a mesh element T𝒯h be fixed. For any integer l1, we define the oblique elliptic projector ϖ𝜿,Tl:H1(T)l(T) such that, for all vH1(T), (ϖ𝜿,Tlv-v,1)T=0 and it holds

(4.1)(𝜿T(ϖ𝜿,Tlv-v),w)T=0for all wl(T).

Proposition 4.1 (Characterization of rTpT+1I¯Tp¯T).

It holds, for all TTh,



For a generic vH1(T), letting v¯T=I¯Tp¯Tv in (3.2), we infer, for all wpT+1(T),


where we have used the fact that (𝜿Tw)pT-1(T)pT(T) and (𝜿Tw)|F𝒏TFpT(F)pF(F) (cf. definition (3.1) of pT) to cancel the projectors in the second line, and an integration by parts to conclude. ∎

We next study the approximation properties of ϖ𝜿,Tl, from which those of (rTpT+1I¯Tp¯T) follow in the light of Proposition 4.1.

Lemma 4.2 (Approximation Properties of ϖκ,Tl).

For all integers l1, all mesh elements TTh, all 0sl, and all vHs+1(T), it holds



By definition (4.1) of ϖ𝜿,Tl, it holds


Hence, using (2.5) with q=1, it is readily inferred that


To prove the second bound in (4.2), use the triangle inequality to infer


For the first term, the multiplicative trace inequality (2.3) combined with (2.5) (with q=1,2) gives


For the second term, we have


where we have used the discrete trace inequality (2.6) in the first line, the triangle inequality in the second line, the estimate (4.3) in the third, and the approximation result (2.5) with q=1 to conclude. To obtain the third bound in (4.2), after recalling that (v-ϖ𝜿,Tlv,1)T=0, we apply the local Poincaré inequality (2.4) to infer


where the conclusion follows from the first bound in (4.2). Finally, to obtain the last bound, we use the multiplicative trace inequality (2.3) to infer


and use the first and third bound in (4.2) to estimate the various terms. ∎

4.2 Consistency of the Stabilization Term

The consistency properties of the stabilization bilinear form sT defined by (3.4) are summarized in the following lemma.

Lemma 4.3 (Consistency of the Stabilization Term).

For all TTh, all 0qpT, and all vHq+2(T), it holds



Let T𝒯h and vHq+2(T) and set, for the sake of brevity,


(cf. Proposition 4.1). For all FT, recalling the definitions of the face residual δTFp¯T (cf. (3.5)) and of the local interpolator I¯Tp¯T (cf. (3.6)), together with the fact that pTpF by definition (3.1), we get


Using the triangle inequality and the L2(F)-stability of πFpF, we infer,


For the first term, the approximation properties (4.2) of ϖ𝜿,TpT+1 (with l=pT+1 and s=q+1) readily yield


For the second term, on the other hand, the discrete trace inequality (2.6) followed by the L2(T)-stability of πTpT and (4.2) (with l=pT+1 and s=q+1) gives


The bound (4.4) follows using (4.6)–(4.7) in the right-hand side of (4.5), squaring the resulting inequality, multiplying it by κF/hT, summing over FT, and using the bound (2.2) on card(T). ∎

4.3 Energy Error Estimate

Proof of Theorem 3.3.

We start by noting the following abstract error estimate:


with nonconformity error


To prove (4.8), it suffices to observe that


where we have used the definition of the a,h-norm in the first line, the linearity of ah in its first argument in the second line, and the discrete problem (3.7) in the third. The conclusion follows dividing both sides by u¯h-u¯^ha,h, using linearity, and passing to the supremum.

We next bound the nonconformity error Eh(v¯h) for a generic vector of DOFs v¯hU¯h,0p¯h. A preliminary step consists in finding a more appropriate rewriting for Eh(v¯h). Observing that f=-(𝜿u) a.e. in Ω, integrating by parts element-by-element, and using the continuity of the normal component of 𝜿u across interfaces together with the strongly enforced boundary conditions in U¯h,0p¯h to insert vF into the second term in parentheses, we infer


Setting, for the sake of conciseness (cf. Proposition 4.1),


and using the definition (3.2) of rTpT+1v¯T with w=uˇT, we have


Subtracting (4.12) from (4.10), and observing that the first terms inside the summations cancel out owing to (4.1), we have


Denote by 𝔗1(T) and 𝔗2(T) the two summands in parentheses. Using the Cauchy–Schwarz inequality followed by the approximation properties (4.2) of uˇT (with l=s=pT+1) and (4.16) below, we have for the first term


For the second term, the Cauchy–Schwarz inequality followed by (4.4) (with q=pT) readily yields


Using (4.14)–(4.15) to estimate the right-hand side of (4.13), applying the Cauchy–Schwarz inequality, and passing to the supremum yields (3.10). To prove (3.11), it suffices to observe that, inserting u¯^h and using the triangle inequality, we obtain


and (3.11) follows using the estimates (4.2), (4.4), and (3.10) to bound the terms in the right-hand side. ∎

Remark 4.4 (Nonconformity Error).

The name nonconformity error for the quantity (4.9) is reminiscent of the Finite Element literature. As a matter of fact Eh corresponds to the contribution that appears in the estimate provided by the Second Strang Lemma when nonconforming methods are considered; see, e.g., [23, Lemma 2.25] and, in particular, the third term in the right-hand side of [23, (2.20)].

Proposition 4.5 (Estimate of Boundary Difference Seminorm).

It holds, for all v¯TU¯Tp¯T,



Let T𝒯h, v¯TU¯Tp¯T, and set, for the sake of brevity, vˇT:=rTpT+1v¯T. We have, for all FT,


where we have used the fact that pTpF (cf. (3.1)) to infer that vT|FpF(F) and thus insert πFpF in the first line, added and subtracted (vˇT-πTpTvˇT) in the second line, used the triangle inequality together with the definition (3.5) of the face-based residual δTFp¯T and the L2(F)-stability of πFpF in the third. To conclude, we observe that, if pT=0, then the discrete trace inequality (2.6) followed by Poincaré’s inequality yield


while, if pT1, then


where we have inserted ±πT0vˇT in the first line, used the discrete trace inequality (2.6) in the second line, the L2(T)-optimality of πTpT together with the approximation properties (2.5) (with l=pT and q=s=0) in the third line, and we have concluded observing that (pT+1)/pT2 and using the local Poincaré inequality (2.4) to infer


Plugging the above bounds for vˇT-πTpTvˇTF into (4.17), squaring the resulting inequality, multiplying it by κF/hT, summing over FT, and recalling the bound (2.2) on card(T), we deduce (4.16). ∎

4.4 L2-Error Estimate

Proof of Theorem 3.4.

We let zU solve (3.12) with g=u^h-u and set z¯^h:=I¯hp¯hz and, for all T𝒯h (cf. Proposition 4.1),


For the sake of brevity, we also let e¯h:=u¯^h-u¯hU¯h,0p¯h (recall the definition (3.9) of u¯^h), so that u^T-uT=eT for all T𝒯h. We start by observing that


where we have used the fact that -(𝜿z)=eh a.e. in Ω followed by element-by-element partial integration together with the continuity of the normal component of 𝜿Tz across interfaces and the strongly enforced boundary conditions in U¯h,0p¯h to insert eF into the last term.

In view of adding and subtracting ah(e¯h,z¯^h) to the right-hand side of (4.19), we next provide two useful reformulations of this quantity. First, we have


where we have added the quantity (f,z)-(𝜿u,z)=0 (cf. (1.2)) in the first line, we have passed to the second line using the definition (3.3) of aT (with u¯T=u¯^T and v¯T=z¯^T) together with the discrete problem (3.7) to infer


and we have concluded using the definitions (4.1) of ϖ𝜿,TpT+1 (together with (4.11) and (4.18)) and (2.1) of πTpT and πT1-ΔT. Second, using the definition (3.2) of rTpT+1 (with v¯T=e¯T), we obtain


where we have additionally used the fact that zˇT=ϖ𝜿,TpT+1z (cf. (4.18)) together with the definition (4.1) of ϖ𝜿,TpT+1 to replace zˇT by z in the first term in parentheses.

Thus, adding (4.20) and subtracting (4.21) from (4.19), we obtain after rearranging




Using the Cauchy–Schwarz inequality, the approximation properties (4.2) of zˇT (with l=pT+1 and s=1) together with the consistency properties (4.4) of sT (with q=0) for the first factor, and the bound (4.16) for the second factor, we get


For the second term, the Cauchy–Schwarz inequality followed by the approximation properties (4.2) of uˇT (with q=pT) and zˇT (with q=1), and the consistency properties (4.4) of sT (with q=pT and q=0 for the first and second factor, respectively) yield


Finally, for the third term we have, when pT=0,


while, when pT1,


where we have used the optimality of πTpT in the L2(T)-norm to pass to the second line and the approximation properties (2.5) of ΠTpT to conclude. Using (4.23)–(4.26) to bound the right-hand side of (4.22), and recalling the energy error estimate (3.10) and elliptic regularity (3.13) concludes the proof. ∎

Funding source: Agence Nationale de la Recherche

Award Identifier / Grant number: project HHOMM (ANR-15-CE40-0005)

Funding statement: We acknowledge the support of Agence Nationale de la Recherche through project HHOMM (ANR-15-CE40-0005).

A Proof of Lemma 2.3

Let K^d be a L-Lipschitz set (that is, such that its boundary can be locally parametrized by means of L-Lipschitz functions) with diam(K^)=1, and fix r0>1 and a d-cube R(r0) containing K^. In [6, Lemma 3.1], the following is shown: Suppose that K^ is a cube, then given a function vHs+1(K^), there exists a polynomial ΠK^lvl(K^) that satisfies


for qs+1. Actually, by analyzing the proof of that result, we can reinforce that statement not only for cubes, but for sets K^ which are regular enough, in the sense that there has to exist an extension operator E:Hs+1(K^)Hs+1(R(2r0)) such that

(A.2)E(v)s+1,R(2r0)Cvs+1,R(K^),E(v)=0 on R(2r0)R(32r0).

So the problem is to show the existence of such an extension (in any dimension d1). The minimal assumptions such that this is doable are studied in [32, Theorems 5 and 5, p. 181]. Namely, by means of a careful inspection of [32, Theorems 5 and 5] and their proofs, and in particular formulas (25), (30) and the end of the proof of Theorem 5 at p. 192, we get that the constant C in (A.2) depends on the Lipschitz constant L and on the (minimal) number of L-Lipschitz coverings of K^, that is, the number of open sets which cover K^ and in each of whom K^ can be parametrized by means of an L-Lipschitz function.[1] Thus, to get the hp-estimate (2.5), we have to show two facts:

  1. first, that (2.5) follows from (A.1) whenever we replace K^ with an element T of the mesh;

  2. second, that (A.1) actually holds with a uniform bounding constant for the elements of the mesh.

(i) Proof of (2.5) for regular elements. Assume for the moment that any T𝒯h satisfies the regularity assumptions needed to have the extension (A.2), and thus (A.1), with T^:=T/hT in place of K^. We additionally assume, without loss of generality, that the barycenter of T (and thus of T^) is . Then, by homogeneity, we get that, for every fHr(T), letting f^(𝒙):=f(𝒙/λ),


Thus, setting r=s+1, λ=hT, and f=v-𝔮 with 𝔮 a generic polynomial of degree l, we get by (A.1) (applied to v-𝔮 and ΠTl(v-𝔮) in place of v and ΠTlv, respectively),


Now, [22, Theorem 3.2] states that


for any il+1. Thus by the definition of the Hq-norm we conclude that


Notice that we exploited again (A.3) to scale back on T, and the fact that hT<1.

(ii) Proof of regularity under Assumption 2.1. To conclude the proof, we are left to show that Assumption 2.1 entails a uniform bound only in terms of ϱ for the Lipschitz constant of every element T𝒯h. To this aim, consider 𝒙T. Then, 𝒙S for some (convex) element of the submesh S𝔗h contained in T. Since the sets S form a submesh of 𝒯h, a uniform bound on the Lipschitz regularity of S would imply a bound on the Lipschitz regularity of T. Thus, we focus on the regularity of S. Since S is convex, we can cover S by means of 2(d+1) open sets Ui, such that SUi admits a local convex (and thus Lipschitz) parametrization ϕi, i.e., there exists an orthogonal coordinate system such that SUi is the graph of a Lipschitz function ϕi:Iid-1. This bound on the number of open sets Ui is crucial to get [32, Theorem 5] to work (clearly, thanks to (2.2), the bound on the number of Lipschitz coverings of T is bounded by a constant 2dN=c(d,ϱ)). We claim that each ϕi is 1/ϱ-Lipschitz.

Suppose that 𝒙Ui=:U and set ϕ:=ϕi. Up to a rotation and a rescaling, we can suppose that 𝒙= and ϕ(𝒙)=ϕ()=0. Let now rs be the inradius of S and hS be its diameter. By Assumption 2.1, we know that hsrS1ϱ. Let Brs be a ball contained in S of radius rS. Up to a further rotation of center 𝒙= of the coordinate system, we can suppose that BrS is centered on the xd-axis. In place of ϕ:I, it is useful to consider its Lipschitz extension ϕ~:d-1 defined by, denoting by || the usual Euclidian norm,


We know that ϕ~ is Lipschitz on d-1 and that Lip(ϕ~)=Lip(ϕ) (see for instance [2, Proposition 2.12]). Moreover, it is clear that ϕ~ is convex on d-1. The fact that BrSS and it is centered on the xd-axis (without loss of generality, we can suppose that its center is 𝝃=(,ξd) with ξd>0) translates into the fact that BrS is contained in the epigraph of ϕ~ and its center has distance from d at most hS. Let now 𝒑ϕ~(0), where ϕ~ is the subdifferential of ϕ~. Then, for every 𝒚d-1 we have


By choosing 𝒚=λ𝒑, with λ0, we get the inequality


Since the epigraph of ϕ~ contains BrS, which is centered at a height less than hS on the xd-axis, and by the convexity of ϕ~, we have that the truncated cone


is contained in the epigraph of ϕ~ (see Figure 4). Then from (A.4) we get


and so, by Assumption 2.1,


Let now 𝒚d-1. Then


Since 𝒙 is arbitrary, this shows that Lip(ϕ)=Lip(ϕ~)=ϱ-1, and so that S is ϱ-1-Lipschitz.

Figure 4 Illustration for point (ii) in the proof of Lemma 2.3.

Figure 4

Illustration for point (ii) in the proof of Lemma 2.3.


[1] J. Aghili, S. Boyaval and D. A. Di Pietro, Hybridization of mixed high-order methods on general meshes and application to the Stokes equations, Comput. Methods Appl. Math. 15 (2015), no. 2, 111–134. 10.1515/cmam-2015-0004Search in Google Scholar

[2] L. Ambrosio, N. Fusco and D. Pallara, Functions of Bounded Variation and Free Discontinuity Problems, Oxford Math. Monogr., Oxford University Press, Oxford, 2000. Search in Google Scholar

[3] P. F. Antonietti, S. Giani and P. Houston, hp-version composite discontinuous Galerkin methods for elliptic problems on complicated domains, SIAM J. Sci. Comput. 35 (2013), no. 3, A1417–A1439. 10.1137/120877246Search in Google Scholar

[4] B. Ayuso de Dios, K. Lipnikov and G. Manzini, The nonconforming virtual element method, ESAIM Math. Model Numer. Anal. (M2AN) 50 (2016), no. 3, 879–904. 10.1051/m2an/2015090Search in Google Scholar

[5] I. Babuška and M. Suri, The h-p version of the finite element method with quasi-uniform meshes, RAIRO Modél. Math. Anal. Numér. 21 (1987), no. 2, 199–238. 10.1051/m2an/1987210201991Search in Google Scholar

[6] I. Babuška and M. Suri, The optimal convergence rate of the p-version of the finite element method, SIAM J. Numer. Anal. 24 (1987), no. 4, 750–776. 10.1137/0724049Search in Google Scholar

[7] I. Babuška, B. A. Szabo and I. N. Katz, The p-version of the finite element method, SIAM J. Numer. Anal. 18 (1981), no. 3, 515–545. 10.1137/0718033Search in Google Scholar

[8] L. Beirão da Veiga, A. Chernov, L. Mascotto and A. Russo, Basic principes of hp virtual elements on quasiuniform meshes, Math. Models Methods Appl. Sci. 26 (2016), no. 8, 1567–1598. 10.1142/S021820251650038XSearch in Google Scholar

[9] A. Cangiani, E. H. Georgoulis and P. Houston, hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes, Math. Models Methods Appl. Sci. 24 (2014), no. 10, 2009–2041. 10.1007/978-3-319-67673-9Search in Google Scholar

[10] P. Castillo, B. Cockburn, D. Scötzau and C. Schwab, Optimal a priori error estimates for the hp-version of the local discontinuous Galerkin method for convection-diffusion problems, Math. Comp. 71 (2001), no. 238, 455–478. 10.1090/S0025-5718-01-01317-5Search in Google Scholar

[11] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, Classics Appl. Math. 40, Society for Industrial and Applied Mathematics, Philadelphia, 2002. 10.1137/1.9780898719208Search in Google Scholar

[12] B. Cockburn, D. A. Di Pietro and A. Ern, Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods, ESAIM Math. Model. Numer. Anal. (M2AN) 50 (2016), no. 3, 635–650. 10.1051/m2an/2015051Search in Google Scholar

[13] B. Cockburn, J. Gopalakrishnan and R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), no. 2, 1319–1365. 10.1137/070706616Search in Google Scholar

[14] D. A. Di Pietro and J. Droniou, A hybrid high-order method for Leray–Lions elliptic equations on general meshes, Math. Comp. 86 (2016), no. 307, 2159–2191. 10.1090/mcom/3180Search in Google Scholar

[15] D. A. Di Pietro and J. Droniou, Ws,p-approximation properties of elliptic projectors on polynomial spaces, with application to the error analysis of a hybrid high-order discretisation of Leray–Lions problems, Math. Models Methods Appl. Sci. 27 (2017), no. 5, 879–908. 10.1142/S0218202517500191Search in Google Scholar

[16] D. A. Di Pietro and A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Math. Appl. (Berlin) 69, Springer, Berlin, 2012. 10.1007/978-3-642-22980-0Search in Google Scholar

[17] D. A. Di Pietro and A. Ern, A hybrid high-order locking-free method for linear elasticity on general meshes, Comput. Methods Appl. Mech. Engrg. 283 (2015), 1–21. 10.1016/j.cma.2014.09.009Search in Google Scholar

[18] D. A. Di Pietro, A. Ern and S. Lemaire, An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators, Comput. Methods Appl. Math. 14 (2014), no. 4, 461–472. 10.1515/cmam-2014-0018Search in Google Scholar

[19] D. A. Di Pietro, A. Ern and S. Lemaire, A review of hybrid high-order methods: Formulations, computational aspects, comparison with other methods, Building Bridges: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations (Durham 2014), Lecture Notes in Comput. Sci. 114, Springer, Cham (2016), 205–236. 10.1007/978-3-319-41640-3_7Search in Google Scholar

[20] D. A. Di Pietro and S. Lemaire, An extension of the Crouzeix–Raviart space to general meshes with application to quasi-incompressible linear elasticity and Stokes flow, Math. Comp. 84 (2015), no. 291, 1–31. 10.1090/S0025-5718-2014-02861-5Search in Google Scholar

[21] D. A. Di Pietro and R. Tittarelli, An introduction to hybrid high-order methods, Lectures from the Fall 2016 Thematic Quarter at Institut Henri Poincaré, Springer, to appear. 10.1007/978-3-319-94676-4_4Search in Google Scholar

[22] T. Dupont and R. Scott, Polynomial approximation of functions in Sobolev spaces, Math. Comp. 34 (1980), no. 150, 441–463. 10.1090/S0025-5718-1980-0559195-7Search in Google Scholar

[23] A. Ern and J.-L. Guermond, Theory and Practice Of Finite Elements, Appl. Math. Sci. 159, Springer, New York, 2004. 10.1007/978-1-4757-4355-5Search in Google Scholar

[24] E. H. Georgoulis and E. Süli, Optimal error estimates for the hp-version interior penalty discontinuous Galerkin finite element method, IMA J. Numer. Anal. 25 (2005), 205–220. 10.1093/imanum/drh014Search in Google Scholar

[25] S. Giani and P. Houston, hp-adaptive composite discontinuous Galerkin methods for elliptic problems on complicated domains, Numer. Methods Partial Differential Equations 30 (2014), no. 4, 1342–1367. 10.1002/num.21872Search in Google Scholar

[26] R. Herbin and F. Hubert, Benchmark on discretization schemes for anisotropic diffusion problems on general grids, Finite Volumes for Complex Applications V, John Wiley & Sons, New York (2008), 659–692. Search in Google Scholar

[27] C. Le Potier, A finite volume method for the approximation of highly anisotropic diffusion operators on unstructured meshes, Finite Volumes for Complex Applications IV, ISTE, London (2005), 401–412. Search in Google Scholar

[28] I. Perugia and D. Schötzau, A hp-analysis of the local discontinuous Galerkin method for diffusion problems, J. Sci. Comput. 17 (2002), no. 1–4, 561–571. 10.1023/A:1015118613130Search in Google Scholar

[29] B. Rivière, M. F. Wheeler and V. Girault, Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems. Part I, Comput. Geosci. 3 (1999), 337–360. 10.1023/A:1011591328604Search in Google Scholar

[30] C. Schwab, p- and hp-FEM. Theory and Application to Solid and Fluid Mechanics, Oxford University Press, Oxford, 1998. Search in Google Scholar

[31] B. Stamm and T. P. Wihler, hp-optimal discontinuous Galerkin methods for linear elliptic problems, Math. Comp. 79 (2010), no. 272, 2117–2133. 10.1090/S0025-5718-10-02335-5Search in Google Scholar

[32] E. M. Stein, Singular Integrals and Differentiability Properties of Functions, Princeton Math. Ser. 30, Princeton University Press, Princeton, 1970. 10.1515/9781400883882Search in Google Scholar

[33] C. Talischi, G. H. Paulino, A. Pereira and I. F. M. Menezes, PolyMesher: A general-purpose mesh generator for polygonal elements written in Matlab, Struct. Multidiscip. Optim. 45 (2012), no. 3, 309–328. 10.1007/s00158-011-0706-zSearch in Google Scholar

[34] M. Vohralík, A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations, SIAM J. Numer. Anal. 45 (2007), no. 4, 1570–1599. 10.1137/060653184Search in Google Scholar

Received: 2017-4-4
Revised: 2017-4-6
Accepted: 2017-4-7
Published Online: 2017-6-17
Published in Print: 2017-7-1

© 2017 Walter de Gruyter GmbH, Berlin/Boston

Scroll Up Arrow