An introduction to Hybrid High-Order methods

An introduction to Hybrid High-Order methods

Daniele A. Di Pietro and Roberta Tittarelli Daniele A. Di Pietro Roberta Tittarelli Institut Montpelliérain Alexander Grothendieck, CNRS, Univ. de Montpellier. 1, 2

This chapter provides an introduction to Hybrid High-Order (HHO) methods. These are new generation numerical methods for PDEs with several advantageous features: the support of arbitrary approximation orders on general polyhedral meshes, the reproduction at the discrete level of relevant continuous properties, and a reduced computational cost thanks to static condensation and compact stencil. After establishing the discrete setting, we introduce the basics of HHO methods using as a model problem the Poisson equation. We describe in detail the construction, and prove a priori convergence results for various norms of the error as well as a posteriori estimates for the energy norm. We then consider two applications: the discretization of the nonlinear -Laplace equation and of scalar diffusion-advection-reaction problems. The former application is used to introduce compactness analysis techniques to study the convergence to minimal regularity solution. The latter is used to introduce the discretization of first-order operators and the weak enforcement of boundary conditions. Numerical examples accompany the exposition.

1 Introduction

This chapter provides an introduction to Hybrid High-Order (HHO) methods. The material is closely inspired by a series of lectures given by the first author at Institut Henri Poincaré in September 2016 within the thematic quarter Numerical Methods for PDEs (see

HHO methods, introduced in Di-Pietro.Ern.ea:14 ; Di-Pietro.Ern:15 , are discretization methods for Partial Differential Equations (PDEs) with relevant features that set them apart from classical techniques such as finite elements or finite volumes. These include, in particular:

  1. The support of general polytopal meshes in arbitrary space dimension, paving the way to a seamless treatment of complex geometric features and unified 1d-2d-3d implementations;

  2. The possibility to select the approximation order which, possibly combined with adaptivity, leads to a reduction of the simulation cost for a given precision or better precision for a given cost;

  3. The compliance with the physics, including robustness with respect to the variations of physical coefficients and reproduction at the discrete level of key continuous properties such as local balances and flux continuity;

  4. A reduced computational cost thanks to their compact stencil along with the possibility to perform static condensation.

As of today, HHO methods have been successfully applied to the discretization of several linear and nonlinear problems of engineering interest including: variable diffusion Di-Pietro.Ern.ea:14 ; Di-Pietro.Ern.ea:16 ; Di-Pietro.Ern:16 , quasi incompressible linear elasticity Di-Pietro.Ern:15 ; Di-Pietro.Ern:15*1 , locally degenerate diffusion-advection-reaction Di-Pietro.Droniou.ea:15 , poroelasticity Boffi.Botti.ea:16 , creeping flows Aghili.Boyaval.ea:15 possibly driven by volumetric forces with large irrotational part Di-Pietro.Ern.ea:16*1 , electrostatics Di-Pietro.Specogna:16 , phase separation problems governed by the Cahn–Hilliard equation Chave.Di-Pietro.ea:16 , Leray–Lions type elliptic problems Di-Pietro.Droniou:16 ; Di-Pietro.Droniou:17 . More recent applications also include steady incompressible flows governed by the Navier–Stokes equations Di-Pietro.Krell:16 and nonlinear elasticity Botti.Di-Pietro.ea:16 . Generalizations of HHO methods and comparisons with other (new generation or classical) discretization methods for PDEs can be found in Cockburn.Di-Pietro.ea:16 ; Boffi.Di-Pietro:16 . Implementation tools based on advanced programming techniques have been recently discussed in Cicuttin.Di-Pietro.ea:17 .

Discretization methods that support polytopal meshes and, possibly, arbitrary approximation orders have experienced a vigorous development over the last decade. Novel approaches to the analysis and the design have been developed borrowing ideas from other branches of mathematics (such as topology and geometry), or expanding past their initial limits the original ideas underlying finite element or finite volume methods. A brief state-of-the-art is provided in what follows.

Several lowest-order methods for diffusive problems have been proposed to circumvent the strict conditions of mesh-data compliance required for the consistency of classical (two-points) finite volume schemes; see Droniou:14 for a comprehensive review. We mention here, in particular, the Mixed and Hybrid Finite Volume methods of Droniou.Eymard:06 ; Eymard.Gallouet.ea:10 . These methods possess local conservation properties on the primal mesh, and enable an explicit identification of equilibrated numerical fluxes. Their relation with the lowest-order version of HHO methods has been studied in (Di-Pietro.Ern.ea:14, , Section 2.5) for pure diffusion and in (Di-Pietro.Droniou.ea:15, , Section 5.4) for advection-diffusion-reaction. Other families of lowest-order methods have been obtained by reproducing at the discrete level salient features of the continuous problem. Mimetic Finite Difference methods are derived by emulating the Stokes theorems to formulate counterparts of differential operators and of -products; cf. Brezzi.Lipnikov.ea:05 and Droniou.Eymard.ea:10 for a study of their relation with Mixed and Hybrid Finite Volume methods. In the Discrete Geometric Approach of Codecasa.Specogna.ea:10 as well as in Compatible Discrete Operators Bonelle.Ern:14 , formal links with the continuous operators are expressed in terms of Tonti diagrams. To different extents, the aforementioned methods owe to the seminal ideas of Whitney on geometric integration Whitney:57 . A different approach to lowest-order schemes on general meshes consists in extending classical properties of nonconforming and penalized finite elements as in the Cell Centered Galerkin Di-Pietro:12 and generalized Crouzeix–Raviart Di-Pietro.Lemaire:15 methods. We also cite here Vohralik.Wohlmuth:13 concerning the use of classical mixed finite elements on polyhedral meshes (see, in particular, Section 7 therein). Further investigations have recently lead to unifying frameworks that encompass the above (and other) methods. We mention, in particular, the Gradient Schemes discretizations of Droniou.Eymard.ea:13 . Finally, the methods discussed here can often be regarded as lowest-order versions of more recent technologies.

Methods that support the possibility to increase the approximation order have received a considerable amount of attention over the last few years. High-order discretizations on general meshes that are possibly physics-compliant can be obtained by the discontinuous Galerkin approach; cf., e.g., Arnold.Brezzi.ea:02 ; Di-Pietro.Ern:12 and also Bassi.Botti.ea:12 . Discontinuous Galerkin methods, however, have some practical limitations. For problems in incompressible fluid mechanics, e.g., a key ingredient for inf-sup stability is a reduction map that can play the role of a Fortin interpolator. Unfortunately, such an interpolator is not available for discontinuous Galerkin methods on non-standard elements. Additionally, in particular for modal implementations on general meshes, the number of unknowns can become unbearably large. This has motivated the introduction of Hybridizable Discontinuous Galerkin methods Castillo.Cockburn.ea:00 ; Cockburn.Gopalakrishnan.ea:09 , which mainly focus on standard meshes (the extension to general meshes is possible in some cases); see also the very recent -decomposition techniques Cockburn.Fu:17 . High-order discretization methods that support general meshes also include Virtual Element methods; cf. Beirao-da-Veiga.Brezzi.ea:13*1 for an introduction. In short, Virtual Element methods are finite element methods where explicit expressions for the basis functions are not available at each point, and computable approximations thereof are used instead. This provides the extra flexibility required, e.g., to handle polyhedral elements. Links between HHO and the nonconforming Virtual Element method have been pointed out in (Cockburn.Di-Pietro.ea:16, , Section 2.4); see also Boffi.Di-Pietro:16 concerning HHO and mixed Virtual Element methods.

We next describe in detail the content of this chapter. We start in Section 2 by presenting the discrete setting: we introduce the notion of polytopal mesh (Section 2.1), formulate assumptions on the way meshes are refined that are suitable to carry out a -convergence analysis (Section 2.2), introduce the local polynomial spaces (Section 2.3) and projectors (Section 2.4) that lie at the heart of the HHO construction.

In Section 3 we present the basic principles of HHO methods using as a model problem the Poisson equation. While the material in this section is mainly adapted from Di-Pietro.Ern.ea:14 , some results are new and the arguments have been shortened or made more elegant. In Section 3.1 we introduce the local space of degrees of freedom (DOFs) and discuss the main ingredients upon which HHO methods rely, namely:

  1. Reconstructions of relevant quantities obtained by solving small, embarassingly parallel problems on each element;

  2. High-order stabilization terms obtained by penalizing cleverly designed residuals.

In Section 3.2 we show how to combine these ingredients to formulate local contributions, which are then assembled element-by-element as in standard finite elements. The construction is conceived so that only face-based DOFs are globally coupled, which paves the way to efficient practical implementations where element-based DOFs are statically condensed in a preliminary step. In Sections 3.3 and 3.4 we discuss, respectively, optimal a priori estimates for various norms and seminorms of the error, and residual-based a posteriori estimates for the energy-norm of the error. Finally, some numerical examples are provided in Section 3.5 to demonstrate the theoretical results.

In Section 4 we consider the HHO discretization of the -Laplace equation. The material is inspired by Di-Pietro.Droniou:16 ; Di-Pietro.Droniou:17 , where more general Leray–Lions operators are considered. When dealing with nonlinear problems, regularity for the exact solution is often difficult to prove and can entail stringent assumptions on the data. For this reason, the -convergence analysis can be carried out in two steps: in a first step, convergence to minimal regularity solutions is proved by a compactness argument; in a second step, convergence rates are estimated for smooth solutions (and smooth data). Convergence by compactness typically requires discrete counterparts of functional analysis results relevant for the study of the continuous problem. In our case, two sets of discrete functional analysis results are needed: discrete Sobolev embeddings (Section 4.1) and compactness for sequences of HHO functions uniformly bounded in a -like seminorm (Section 4.2). The interest of both results goes beyond the specific method and problem considered here. As an example, in Di-Pietro.Krell:16 they are used for the analysis of a HHO discretization of the steady incompressible Navier–Stokes equations. The HHO method for the -Laplacian stated in Section 4.3 is designed according to similar principles as for the Poisson problem. Convergence results are stated in Section 4.4, and numerical examples are provided in Section 4.5.

Following Di-Pietro.Droniou.ea:15 , in Section 5 we extend the HHO method to diffusion-advection-reaction problems. In this context, a crucial property from the numerical point of view is robustness in the advection-dominated regime. In Section 5.1 we modify the diffusive bilinear form introduced in Section 3.2 to incorporate weakly enforced boundary conditions. The weak enforcement of boundary conditions typically improves the behaviour of the method in the presence of boundary layers, since the discrete solution is not constrained to a fixed value on the boundary. In Section 5.2 we introduce the HHO discretization of first-order terms based on two novel ingredients: a local advective derivative reconstruction and an upwind penalty term. The former is used to formulate the consistency terms, while the role of the latter is to confer suitable stability properties to the advective-reactive bilinear form. The HHO discretization is finally obtained in Section 5.3 combining the diffusive and advective-reactive contributions, and its stability with respect to an energy-like norm including an advective derivative contribution is studied. In Section 5.4 we state an energy-norm error estimate which accounts for the dependence of the error contribution of each mesh element on a local Péclet number. A numerical illustration is provided in Section 5.5.

2 Discrete setting

Let , , denote a bounded connected open polyhedral domain with Lipschitz boundary and outward normal . We assume that does not have cracks, i.e., it lies on one side of its boundary. In what follows, we introduce the notion of polyhedral mesh of , formulate assumptions on the way meshes are refined that enable to prove useful geometric and functional results, and introduce functional spaces and projectors that will be used in the construction and analysis of HHO methods.

2.1 Polytopal mesh

The following definition enables the treatment of meshes as general as the ones depicted in Fig. 1.

(a) Matching triangular
(b) Nonconforming
(c) Polygonal
(d) Agglomerated
Figure 1: Examples of polytopal meshes in two and three space dimensions. The triangular and nonconforming meshes are taken from the FVCA5 benchmark Herbin.Hubert:08 , the polygonal mesh family from (Di-Pietro.Lemaire:15, , Section 4.2.3), and the agglomerated polyhedral mesh from Di-Pietro.Specogna:16 .
Definition 1 (Polytopal mesh)

A polytopal mesh of is a couple where:

(i) The set of mesh elements is a finite collection of nonempty disjoint open polytopes with boundary and diameter such that the meshsize satisfies and it holds that (ii) The set of mesh faces is a finite collection of disjoint subsets of such that, for any , is an open subset of a hyperplane of , the -dimensional Hausdorff measure of is strictly positive, and the -dimensional Hausdorff measure of its relative interior is zero. Moreover, (a) for each , either there exist two distinct mesh elements such that and is called an interface or there exists one mesh element such that and is called a boundary face; (b) the set of faces is a partition of the mesh skeleton, i.e.,

Interfaces are collected in the set and boundary faces in , so that . For any mesh element ,

denotes the set of faces contained in . Similarly, for any mesh face ,

is the set of mesh elements sharing . Finally, for all , is the unit normal vector to pointing out of .

Figure 2: Treatment of a nonconforming junction (red) as multiple coplanar faces. Gray elements are pentagons with two coplanar faces, white elements are squares.
Remark 1 (Nonconforming junctions)

Meshes including nonconforming junctions such as the one depicted in Fig. 2 are naturally supported provided that each face containing hanging nodes is treated as multiple coplanar faces.

2.2 Regular mesh sequences

When studying the convergence of HHO methods with respect to the meshsize , one needs to make assumptions on how the mesh is refined. The ones provided here are closely inspired by (Di-Pietro.Ern:12, , Chapter 1), and refer to the case of isotropic meshes with non-degenerate faces. Isotropic means here that we do not consider the case of elements that become more and more stretched when refining. Non-degenerate faces means, on the other hand, that the diameter of each mesh face is uniformly comparable to that of the element(s) it belongs to; see (2) below.

Definition 2 (Matching simplicial submesh)

Let be a polytopal mesh of . We say that is a matching simplicial submesh of if (i) is a matching simplicial mesh of ; (ii) for all simplices , there is only one mesh element such that ; (iii) for all , the set collecting the simplicial faces of , there is at most one face such that .

If itself is matching simplicial and collects the corresponding simplicial faces, we can simply take , so that . The notion of regularity for refined mesh sequences is made precise by the following

Definition 3 (Regular mesh sequence)

Denote by a countable set of meshsizes having as its unique accumulation point. A sequence of refined meshes is said to be regular if there exists a real number such that, for all , there exists a matching simplicial submesh of and (i) for all simplices of diameter and inradius , ; (ii) for all mesh elements and all simplices such that , .

Remark 2 (Role of the simplicial submesh)

The simplicial submesh introduced in Definition 3 is merely a theoretical tool, and needs not be constructed in practice.

Geometric bounds on regular mesh sequences can be proved as in (Di-Pietro.Ern:12, , Section 1.4.2) (the definition of mesh face is slightly different therein since planarity is not required, but the proofs are based on the matching simplicial submesh and one can check that they carry out unchanged). We recall here, in particular, that the number of faces of one mesh element is uniformly bounded: There is such that


Morevover, according to (Di-Pietro.Ern:12, , Lemma 1.42), for all , all , and all


Discrete functional analysis results for arbitrary-order methods on regular mesh sequences can be found in (Di-Pietro.Ern:12, , Chapter 1) and Di-Pietro.Droniou:16 ; Di-Pietro.Droniou:17 . We also refer the reader to Eymard.Gallouet.ea:00 for a first theorization of discrete functional analysis in the context of lowest-order finite volume methods, as well as to the subsequent extensions of Eymard.Gallouet.ea:10 ; Droniou.Eymard.ea:17 .

Throughout the rest of this work, it is tacitly understood that we work on regular mesh sequences.

2.3 Local and broken spaces

Throughout the rest of this chapter, for any , we denote by and the standard -product and norm, with the convention that the subscript is omitted whenever . The same notation is used for the vector-valued space .

Let now the set be a mesh element or face. For an integer , we denote by the space spanned by the restriction to of scalar-valued, -variate polynomials of total degree . We note the following trace inequality (see (Di-Pietro.Ern:12, , Lemma 1.46)): There is a real number only depending on , , and such that, for all , all , all , and all ,


At the global level, we define the broken polynomial space

Functions in belong to the broken Sobolev space

We denote by the usual broken gradient operator such that, for all ,

2.4 Projectors on local polynomial spaces

Projectors on local polynomial spaces play a key role in the design and analysis of HHO methods.

-orthogonal projector

Let denote a mesh element or face. The -orthogonal projector (in short, -projector) is defined as follows: For all , is the unique polynomial in that satisfies


Existence and uniqueness of follow from the Riesz representation theorem in for the standard -inner product. Moreover, we have the following characterization:

In what follows, we will also need the vector-valued -projector denoted by and obtained by applying component-wise. The following -boundedness result is a special case of (Di-Pietro.Droniou:16, , Corollary 3.7): For any , there exists a real number depending only on , , , and such that, for all , all , and all ,


At the global level, we denote by the -projector on the broken polynomial space such that, for all ,

Elliptic projector

For any mesh element , we also define the elliptic projector as follows: For all , is a polynomial in that satisfies

By the Riesz representation theorem in for the -inner product, this relation defines a unique element , and thus a polynomial up to an additive constant. This constant is fixed by writing

Observing that (6a) is trivially verified when , it follows from (6b) that . Finally, the following characterization holds:

Approximation properties

On regular mesh sequences, both and have optimal approximation properties in , as summarized by the following result (for a proof, see Theorem 1, Theorem 2, and Lemma 13 in Di-Pietro.Droniou:16 ): For any and , there exists a real number depending only on , , , , and such that, for all , all , and all ,

and, if ,
where .

3 Basic principles of Hybrid High-Order methods

To fix the main ideas and notation, we study in this section the HHO discretization of the Poisson problem: Find such that


where is a given volumetric source term. More general boundary conditions can replace (8b), but we restrict the discussion to the homogeneous Dirichlet case for the sake of simplicity. A detailed treatment of more general boundary conditions including also variable diffusion coefficients can be found in Di-Pietro.Ern.ea:16 .

The starting point to devise a HHO discretization is the following weak formulation of problem (8): Find such that


where the bilinear form is such that


In what follows, the quantities and will be referred to, respectively, as the potential and the flux.

3.1 Local construction

Throughout this section, we fix a polynomial degree and a mesh element . We introduce the local ingredients underlying the HHO construction: the DOFs, the potential reconstruction operator, and the discrete counterpart of the restriction to of the global bilinear form defined by (10).

Computing the local elliptic projection from -projections

Consider a function . We note the following integration by parts formula, valid for all :


Specializing (11) to , we obtain

where we have used (6) to insert into the left-hand side and (4) to insert and into the right-hand side after observing that and for all . Moreover, recalling (6b) and using the definition (4) of the -projector, we infer that

The relations (12) show that computing the elliptic projection does not require a full knowledge of the function . All that is required is

  1. , the -projection of on the polynomial space . Clearly, one could also choose instead, which has the advantage of not requiring a special treatment of the case ;

  2. for all , , the -projection of the trace of on on the polynomial space .

Local space of degrees of freedom

The remark at the end of the previous section motivates the introduction of the following space of DOFs (see Fig. 3):

Figure 3: DOFs in for .

Observe that naming space of DOFs involves a shortcut: the actual DOFs can be chosen in several equivalent ways (polynomial moments, point values, etc.), and the specific choice does not affect the following discussion. For a generic vector of DOFs in , we use the underlined notation . On , we define the -like seminorm such that, for all ,


where denotes the diameter of . The negative power of in the second term ensures that both contributions have the same scaling. The DOFs corresponding to a smooth function are obtained via the reduction map such that


Potential reconstruction operator

Inspired by formula (12), we introduce the potential reconstruction operator such that, for all ,


Notice that is a polynomial function on one degree higher than the element-based DOFs . By definition, for all it holds that


i.e., the composition of the potential reconstruction operator with the reduction map gives the elliptic projector on . An immediate consequence of (17) together with (7) is that has optimal approximation properties in .

Local contribution

We approximate the restriction to of the continuous bilinear form defined by (10) by the discrete bilinear form such that


where the first term in the right-hand side is the usual Galerkin contribution, while the second is a stabilization contribution for which we consider the following design conditions, originally proposed in Boffi.Di-Pietro:16 :

Assumption 1 (Local stabilization bilinear form )

The local stabilization bilinear form satisfies the following properties:

  1. Symmetry and positivity. is symmetric and positive semidefinite;

  2. Stability. There is a real number independent of and of , but possibly depending on , , and , such that

  3. Polynomial consistency. For all and all , it holds that


These requirements suggest that can be obtained penalizing in a least square sense residuals that vanish for reductions of polynomial functions in . Paradigmatic examples of such residuals are provided by the operators and, for all , such that, for all ,


To check that vanishes when with , we observe that

where we have used the definition of in the first equality, the relation (17) to replace by and the fact that to cancel from the second term in parentheses, and the fact that leaves polynomials of total degree up to unaltered as a projector to conclude. A similar argument shows that for all whenever .

Accounting for dimensional homogeneity with the Galerkin term, one possible expression for is thus


This choice, inspired by the Virtual Element literature Beirao-da-Veiga.Brezzi.ea:13 , differs from the original HHO stabilization of Di-Pietro.Ern.ea:14 , where the following expression is considered instead:


In this case, only quantities at faces are penalized. Both of the above expressions match the design conditions (S1)–(S3) and are essentially equivalent in terms of implementation. A detailed proof for as in (23) can be found in (Di-Pietro.Ern.ea:14, , Lemma 4). Yet another example of stabilization bilinear form used in the context of HHO methods is provided by (Aghili.Boyaval.ea:15, , Eq. (3.24)). This expression results from the hybridization of the Mixed High-Order method of Di-Pietro.Ern:16 .

Remark 3 (Original HDG stabilization)

The following stabilization bilinear form is used in the original Hybridizable Discontinuous Galerkin (HDG) method of Castillo.Cockburn.ea:00 ; Cockburn.Gopalakrishnan.ea:09 :

While this choice obviously satisfies the properties (S1)-(S2), it fails to satisfy (S3) (it is only consistent for polynomials of degree up to ). As a result, up to one order of convergence is lost with respect to the estimates of Theorems 3.1 and 3.2 below. For a discussion including fixes that restore optimal orders of convergence in HDG see Cockburn.Di-Pietro.ea:16 .

Consistency properties of the stabilization for smooth functions

In the following proposition we study the consistency properties of when its arguments are reductions of a smooth function. We give a detailed proof since this result is a new extension of the bound in (Di-Pietro.Ern.ea:14, , Theorem 8) (see, in particular, Eq. (45) therein) to more general stabilization bilinear forms.

Proposition 1 (Consistency of )

Let denote a stabilization bilinear form satisfying assumptions (S1)–(S3). Then, there is a real number independent of , but possibly depending on , , and , such that, for all and all , it holds that


We set, for the sake of brevity, and abridge as the inequality with multiplicative constant having the same dependencies as in (24). Using (S2) and (S3) we infer that


Recalling (14), we have that


Using the -boundedness of resulting from (5) with , and followed by the optimal approximation properties (7a) of (with , , , and ), it is inferred that


On the other hand, for all it holds that


where we have used the -boundedness of together with (2) and the discrete trace inequality (3) in the first line, a local Poincaré inequality resulting from (7a) with , , , and to pass to the second line, and the optimal approximation properties of expressed by (7a) with , , , and to conclude. Plugging (27) and (28) into (26), recalling that (see (1)), and using the resulting bound to estimate (25), (24) follows.∎

3.2 Discrete problem

We now show how to formulate the discrete problem from the local contributions introduced in the previous section.

Global spaces of degrees of freedom

We define the following global space of DOFs with single-valued interface unknowns:

Notice that single-valued means here that interface values match from one element to the adjacent one. For a generic element , we use the underlined notation and, for all , we denote by its restriction to . We also define the broken polynomial function such that

The DOFs corresponding to a smooth function are obtained via the reduction map such that

We define on the seminorm such that, for all ,


with local seminorm defined by (14). To account for the homogeneous Dirichlet boundary condition (8b) in a strong manner, we introduce the subspace

We recall the following discrete Poincaré inequality proved in (Di-Pietro.Droniou:16, , Proposition 5.4): There exists a real number independent of , but possibly depending on , , and , such that, for all ,

Proposition 2 (Norm )

The map defines a norm on .


The seminorm property being evident, it suffices to prove that, for all , . Let be such that . By (30), we have , hence for all . From the definition (14) of the norm , we also have that for all and all , hence . Since any mesh face belongs to the set for at least one mesh element , this concludes the proof.∎

Global bilinear form

We define the global bilinear forms and by element-by-element assembly setting, for all ,

Lemma 1 (Properties of )

The bilinear form enjoys the following properties:

  1. Stability. For all it holds with as in (19) that

  2. Consistency. There is a real number independent of , but possibly depending on , , and , such that, for all