documentation pending 
1. HDG methods for diffusionreaction
Most of the material presented in this section can be found in (Cockburn 2009).
We provide a first characterization of HDG methods for the following secondorder elliptic model problem:
Here \(\Omega\subset\mathbb R^n\) is a polyhedral domain \((n\geq 2)\), \(\partial\Omega = \Gamma = \Gamma_D \cup \Gamma_N\), \(d(\mathbf x)\) is a scalar nonnegative function, \(\Lambda(\mathbf x)\) is a matrix valued function that is symmetric and uniformly positive definite on \(\Omega\), and \(f\in L^2(\Omega)\). If \(d = 0\), we get the Darcy model presented in [darcymodel]. These assumptions can be generalized but we take them to simplify the discussion.
The Darcy problem presented in [stabilizedformulations] fits into this framework. Here, rather than using stabilized Galerkin formulations, the idea is to insert a Lagrange multiplier to handle the continuity of the normal components of the approximated flux \(\mathbf{u}_h\). In other words, the requirement \(\mathbf{u}_h\in H(\text{div},\Omega)\) will be written in weak form initially, and then recovered in strong form later. The introduction of a Lagrange multipliers leads to a formulation with three fields (Boffi 2013). We will show how the interior fields can be eliminated in order to build a discrete system that has lost the saddle point structure and only contains degrees of freedom on the faces.
2. The general structure of the methods
2.1. Notation
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 \((n1)\)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 \((n1)\)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 doublevalued 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^{n1}\). Also, we will use the notation
2.2. Reaching the formulation
2.2.1. A characterization of the exact solution
Two functions \(\mathbf u, p\) are exact solutions of \((1)(4)\) if and only if they satisfy the local problems \(\) \begin{align*} \Lambda\mathbf u + \nabla p &= \mathbf 0 & &\text{in }K,\\ \nabla\cdot\mathbf u + d p &= f & &\text{in K}, \end{align*} \(\) the transmission conditions \(\) \begin{align*} [[ \hat p ]] &= 0 & &\text{if }F\in\mathcal E^o_h,\\ [[ \hat{\mathbf u}]] &= 0 & &\text{if }F\in\mathcal E^o_h, \end{align*} \(\) and the boundary conditions \(\) \begin{align} \hat p &= h_D & &\text{if }F\in\Gamma_D,\\ \hat{\mathbf u}\cdot\mathbf n &= h_N & &\text{if }F\in\Gamma_N, \end{align} \(\) where \(\hat p\) and \(\hat{\mathbf u}\) are the traces of \(p\) and \(\mathbf u\), respectively. The transmission conditions state that the normal component of \(\hat{\mathbf u}\) across interelement boundary must be continuous. We can obtain \((\mathbf u, p)\) in \(K\) in terms of \(\hat p\) on \(\partial K\) and \(f\) by solving the following local Dirichlet problem \(\) \begin{align*} \Lambda\mathbf u + \nabla p &= \mathbf 0 & &\text{in }K, & &(5)\\ \nabla\cdot\mathbf u + d p &= f & &\text{in K}, & &(6)\\ p &= \hat p & &\text{on }\partial K. & &(7) \end{align*} \(\) The function \(\hat p\) can now be determined as the solution, on each \(F\in\mathcal E_h\), of the global equations \(\) \begin{align*} [[ \hat{\mathbf u} ]]_F &= 0 & &\text{if }F\in\mathcal E^o_h, & &(8)\\ [[ \hat{\mathbf u} ]]_F &= h_N & &\text{if }F\in\Gamma_N, & & (9)\\ \hat p_F &= h_D & &\text{if }F\in\Gamma_D. & &(10) \end{align*} \(\)
Hybrid methods are obtained by constructing discrete versions of \((5)(10)\). In this way, the globally coupled degrees of freedom will be those of the corresponding global formulations.
2.2.2. The local solvers: a weak formulation on each element
From \((5)(7)\), on the element \(K\in\Omega_h\), we define \((\mathbf u_h, p_h)\) in terms of \((\hat p_h, f)\) as the element of \(\mathbf V(K)\times W(K)\) such that
for all \((\mathbf v, w) \in \mathbf V(K)\times W(K)\), where \(\hat{\mathbf u}_h\) is the numerical trace of the flux and depends, in general, on \(\mathbf u_h, p_h\), and \(\hat p_h\). Different methods will correspond to different choices of \(\hat{\mathbf u}_h\).
2.2.3. The global problem
From \((8)(10)\), We determine \(\hat p_h\in M_h\) by requiring that
By solving \((11), (12)\) for \((\mathbf u_h, p_h)\) in terms of \((\hat p_h, f)\) at each element and plugging the results into \((13)(15)\), we get a system whose globally coupled degrees of freedom are those of the numerical trace \(\hat p_h\). This procedure corresponds to performing static condensation on the full discrete global system written in terms of \(\mathbf u_h, p_h, \hat p_h\).
If the (extension by zero to \(\mathcal E_h\) of the) function \([[ \hat{\mathbf u}_h ]]_{\mathcal E_h^o}\) belongs to the space \(M_h\), then condition \((13)\) is stating that \([[ \hat{\mathbf u}_h ]]_{\mathcal E_h^o} = 0\) pointwise, that is, the normal component of the numerical trace \(\hat{\mathbf u}_h\) is single valued. This means that the function \(\hat{\mathbf u}_h\) is a conservative numerical flux (\(\hat{\mathbf u}_h\in H(\text{div},\Omega)\)).
2.2.4. Summary
The approximate solution \((\mathbf u_h, p_h, \hat p_h)\) is the element of the space \(\mathbf V_h\times W_h\times M_h\) satisfying the equations
where the local spaces \(\mathbf V(K), W(K), M(F)\), as well as the numerical trace \(\hat{\mathbf q}_h\), need to be specified.
3. Examples of hybridizable methods
In this section we give som examples of methods fitting the general structure described in the previous section. The first three methods use the same local solver in all the elements \(K\) of the mesh \(\Omega_h\) and assume that \(\Omega_h\) is a conforming simplicial mesh. The fourth example is a class of methods employing different local solvers in different parts of the domain, which can easily deal with nonconforming meshes. To define each method, we have only to specify:

the numerical trace of the flux \(\hat{\mathbf u}_h\);

the local spaces \(\mathbf V(K), W(K)\);

the space of approximate traces \(M_h\).
3.1. The RTH method
This method is obtained by using the RaviartThomas method to define the local solvers. The three ingredients of the RTH method are:

\(\hat{\mathbf u}_h = \mathbf u_h\) on \(\partial K\), for each \(K\in\Omega_h\).

\(\mathbf V(K) = [P_k(K)]^n + \mathbf x P_k(K),\quad W(K) = P_k(K),\quad k\geq 0\).

\(M_h = \{\mu\in L^2(\mathcal E_h) : \mu_F\in P_k(F)\quad\forall\,F\in\mathcal E_h\}\).
The accuracy of the RTH 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 singlevalued, and \(\mathbf u_h\in H(\text{div},\Omega)\).
3.2. The BDMH method
This method is obtained by using the BrezziDouglasMarini method to define the local solvers. The three ingredients of the BDMH method are:

\(\hat{\mathbf u}_h = \mathbf u_h\) on \(\partial K\), for each \(K\in\Omega_h\).

\(\mathbf V(K) = [P_k(K)]^n,\quad W(K) = P_{k1}(K),\quad k\geq 1\).

Same \(M_h\) of the RTH method.
Everything said about the RTH method in the previous subsection applies to the BDMH method.
3.3. The HDG method
The spaces of RTH and BDMH 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:

For each \(K\in\Omega_h\): \(\hat{\mathbf u}_h = \mathbf u_h + \tau_K(p_h  \hat p_h)\mathbf n\quad\text{on }\partial K,\)
where \(\tau_K\) is a nonnegative function that can vary on \(\partial K\), and \(\tau_K > 0\) on at least one face of \(\partial K\). 
\(\mathbf V(K) = [P_k(K)]^n,\quad W(K) = P_k(K),\quad k\geq 0\).

Same \(M_h\) of the RTH method.
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].
3.3.1. Enhanced accuracy by postprocessing
The approximate solution and flux of the HDG method can be locally postprocessed to enhance their accuracy (Cockburn 2010).

Postprocessing of the scalar variable:
if we look for \(p_h^*:\Omega\to\mathbb R\) such that \(p_h^*_K\in P_{k+1}(K)\) and for all \(K\in\Omega_h\)
then it can be shown that this local postprocessed approximation has one additional order of convergence.

Postprocessing of the flux:
we can obtain a postprocessed flux \(\mathbf u_h^*\) with better conservation properties. Although \(\mathbf u_h^*\) converges at the same order as \(\mathbf u_h\), it is in \(H(\text{div},\Omega)\) and its divergence converges at one higher order than \(\mathbf u_h\). On each \(K\in\Omega_h\), we take \(\mathbf u_h^* :=\mathbf u_h + \boldsymbol\eta_h\) where \(\boldsymbol\eta_h\) is the only element of \([P_k(K)]^n + \mathbf x P_k(K)\) satisfying
3.4. Hybridization in matrix form
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 elementbyelement 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\).
3.4.1. Local solvers
Introduce the matrices related to the local bilinear forms
If \(\hat p_h\in M_h\) is known, equations \((21), (22)\) are uniquely solvable for \(\mathbf u_h, p_h\)and can be solved elementbyelement. Let us represent \(\mathbf u_h_K, p_h_K\), and \(\hat p_h_{\partial K}\) with vectors \(\mathbf u_K, \mathbf p_K\), and \(\mathbf p_{\partial K}\), respectively. Also, let
Then, the matrix representation of the local solutions is
Let us define
The flux prescribed by the HDG method
creates a bilinear form
whose matrix representation is (using \((26)\))
with
3.4.2. Boundary conditions and global solver

Dirichlet boundary conditions. The discrete Dirichlet boundary conditions \((25)\) require finding the projection \(\mathbf{\hat p}_D\) of the function \(h_D\) on the space \(M_h_{\Gamma_D}\).

Neumann boundary conditions. Neumann boundary conditions will appear in the right hand side of the global system.

Assemblying the global solver. The local solvers produce matrices \(D^K\) that need to be assembled to get a global matrix \(\mathbb H\). This matrix collects the fluxes \((27)\) from all the elements, with the result that opposing sign fluxes in internal faces (the normal vector points in different directions) are added. The matrices \(D^K_f\) also have to be assembled to get a global vector \(\mathbf F\). At this point, the global system reads
where \(\mathbf G_N\) is the vector containing the elements of \(\langle h_N, \mu\rangle_{\Gamma_N}, \mu\in M_h_{\Gamma_N}\) in the degrees of freedom corresponding to Neumann faces and zeros everywhere else. What is left is the elimination of Dirichlet degrees of freedom from \((28)\), namely, values of Dirichlet faces are taken from \(\mathbf{\hat p}_D\) and sent to the right hand side of the system, and rows corresponding to Dirichlet degrees of freedom are ignored.
3.5. Orders of accuracy for RTH, BDMH, HDG
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.
Method 
\(\tau\) 
\(\mathbf u_h\) 
\(p_h\) 
\(\overline p_h\) 
\(k\) 
RTH 
\(0\) 
\(k+1\) 
\(k+1\) 
\(k+2\) 
\(\geq 0\) 
BDMH 
\(0\) 
\(k+1\) 
\(k\) 
\(k+2\) 
\(\geq 2\) 
HDG 
\(O(h)\) 
\(k+1\) 
\(k\) 
\(k+2\) 
\(\geq 1\) 
HDG 
\(O(1)\) 
\(k+1\) 
\(k+1\) 
\(k+2\) 
\(\geq 1\) 
HDG 
\(O(1)\) 
\(1\) 
\(1\) 
\(1\) 
\(=0\) 
HDG 
\(O(1/h)\) 
\(k\) 
\(k+1\) 
\(k+1\) 
\(\geq 1\) 
3.6. A class of hybridizable methods well suited for adaptivity
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:

For any simplex \(K\in\Omega_h\), we take
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:

\(([P_{k(K)}(K)]^n + \mathbf x P_{k(K)}(K)) \times P_{k(K)}(K)\), where \(k(K)\geq0\) and \(\tau_K\geq 0\) on \(\partial K\),

\([P_{k(K)}(K)]^n \times P_{k(K)1}(K)\), where \(k(K)\geq1\) and \(\tau_K\geq 0\) on \(\partial K\),

\([P_{k(K)}(K)]^n \times P_{k(K)}(K)\), where \(k(K)\geq0\) and \(\tau_K > 0\) on at least one face \(F\in\partial K\).

The space of approximate traces is

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:

Variable degree approximation spaces on conforming meshes. The RTH, BDMH, and HDM methods considered above use a single local solver in each of the elements \(K\) of the conforming triangulation \(\Omega_h\). A variabledegree version of each of these methods is a particular case of the clas of methods presented here.

Automatic coupling of different methods on conforming meshes. The class presented here allows for the use of different local solvers in different elements \(K\in\Omega_h\), which are then automatically coupled.

Mortaring capabilities (for nonconforming meshes). This class incorporate a mortaring ability thanks to the form that the numerical trace of the flux on \(\partial K\) takes on an interior face \(F\in\mathcal E_h^o\), and thanks to the definition of the stabilization parameter. Let us give an example. If we have a conforming mesh, we can take the first choice of local spaces (2a) and set \(\tau = 0\). The resulting method is nothing but the RTH method. We can easily modify this method to handle nonconforming meshes by simply taking \(\tau_K\in(0,\infty)\) on every \(F\in\mathcal E_h^o\) which is not a face of \(K\), and otherwise, taking \(\tau_K 0\).
For other possible methods, see (Cockburn 2009).
The goal of this section is to test the convergence of the HDG algorithm for the
mixed poisson problem in 2D and 3D domain.
The exact solutions and domain are presented in the following figures
for 2D and 3D.
We have tested three cases:

only Dirichlet conditions:

only Neumann conditions:

with integral boundary condition:
For those tests, we use the following parameters:

\(k=1\)

\(p=\frac{1}{2\pi}atan(y,x)\)

\(u=\frac{1}{2\pi(x^2+y^2)}(y,x)\)

\(f=0\)
On the following domain: