An introduction to Hybrid HighOrder methods
Abstract
This chapter provides an introduction to Hybrid HighOrder (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 diffusionadvectionreaction 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 firstorder operators and the weak enforcement of boundary conditions. Numerical examples accompany the exposition.
1 Introduction
This chapter provides an introduction to Hybrid HighOrder (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 http://tinyurl.com/IHPquarternmpdes).
HHO methods, introduced in DiPietro.Ern.ea:14 ; DiPietro.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:

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

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;

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;

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 DiPietro.Ern.ea:14 ; DiPietro.Ern.ea:16 ; DiPietro.Ern:16 , quasi incompressible linear elasticity DiPietro.Ern:15 ; DiPietro.Ern:15*1 , locally degenerate diffusionadvectionreaction DiPietro.Droniou.ea:15 , poroelasticity Boffi.Botti.ea:16 , creeping flows Aghili.Boyaval.ea:15 possibly driven by volumetric forces with large irrotational part DiPietro.Ern.ea:16*1 , electrostatics DiPietro.Specogna:16 , phase separation problems governed by the Cahn–Hilliard equation Chave.DiPietro.ea:16 , Leray–Lions type elliptic problems DiPietro.Droniou:16 ; DiPietro.Droniou:17 . More recent applications also include steady incompressible flows governed by the Navier–Stokes equations DiPietro.Krell:16 and nonlinear elasticity Botti.DiPietro.ea:16 . Generalizations of HHO methods and comparisons with other (new generation or classical) discretization methods for PDEs can be found in Cockburn.DiPietro.ea:16 ; Boffi.DiPietro:16 . Implementation tools based on advanced programming techniques have been recently discussed in Cicuttin.DiPietro.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 stateoftheart is provided in what follows.
Several lowestorder methods for diffusive problems have been proposed to circumvent the strict conditions of meshdata compliance required for the consistency of classical (twopoints) 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 lowestorder version of HHO methods has been studied in (DiPietro.Ern.ea:14, , Section 2.5) for pure diffusion and in (DiPietro.Droniou.ea:15, , Section 5.4) for advectiondiffusionreaction. Other families of lowestorder 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 lowestorder schemes on general meshes consists in extending classical properties of nonconforming and penalized finite elements as in the Cell Centered Galerkin DiPietro:12 and generalized Crouzeix–Raviart DiPietro.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 lowestorder 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. Highorder discretizations on general meshes that are possibly physicscompliant can be obtained by the discontinuous Galerkin approach; cf., e.g., Arnold.Brezzi.ea:02 ; DiPietro.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 infsup 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 nonstandard 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 . Highorder discretization methods that support general meshes also include Virtual Element methods; cf. BeiraodaVeiga.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.DiPietro.ea:16, , Section 2.4); see also Boffi.DiPietro: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 DiPietro.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:

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

Highorder 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 elementbyelement as in standard finite elements. The construction is conceived so that only facebased DOFs are globally coupled, which paves the way to efficient practical implementations where elementbased 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 residualbased a posteriori estimates for the energynorm 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 DiPietro.Droniou:16 ; DiPietro.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 DiPietro.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 DiPietro.Droniou.ea:15 , in Section 5 we extend the HHO method to diffusionadvectionreaction problems. In this context, a crucial property from the numerical point of view is robustness in the advectiondominated 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 firstorder 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 advectivereactive bilinear form. The HHO discretization is finally obtained in Section 5.3 combining the diffusive and advectivereactive contributions, and its stability with respect to an energylike norm including an advective derivative contribution is studied. In Section 5.4 we state an energynorm 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.
Definition 1 (Polytopal mesh)
A polytopal mesh of is a couple where:
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 .
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 (DiPietro.Ern:12, , Chapter 1), and refer to the case of isotropic meshes with nondegenerate faces. Isotropic means here that we do not consider the case of elements that become more and more stretched when refining. Nondegenerate 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 (DiPietro.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
(1) 
Morevover, according to (DiPietro.Ern:12, , Lemma 1.42), for all , all , and all
(2) 
Discrete functional analysis results for arbitraryorder methods on regular mesh sequences can be found in (DiPietro.Ern:12, , Chapter 1) and DiPietro.Droniou:16 ; DiPietro.Droniou:17 . We also refer the reader to Eymard.Gallouet.ea:00 for a first theorization of discrete functional analysis in the context of lowestorder 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 vectorvalued 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 scalarvalued, variate polynomials of total degree . We note the following trace inequality (see (DiPietro.Ern:12, , Lemma 1.46)): There is a real number only depending on , , and such that, for all , all , all , and all ,
(3) 
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
(4) 
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 vectorvalued projector denoted by and obtained by applying componentwise. The following boundedness result is a special case of (DiPietro.Droniou:16, , Corollary 3.7): For any , there exists a real number depending only on , , , and such that, for all , all , and all ,
(5) 
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
(6a)  
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  
(6b) 
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 DiPietro.Droniou:16 ): For any and , there exists a real number depending only on , , , , and such that, for all , all , and all ,
(7a)  
and, if ,  
(7b)  
where . 
3 Basic principles of Hybrid HighOrder methods
To fix the main ideas and notation, we study in this section the HHO discretization of the Poisson problem: Find such that
(8a)  
(8b) 
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 DiPietro.Ern.ea:16 .
The starting point to devise a HHO discretization is the following weak formulation of problem (8): Find such that
(9) 
where the bilinear form is such that
(10) 
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 :
(11) 
Specializing (11) to , we obtain
(12a)  
where we have used (6) to insert into the lefthand side and (4) to insert and into the righthand side after observing that and for all . Moreover, recalling (6b) and using the definition (4) of the projector, we infer that  
(12b) 
The relations (12) show that computing the elliptic projection does not require a full knowledge of the function . All that is required is

, 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 ;

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):
(13) 
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 ,
(14) 
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
(15) 
Potential reconstruction operator
Inspired by formula (12), we introduce the potential reconstruction operator such that, for all ,
(16a)  
and  
(16b) 
Notice that is a polynomial function on one degree higher than the elementbased DOFs . By definition, for all it holds that
(17) 
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
(18) 
where the first term in the righthand 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.DiPietro:16 :
Assumption 1 (Local stabilization bilinear form )
The local stabilization bilinear form satisfies the following properties:

Symmetry and positivity. is symmetric and positive semidefinite;

Stability. There is a real number independent of and of , but possibly depending on , , and , such that
(19) 
Polynomial consistency. For all and all , it holds that
(20)
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 ,
(21) 
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
(22) 
This choice, inspired by the Virtual Element literature BeiraodaVeiga.Brezzi.ea:13 , differs from the original HHO stabilization of DiPietro.Ern.ea:14 , where the following expression is considered instead:
(23) 
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 (DiPietro.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 HighOrder method of DiPietro.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.DiPietro.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 (DiPietro.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
(24) 
Proof
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
(25) 
Recalling (14), we have that
(26) 
Using the boundedness of resulting from (5) with , and followed by the optimal approximation properties (7a) of (with , , , and ), it is inferred that
(27) 
On the other hand, for all it holds that
(28)  
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 singlevalued interface unknowns:
Notice that singlevalued 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 ,
(29) 
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 (DiPietro.Droniou:16, , Proposition 5.4): There exists a real number independent of , but possibly depending on , , and , such that, for all ,
(30) 
Proposition 2 (Norm )
The map defines a norm on .
Proof
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 elementbyelement assembly setting, for all ,
(31) 
Lemma 1 (Properties of )
The bilinear form enjoys the following properties:

Stability. For all it holds with as in (19) that
(32) 
Consistency. There is a real number independent of , but possibly depending on , , and , such that, for all