Let \(\Omega_h\) be a collection of disjoint elements that partition \(\Omega\). The shape of the elements is not important in this general framework. Moreover, \(\Omega_h\) needs not to be conforming. An interior ''face'' of \(\Omega_h\) is any set \(F\) of positive \((n-1)\)-Lebesgue measure of the form \(F = \partial K^+\cap\partial K^-\) for some two elements \(K^+\) and \(K^-\) of \(\Omega_h\). We say that \(F\) is a boundary face if there is an element \(K\in\Omega_h\) such that \(F = \partial K\cap\Gamma\) and the \((n-1)\)-Lebesgue measure of \(F\) is not zero. Let \(\mathcal E^o_h\) and \(\mathcal E^\partial_h\) denote the set of interior and boundary faces of \(\Omega_h\), respectively. We denote by \(\mathcal E_h\) the union of all the faces in \(\mathcal E^o_h\) and \(\mathcal E^\partial_h\).

The global finite element spaces for the approximated flux \(\mathbf u_h\) and scalar solutions \(p_h\) are

We also need to introduce the spaces

Different choices of the local spaces \(\mathbf V(K), W(K)\), and \(M(F)\) correspond to different hybrid methods.

For any discontinuous (scalar or vector) function \(u\) in \(W_h\) or \(\mathbf V_h\), the trace \(u|_F\) on an interior face \(F = \partial K^+\cap\partial K^-\) is a double value function, whose two branches are denoted by \(u|_{K^+}\) and \(u|_{K^-}\). Here, \(\mathbf n_K\) denotes the unit outward normal of \(K\). For any double-valued vector function \(\mathbf v\), we define the jump of its normal component across an interior face \(F\) by

On any face \(F\) of \(K\) lying on the boundary, we set

A similar notation will be used for scalar functions. For functions \(u\) and \(v\) in \(L^2(D)\), we write \((u, v)_D = \int_D uv\) if \(D\) is a domain of \(\mathbb R^n\), and \(\langle u, v\rangle_D = \int_Duv\) if \(D\) is a domain of \(\mathbf R^{n-1}\). Also, we will use the notation

This method is obtained by using the Raviart-Thomas method to define the local solvers. The three ingredients of the RT-H method are:

The accuracy of the RT-H method is summarized in section [accuracy]. Note that, because \([[ \hat{\mathbf u}_h ]]\) and test functions \(\mu\) belong to the same space <Sayas 2013>, conservativity condition \((13)\) forces

so the normal component of the numerical trace \(\hat{\mathbf u}_h\) is single-valued, and \(\mathbf u_h\in H(\text{div},\Omega)\).

This method is obtained by using the Brezzi-Douglas-Marini method to define the local solvers. The three ingredients of the BDM-H method are:

Everything said about the RT-H method in the previous subsection applies to the BDM-H method.

The spaces of RT-H and BDM-H can be balanced to have equal polynomial degree. Stability is restored using a discrete stabilization (not penalization) function. The resulting method is known as the Hybridizable Discontinuous Galerkin (HDG) method. The HDG method is obtained by using the local DG method to define the local solvers. The three ingredients of the HDG method are:

The function \(\tau\) can be double valued on \(\mathcal E_h^o\), with two branches \(\tau^-=\tau_{K^-}\) and \(\tau^+=\tau_{K^+}\) defined on the face \(F\) shared by the finite elements \(K^+\) and \(K^-\). Note that the numerical trace of the flux \(\hat{\mathbf u}_h\) (but not the flux itself, as \(\tau_K\ne 0\)) is conservative. The accuracy of the HDG method is summarized in section [accuracy].

This section is mainly based on <Fu 2013>. As stated before, the goal of hybridization is the reduction (or static condensation) of the system \((16)-(20)\) to a linear system where only \(\hat p_h\) shows up. The remaining two variables \(\mathbf u_h\) and \(p_h\) will be reconstructed after solving for \(\hat p_h\), in an element-by-element fashion, easy to realize due to the fact that equations \((16)\) and \((17)\) are local or, in other words, the spaces \(\mathbf V_h\) and \(W_h\) are completely discontinous. In this section we will show how to perform static condensation on the linear system obtained by using the HDG method. This procedure can be easily adapted to other hybrid methods. Let us recall that the HDG method looks for an approximate solution \((\mathbf u_h, p_h, \hat p_h)\) in the space \(\mathbf V_h\times W_h\times M_h\) satisfying the equations

for all \((\mathbf v, w, \mu_1, \mu_2, \mu_3)\in\mathbf V_h\times W_h\times M_h^o\times M_h^N\times M_h^D\).

The following table summarizes the effect of the local spaces and the stabilization parameter \(\tau\) on the accuracy of the method on simplexes. We denote by \(\overline p_h|_K\) the integral average of \(p_h\) on \(K\in\Omega_h\). For the HDG method, the superconvergence of \(\overline p_h\) is what allows to get a solution of enhanced accuracy by postprocessing.

We introduce here a class of hybridizable methods able to use different local solvers in different elements and to easily handle nonconforming meshes. To define these methods, we need to specify the numerical fluxes, the local finite element spaces, and the space of approximate traces:

the function \(\tau_K\) is allowed to change on \(\partial K\). . The local space \(\mathbf V(K)\times W(K)\) can be any of the following:

Here, if \(F = \partial K^+\cap\partial K^-\), we set \(k(F) := \max\{k(K^+), k(K^-)\}\). For each element \(K\in\Omega_h\) and each face \(F\in\mathcal E_h\) on \(\partial K\), we take \(\tau_K|_F\in[0,\infty)\) and

Choice \((16)\) allows to deal with the nonconformity of the mesh in a very natural way. Also, the choice \(\tau_K = \infty\) could be allowed provided that the definition of the local solvers is modified as in <Cockburn 2009>.

The main features of this class of methods are:

For other possible methods, see <Cockburn 2009>.