Joubine Aghili, Daniele A. Di Pietro and Berardo Ruffini

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

Accessible
De Gruyter | Published online: June 17, 2017

# Abstract

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 L 2 -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 , d 1 , denote a bounded connected polytopal domain with Lipschitz boundary. We consider the variable diffusion model problem: Find u : Ω such that

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

where 𝜿 is a uniformly positive, symmetric, tensor-valued field on Ω, while f L 2 ( Ω ) 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 u U := H 0 1 ( Ω ) such that

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

where we have used the notation ( , ) for the usual inner products of both L 2 ( Ω ) and L 2 ( Ω ) 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. (i)

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

2. (ii)

elemental DOFs, consisting in d-variate polynomials of degree p T on each mesh element T, where p T 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 p F = p for all mesh faces F.

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

1. (i)

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

2. (ii)

a stabilization term penalizing face-based residuals and polynomially consistent up to degree ( p T + 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. (i)

It is valid for fairly general polytopal meshes.

2. (ii)

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

3. (iii)

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. (iv)

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

5. (v)

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

6. (vi)

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 L 2 -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 d 1 ; 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 ( p T + 1 ) - p T with respect to the polynomial degree p T 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 p T , cf. [24] (on more general meshes, the scaling for the symmetric interior penalty dG method is p T - ( p T - 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 L 2 -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 1 2 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 𝒯 h T ¯ and h = max T 𝒯 h h T ( h T 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. (i)

either there exist distinct T 1 , T 2 𝒯 h such that F = T 1 T 2 (and F is an interface)

2. (ii)

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

Interfaces are collected in h i , boundary faces in h b , and we let h := h i h b . 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 𝒯 h T = F h F . For all T 𝒯 h , the set T := { F h F T } collects the faces lying on the boundary of T and, for all F T , we denote by 𝒏 T F 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. (i)

ϱ h S r S for all h and all simplices S 𝔗 h of diameter h S and inradius r S .

2. (ii)

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

3. (iii)

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

### 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

𝜿 T := 𝜿 | T 0 ( T ) d × d .

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, H q ( X ) denotes the Hilbert space of functions which are in L 2 ( 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 L 2 ( 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 l 0 , 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 L 2 -projector π X l : L 1 ( X ) l ( X ) such that, for all w L 1 ( X ) ,

(2.1) ( π X l w - w , v ) X = 0 for all  v l ( 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,

(2.2) max h , T 𝒯 h card ( T ) N .

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

(2.3) v T 2 C ( v T v T + h T - 1 v T 2 ) ,

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

(2.4) v T C P h T v T ,

where C P = π - 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 h H , all T T h , all integers l 1 , all s 0 , and all v H s + 1 ( T ) , there exists a polynomial Π T l v P l ( T ) satisfying, for all 0 q s + 1 ,

(2.5) v - Π T l v q , T C h T min ( l , s ) - q + 1 l s + 1 - q v s + 1 , T .

### Proof.

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 h H , all T T h , all integers l 1 , and all v P l ( T ) , it holds

(2.6) v T C l h T 1 / 2 v T .

### Proof.

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 d 2 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 = ( p F ) F h h of skeletal polynomial degrees be given. For all T 𝒯 h , we denote by p ¯ T = ( p F ) F T the restriction of p ¯ h to T , and we introduce the following local space of DOFs:

(3.1) U ¯ T p ¯ T := p T ( T ) × ( × F T p F ( F ) ) , p T := min F T p F .

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

r T p T + 1 : U ¯ T p ¯ T p T + 1 ( T )

such that, for all v ¯ T U ¯ T p ¯ T and all w p T + 1 ( T ) ,

(3.2) ( 𝜿 T r T p T + 1 v ¯ T , w ) T = - ( v T , ( 𝜿 T w ) ) T + F T ( v F , 𝜿 T w 𝒏 T F ) F

and

( r T p T + 1 v ¯ T - v T , 1 ) T = 0 .

Note that computing r T p T + 1 v ¯ 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 ¯ T p ¯ T × U ¯ T p ¯ T the local bilinear form a T such that

(3.3) a T ( u ¯ T , v ¯ T ) := ( 𝜿 T r T p T + 1 u ¯ T , r T p T + 1 v ¯ T ) T + s T ( u ¯ T , v ¯ T ) ,
(3.4) s T ( u ¯ T , v ¯ T ) := F T κ F h T ( δ T F p ¯ T u ¯ T , δ T F p ¯ T v ¯ T ) F ,

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

δ T F p ¯ T : U ¯ T p ¯ T p F ( F )

is such that, for all v ¯ T U ¯ T p ¯ T ,

(3.5) δ T F p ¯ T v ¯ T := π F p F ( v F - r T p T + 1 v ¯ T + π T p T r T p T + 1 v ¯ T - v T ) .

The first contribution in a T is in charge of consistency, whereas the second ensures stability by a least-square penalty of the face-based residuals δ T F p ¯ T . This subtle form for δ T F p ¯ T ensures that the residual vanishes when its argument is the interpolate of a function in p T + 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 ¯ h p ¯ h := ( × T 𝒯 h p T ( T ) ) × ( × F h p F ( F ) ) , U ¯ h , 0 p ¯ h := { v ¯ h U ¯ h p ¯ h v F 0  for all  F h b } .

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

v ¯ h = ( ( v T ) T 𝒯 h , ( v F ) F h )

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

I ¯ h p ¯ h : H 1 ( Ω ) U ¯ h p ¯ h

such that, for all v H 1 ( Ω ) ,

(3.6) I ¯ h p ¯ h v := ( ( π T p T v ) T 𝒯 h , ( π F p F v ) F h ) ,

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

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

(3.7) a h ( u ¯ h , v ¯ h ) = l h ( v ¯ h ) for all  v ¯ h U ¯ h , 0 p ¯ h ,

where the global bilinear form a h on U ¯ h p ¯ h × U ¯ h p ¯ h and the linear form l h on U ¯ h p ¯ h are assembled element-wise setting

a h ( u ¯ h , v ¯ h ) := T 𝒯 h a T ( u ¯ T , v ¯ T ) , l h ( v ¯ h ) := T 𝒯 h ( f , v T ) T .

### 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

(3.8) N dof = F h i ( p F + d - 1 p F ) .

### 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 F h i , we denote by [ ] F the usual jump operator (the sign is irrelevant), which we extend to boundary faces F h b setting [ φ ] F := φ . Let

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

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

𝔘 T p ¯ T := { 𝔳 T H 1 ( T ) ( 𝜿 T 𝔳 T ) p T ( T )  and  𝜿 T 𝔳 T | F 𝒏 T F p F ( F )  for all  F T } .

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

I ¯ T p ¯ T : 𝔘 T p ¯ T U ¯ T p ¯ T

is an isomorphism. Thus, the triplet ( T , 𝔘 T p ¯ T , I ¯ T p ¯ 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 , 0 p ¯ h such that

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

where

𝔞 h ( 𝔲 h , 𝔳 h ) := a h ( I ¯ h p ¯ h 𝔲 h , I ¯ h p ¯ h 𝔳 h ) , l h ( 𝔳 h ) := l h ( I ¯ h p ¯ h 𝔳 h ) ,

and it can be proved that 𝔲 h is the unique element of 𝔘 h , 0 p ¯ h such that u ¯ h = I ¯ h p ¯ 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 ¯ T p ¯ T by the bilinear forms a T and s T , respectively, and by a , h the seminorm defined by the bilinear form a h on U ¯ h p ¯ h . We also introduce the penalty seminorm | | s , h such that, for all v ¯ h U ¯ h p ¯ h ,

| v ¯ h | s , h 2 := T 𝒯 h | v ¯ T | s , T 2 .

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

r h p ¯ h : U ¯ h p ¯ h L 2 ( Ω )

such that, for all v ¯ h U ¯ h p ¯ h ,

( r h p ¯ h v ¯ h ) = | T r T p T + 1 v ¯ T for all  T 𝒯 h .

Finally, for the sake of conciseness, throughout the rest of the paper we denote by a b the inequality a C b 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 u U and u ¯ h U ¯ h , 0 p ¯ h denote the unique solutions of problems (1.2) and (3.7), respectively, and set

(3.9) u ¯ ^ h := I ¯ h p ¯ h u .

Assuming that u | T H p T + 2 ( T ) for all T T h , it holds

(3.10) u ¯ h - u ¯ ^ h a , h ( T 𝒯 h κ ¯ T λ 𝜿 , T h T 2 ( p T + 1 ) ( p T + 1 ) 2 p T u p T + 2 , T 2 ) 1 / 2 .

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

(3.11) 𝜿 1 / 2 h ( u - r h p ¯ h u ¯ h ) 2 + | u ¯ h | s , h 2 T 𝒯 h κ ¯ T λ 𝜿 , T h T 2 ( p T + 1 ) ( p T + 1 ) 2 p T u p T + 2 , T 2 .

### Proof.

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 ( p T + 1 ) - p T (i.e., half a power more than discontinuous Galerkin methods based on polynomials of degree p T , 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

𝜿 1 / 2 h ( u - r h p ¯ h u ¯ ^ h ) 2 + | u ¯ ^ h | s , h 2 ,

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 1 2 ) is observed with respect to its local anisotropy ratio.

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

u h := | T u T and u ^ h := | T u ^ T = π T p T u for all  T 𝒯 h .

To this end, we need elliptic regularity in the following form: For all g L 2 ( Ω ) , the unique element z U such that

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

satisfies the a priori estimate

(3.13) z 2 κ ¯ - 1 g L 2 ( Ω ) , κ ¯ := min T 𝒯 h κ ¯ T .

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 vs. N dof

### (b)

u h - u ^ h vs. N dof

### (c)

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

### (d)

u h - u ^ h vs. p

Figure 2

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

### (a)

u ¯ h - u ¯ ^ h a , h vs. N dof

### (b)

u h - u ^ h vs. N dof

### (c)

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

### (d)

u h - u ^ h vs. p

### Theorem 3.4 ( L 2 -Error Estimate).

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

κ ¯ u h - u ^ h κ ¯ 1 / 2 λ 𝜿 h ( T 𝒯 h λ 𝜿 , T κ ¯ T h T 2 ( p T + 1 ) ( p T + 1 ) 2 p T u p T + 2 , T 2 ) 1 / 2 + ( T 𝒯 h h T 2 ( p T + 2 ) ( p T + Δ T ) 2 ( p T + 2 ) f p T + Δ T 2 ) 1 / 2

with λ 𝛋 := max T T h λ 𝛋 , T and κ ¯ := max T T h κ ¯ 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 ( π x 1 ) sin ( π x 2 ) and right-hand side f chosen accordingly. We consider two values for the diffusion coefficients:

𝜿 1 = 𝑰 2 , 𝜿 2 = ( ( x 2 - x ¯ 2 ) 2 + ϵ ( x 1 - x ¯ 1 ) 2 - ( 1 - ϵ ) ( x 1 - x ¯ 1 ) ( x 2 - x ¯ 2 ) - ( 1 - ϵ ) ( x 1 - x ¯ 1 ) ( x 2 - x ¯ 2 ) ( x 1 - x ¯ 1 ) 2 + ϵ ( x 2 - x ¯ 2 ) 2 ) ,

where 𝑰 2 denotes the identity matrix of dimension 2, 𝒙 ¯ := - ( 0.1 , 0.1 ) , and ϵ = 1 10 - 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 L 2 -errors as a function of the number of skeletal DOFs N dof (cf. (3.8)) when p F = p for all F h 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. N dof 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].

Triangular

Cartesian

Refined

Staggered

Hexagonal

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 r T p T + 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 l 1 , we define the oblique elliptic projector ϖ 𝜿 , T l : H 1 ( T ) l ( T ) such that, for all v H 1 ( T ) , ( ϖ 𝜿 , T l v - v , 1 ) T = 0 and it holds

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

### Proposition 4.1 (Characterization of r T p T + 1 ∘ I ¯ T p ¯ T ).

It holds, for all T T h ,

r T p T + 1 I ¯ T p ¯ T = ϖ 𝜿 , T p T + 1 .

### Proof.

For a generic v H 1 ( T ) , letting v ¯ T = I ¯ T p ¯ T v in (3.2), we infer, for all w p T + 1 ( T ) ,

( 𝜿 T ( r T p T + 1 I ¯ T p ¯ T ) v , w ) T = - ( π T p T v , ( 𝜿 T w ) ) T + F T ( π F p F v , 𝜿 T w 𝒏 T F ) F
= - ( v , ( 𝜿 T w ) ) T + F T ( v , 𝜿 T w 𝒏 T F ) F
= ( 𝜿 T v , w ) T ,

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

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

### Lemma 4.2 (Approximation Properties of ϖ κ , T l ).

For all integers l 1 , all mesh elements T T h , all 0 s l , and all v H s + 1 ( T ) , it holds

𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T + h T 1 / 2 l 𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T + κ ¯ T 1 / 2 h T v - ϖ 𝜿 , T l v T + κ ¯ T 1 / 2 h T 1 / 2 v - ϖ 𝜿 , T l v T
(4.2) κ ¯ T 1 / 2 h T min ( l , s ) l s v s + 1 , T .

### Proof.

By definition (4.1) of ϖ 𝜿 , T l , it holds

(4.3) 𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T 2 = min w l ( T ) 𝜿 T 1 / 2 ( v - w ) T 2 κ ¯ T ( v - Π T l v ) T 2 .

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

𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T κ ¯ T 1 / 2 h T min ( l , s ) l s v s + 1 , T .

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

𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T 𝜿 T 1 / 2 ( v - Π T l v ) T + 𝜿 T 1 / 2 ( Π T l v - ϖ 𝜿 , T l v ) T = : 𝔗 1 + 𝔗 2 .

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

𝔗 1 κ ¯ T 1 / 2 h T min ( l , s ) - 1 / 2 l s - 1 / 2 v s + 1 , T .

For the second term, we have

𝔗 2 l h T 1 / 2 𝜿 T 1 / 2 ( Π T l v - ϖ 𝜿 , T l v ) T
l h T 1 / 2 ( 𝜿 T 1 / 2 ( Π T l v - v ) T + 𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T )
κ ¯ T 1 / 2 l h T 1 / 2 ( Π T l v - v ) T
κ ¯ T 1 / 2 h T min ( l , s ) - 1 / 2 l s - 1 v s + 1 , T ,

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 - ϖ 𝜿 , T l v , 1 ) T = 0 , we apply the local Poincaré inequality (2.4) to infer

v - ϖ 𝜿 , T l v T h T κ ¯ T 1 / 2 𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T κ ¯ T 1 / 2 κ ¯ T 1 / 2 h T min ( l , s ) + 1 l s v s + 1 , T ,

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

v - ϖ 𝜿 , T l v T 2 κ ¯ T - 1 / 2 v - ϖ 𝜿 , T l v T 𝜿 T 1 / 2 ( v - ϖ 𝜿 , T l v ) T + h T - 1 v - ϖ 𝜿 , T l v T 2 ,

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 s T defined by (3.4) are summarized in the following lemma.

### Lemma 4.3 (Consistency of the Stabilization Term).

For all T T h , all 0 q p T , and all v H q + 2 ( T ) , it holds

(4.4) | I ¯ T p ¯ T v | s , T κ ¯ T 1 / 2 λ 𝜿 , T 1 / 2 h T min ( p T , q ) + 1 ( p T + 1 ) q v q + 2 , T .

### Proof.

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

v ˇ T := ( r T p T + 1 I ¯ T p ¯ T ) v = ϖ 𝜿 , T p T + 1 v

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

δ T F p ¯ T I ¯ T p ¯ T v = π F p F ( v ˇ T - v ) - π T p T ( v ˇ T - v ) .

Using the triangle inequality and the L 2 ( F ) -stability of π F p F , we infer,

(4.5) δ T F p ¯ T I ¯ T p ¯ T v F v ˇ T - v F + π T p T ( v ˇ T - v ) F = : 𝔗 1 + 𝔗 2 .

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

(4.6) 𝔗 1 λ 𝜿 , T 1 / 2 h T min ( p T , q ) + 3 / 2 ( p T + 1 ) q v q + 2 , T .

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

(4.7) 𝔗 2 ( p T + 1 ) h T 1 / 2 π T p T ( v ˇ T - v ) T ( p T + 1 ) h T 1 / 2 v ˇ T - v T λ 𝜿 , T 1 / 2 h T min ( p T , q ) + 3 / 2 ( p T + 1 ) q v q + 2 , T .

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 / h T , summing over F T , and using the bound (2.2) on card ( T ) . ∎

### Proof of Theorem 3.3.

We start by noting the following abstract error estimate:

(4.8) u ¯ h - u ¯ ^ h a , h sup v ¯ h U ¯ h , 0 p ¯ h , v ¯ h a , h = 1 E h ( v ¯ h ) ,

with nonconformity error

(4.9) E h ( v ¯ h ) := l h ( v ¯ h ) - a h ( u ¯ ^ h , v ¯ h ) .

To prove (4.8), it suffices to observe that

u ¯ h - u ¯ ^ h a , h 2 = a h ( u ¯ h - u ¯ ^ h , u ¯ h - u ¯ ^ h )
= a h ( u ¯ h , u ¯ h - u ¯ ^ h ) - a h ( u ¯ ^ h , u ¯ h - u ¯ ^ h )
= l h ( u ¯ h - u ¯ ^ h ) - a h ( u ¯ ^ h , u ¯ h - u ¯ ^ h ) ,

where we have used the definition of the a , h -norm in the first line, the linearity of a h 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 ¯ ^ h a , h , using linearity, and passing to the supremum.

We next bound the nonconformity error E h ( v ¯ h ) for a generic vector of DOFs v ¯ h U ¯ h , 0 p ¯ h . A preliminary step consists in finding a more appropriate rewriting for E h ( 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 , 0 p ¯ h to insert v F into the second term in parentheses, we infer

(4.10) l h ( v ¯ h ) = T 𝒯 h ( ( 𝜿 T u , v T ) T + F T ( 𝜿 T u 𝒏 T F , v F - v T ) F ) .

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

(4.11) u ˇ T := r T p T + 1 u ¯ ^ T = ϖ 𝜿 , T p T + 1 u ,

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

(4.12) a h ( u ¯ ^ h , v ¯ h ) = T 𝒯 h ( ( 𝜿 T u ˇ T , v T ) T + F T ( 𝜿 T u ˇ T 𝒏 T F , v F - v T ) F + s T ( u ¯ ^ T , v ¯ T ) ) .

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

(4.13) E h ( v ¯ h ) = T 𝒯 h ( F T ( 𝜿 T ( u ˇ T - u ) 𝒏 T F , v F - v T ) F + s T ( u ¯ ^ T , v ¯ T ) ) .

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 = p T + 1 ) and (4.16) below, we have for the first term

| 𝔗 1 ( T ) | h T 1 / 2 𝜿 T 1 / 2 ( u ˇ T - u ) T × ( F T κ F h T v F - v T F 2 ) 1 / 2
(4.14) κ ¯ T 1 / 2 λ 𝜿 , T 1 / 2 h T p T + 1 ( p T + 1 ) p T u p T + 2 , T v ¯ T a , T .

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

(4.15) | 𝔗 2 ( T ) | κ ¯ T 1 / 2 λ 𝜿 , T 1 / 2 h T p T + 1 ( p T + 1 ) p T u p T + 2 , T | v ¯ T | s , T κ ¯ T 1 / 2 λ 𝜿 , T 1 / 2 h T p T + 1 ( p T + 1 ) p T u p T + 2 , T v ¯ T a , T .

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

𝜿 1 / 2 h ( u - r h p ¯ h u ¯ h ) 2 + | u ¯ h | s , h 2 𝜿 1 / 2 h ( u - r h p ¯ h u ¯ ^ h ) 2 + | u ¯ ^ h | s , h 2 + u ¯ h - u ¯ ^ h a , h 2 .

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 E h 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 ¯ T U ¯ T p ¯ T ,

(4.16) F T κ F h T v F - v T F 2 λ 𝜿 , T v ¯ T a , T 2 .

### Proof.

Let T 𝒯 h , v ¯ T U ¯ T p ¯ T , and set, for the sake of brevity, v ˇ T := r T p T + 1 v ¯ T . We have, for all F T ,

v F - v T F = π F p F ( v F - v T ) F
= π F p F ( v F - v ˇ T + π T p T v ˇ T - v T + v ˇ T - π T p T v ˇ T ) F
(4.17) δ T F p ¯ T v ¯ T F + v ˇ T - π T p T v ˇ T F ,

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

v ˇ T - π T 0 v ˇ T F h T 1 / 2 κ ¯ T - 1 / 2 𝜿 v ˇ T T

while, if p T 1 , then

v ˇ T - π T p T v ˇ T F = v ˇ T - π T 0 v ˇ T - π T p T ( v ˇ T - π T 0 v ˇ T ) F
p T + 1 h T 1 / 2 v ˇ T - π T 0 v ˇ T - π T p T ( v ˇ T - π T 0 v ˇ T ) T
p T + 1 h T 1 / 2 h T p T v ˇ T - π T 0 v ˇ T 1 , T
h T 1 / 2 κ ¯ T - 1 / 2 𝜿 T 1 / 2 v ˇ T T ,

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

v ˇ T - π T 0 v ˇ T 1 , T h T 1 / 2 v ˇ T T .

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

### Proof of Theorem 3.4.

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

(4.18) z ˇ T := r T p T + 1 z ¯ ^ T = ϖ 𝜿 , T p T + 1 z .

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

(4.19) e h 2 = - ( ( 𝜿 z ) , e h ) = T 𝒯 h ( ( 𝜿 T z , e T ) T + F T ( 𝜿 T z 𝒏 T F , e F - e T ) F ) ,

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

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

a h ( e ¯ h , z ¯ ^ h ) = a h ( u ¯ ^ h , z ¯ ^ h ) - a h ( u ¯ h , z ¯ ^ h ) + ( f , z ) - ( 𝜿 u , z )
= T 𝒯 h ( ( 𝜿 T u ˇ T , z ˇ T ) T - ( 𝜿 T u , z ) T + s T ( u ¯ ^ T , z ¯ ^ T ) + ( f , z - π T p T z ) T )
(4.20) = T 𝒯 h ( ( 𝜿 T ( u ˇ T - u ) , ( z ˇ T - z ) ) T + s T ( u ¯ ^ T , z ¯ ^ T ) + ( f - π T p T f , z - π T 1 - Δ T z ) T ) ,

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 a T (with u ¯ T = u ¯ ^ T and v ¯ T = z ¯ ^ T ) together with the discrete problem (3.7) to infer

a h ( u ¯ h , z ¯ ^ h ) = ( f , z ^ h ) = T 𝒯 h ( f , π T p T z ) T ,

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

(4.21) a h ( e ¯ h , z ¯ ^ h ) = T 𝒯 h ( ( 𝜿 T z , e T ) T + F T ( 𝜿 T z ˇ T , e F - e T ) F + s T ( z ¯ ^ T , e ¯ T ) ) ,

where we have additionally used the fact that z ˇ T = ϖ 𝜿 , T p T + 1 z (cf. (4.18)) together with the definition (4.1) of ϖ 𝜿 , T p T + 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

(4.22) e h 2 = T 𝒯 h ( 𝔗 1 ( T ) + 𝔗 2 ( T ) + 𝔗 3 ( T ) ) ,

with

𝔗 1 ( T ) := F T ( 𝜿 T ( z - z ˇ T ) 𝒏 T F , e F - e T ) F + s T ( z ¯ ^ T , e ¯ T ) ,
𝔗 2 ( T ) := ( 𝜿 T ( u ˇ T - u ) , ( z ˇ T - z ) ) T + s T ( u ¯ ^ T , z ¯ ^ T ) ,
𝔗 3 ( T ) := ( f - π T p T f , z - π T 1 - Δ T z ) T .

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

| 𝔗 1 ( T ) | ( h T 𝜿 T 1 / 2 ( z - z ˇ T ) T 2 + | z ¯ ^ T | s , T 2 ) 1 / 2 × ( F T κ F h T e F - e T F 2 + | e ¯ T | s , T 2 ) 1 / 2
(4.23) κ ¯ T 1 / 2 λ 𝜿 , T h T e ¯ T a , T z 2 , T .

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

| 𝔗 2 ( T ) | ( 𝜿 1 / 2 ( u ˇ T - u ) T 2 + | u ¯ ^ T | s , T 2 ) 1 / 2 × ( 𝜿 1 / 2 ( z ˇ T - z ) T 2 + | z ¯ ^ T | s , T 2 ) 1 / 2
(4.24) κ ¯ T λ 𝜿 , T h T p T + 2 ( p T + 1 ) p T u p T + 2 , T z 2 , T .

Finally, for the third term we have, when p T = 0 ,

(4.25) | 𝔗 3 ( T ) | f - π T 0 f T z - π T 0 z T h T 2 f 1 , T z 1 , T h T 2 f 1 , T z 2 , T ,

while, when p T 1 ,

| 𝔗 3 ( T ) | f - π T p T f T z - π T 1 z T
f - Π T p T f T z - Π T 1 z T
(4.26) h T p T + 2 p T p T + 2 f p T z 2 , T h T p T + 2 p T p T + 2 f p T , T z 2 , T ,

where we have used the optimality of π T p T in the L 2 ( T ) -norm to pass to the second line and the approximation properties (2.5) of Π T p T 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)

## 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 r 0 > 1 and a d-cube R ( r 0 ) containing K ^ . In [6, Lemma 3.1], the following is shown: Suppose that K ^ is a cube, then given a function v H s + 1 ( K ^ ) , there exists a polynomial Π K ^ l v l ( K ^ ) that satisfies

(A.1) v - Π K ^ l v q , K ^ 1 l s + 1 - q v s + 1 , K ^

for q s + 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 : H s + 1 ( K ^ ) H s + 1 ( R ( 2 r 0 ) ) such that

(A.2) E ( v ) s + 1 , R ( 2 r 0 ) C v s + 1 , R ( K ^ ) , E ( v ) = 0  on  R ( 2 r 0 ) R ( 3 2 r 0 ) .

So the problem is to show the existence of such an extension (in any dimension d 1 ). 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. (i)

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

2. (ii)

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 / h T 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 f H r ( T ) , letting f ^ ( 𝒙 ) := f ( 𝒙 / λ ) ,

(A.3) f r , T λ d / 2 - r f ^ r , T λ .

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

v - Π T l v q , T = ( v - 𝔮 ) - Π T l ( v - 𝔮 ) q , T h T d / 2 - q l s + 1 - q v ^ - 𝔮 ^ s + 1 , T ^ .

Now, [22, Theorem 3.2] states that

| f - Π T ^ l f | i , T ^ | f | l , T ^

for any i l + 1 . Thus by the definition of the H q -norm we conclude that

v - Π T l v q , T h T d / 2 - q l s + 1 - q ( i = min ( l , s ) + 1 s + 1 | v ^ | i , T ^ 2 ) 1 / 2 h T min ( l , s ) - q + 1 l s + 1 - q v s + 1 , T .

Notice that we exploited again (A.3) to scale back on T, and the fact that h T < 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 U i , such that S U i admits a local convex (and thus Lipschitz) parametrization ϕ i , i.e., there exists an orthogonal coordinate system such that S U i is the graph of a Lipschitz function ϕ i : I i d - 1 . This bound on the number of open sets U i 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 2 d N = c ( d , ϱ ) ). We claim that each ϕ i is 1 / ϱ -Lipschitz.

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

ϕ ~ ( 𝒙 ) := inf { ϕ ( 𝒚 ) + Lip ( ϕ ) | 𝒚 - 𝒙 | 𝒚 d - 1 } ,

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 B r S S and it is centered on the x d -axis (without loss of generality, we can suppose that its center is 𝝃 = ( , ξ d ) with ξ d > 0 ) translates into the fact that B r S is contained in the epigraph of ϕ ~ and its center has distance from d at most h S . Let now 𝒑 ϕ ~ ( 0 ) , where ϕ ~ is the subdifferential of ϕ ~ . Then, for every 𝒚 d - 1 we have

ϕ ~ ( 𝒚 ) 𝒑 𝒚 .

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

(A.4) ϕ ~ ( λ 𝒑 ) λ | 𝒑 | | 𝒑 | .

Since the epigraph of ϕ ~ contains B r S , which is centered at a height less than h S on the x d -axis, and by the convexity of ϕ ~ , we have that the truncated cone

C = { ( 𝒙 , x n ) d - 1 × h S x n h S r S | 𝒙 | }

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

| 𝒑 | 2 ϕ ~ ( 𝒑 ) h S r S | 𝒑 |

and so, by Assumption 2.1,

| 𝒑 | h S r S 1 ϱ .

Let now 𝒚 d - 1 . Then

| ϕ ~ ( 𝒚 ) - ϕ ~ ( 𝒙 ) | | 𝒑 ( 𝒚 - 𝒙 ) | | 𝒑 | | 𝒚 - 𝒙 | ϱ - 1 | 𝒚 - 𝒙 | .

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.

### References

[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. Search 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. Search 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. Search 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. Search 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. Search 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. Search 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. Search 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. Search 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. Search 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. Search 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. Search 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. Search 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. Search in Google Scholar

[15] D. A. Di Pietro and J. Droniou, W s , 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. Search in Google Scholar

[16] D. A. Di Pietro and A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Math. Appl. (Berlin) 69, Springer, Berlin, 2012. Search 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. Search 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. Search 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. Search 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. Search 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. Search in Google Scholar

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

[23] A. Ern and J.-L. Guermond, Theory and Practice Of Finite Elements, Appl. Math. Sci. 159, Springer, New York, 2004. Search 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. Search 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. Search 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. Search 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. Search 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. Search in Google Scholar

[32] E. M. Stein, Singular Integrals and Differentiability Properties of Functions, Princeton Math. Ser. 30, Princeton University Press, Princeton, 1970. Search 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. Search 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. Search in Google Scholar