HDG methods for diffusion-reaction

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 second-order elliptic model problem:

\[\begin{align*} \Lambda\mathbf{u} + \nabla p &= \mathbf 0 & &\text{in }\Omega & &(1)\\ \nabla\cdot\mathbf u + d p &= f & &\text{in }\Omega & &(2)\\ p &= h_D & &\text{on }\Gamma_D & &(3)\\ \mathbf u\cdot\mathbf n &= h_N & &\text{on }\Gamma_N & &(4) \end{align*}\]

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 [darcy-model]. These assumptions can be generalized but we take them to simplify the discussion.

The Darcy problem presented in [stabilized-formulations] 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.

1. The general structure of the methods

1.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 \((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

\[\begin{align*} \mathbf V_h &= \left\{\mathbf v : \Omega\to\mathbb R^n : \mathbf v|_K \in\mathbf V(K)\quad\forall\,K\in\Omega_h\right\},\\ W_h &= \left\{w : \Omega\to\mathbb R : w|_K \in W(K)\quad\forall\,K\in\Omega_h\right\}. \end{align*}\]

We also need to introduce the spaces

\[\begin{align*} M_h &= \{\mu : \mathcal E_h\to \mathbb R : \mu|_F\in M(F)\quad\forall\,F\in\mathcal E_h \},\\ M_h^o &= \{\mu\in M_h : \mu|_\Gamma = 0 \},\\ M_h^D &= \{\mu\in M_h : \mu|_F = 0\quad\forall\,F\in\mathcal E_h\setminus\Gamma_D \},\\ M_h^N &= \{\mu\in M_h : \mu|_F = 0\quad\forall\,F\in\mathcal E_h\setminus\Gamma_N \}. \end{align*}\]

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

\[[[ \mathbf v ]]_F = \mathbf v_{K^+} \cdot\mathbf n_{K^+} + \mathbf v_{K^-}\cdot\mathbf n_{K^-}.\]

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

\[[[ \mathbf v ]]_F = \mathbf v_K \cdot\mathbf n_K.\]

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

\[(v, w)_{\Omega_h} = \sum_{K\in\Omega_h}(v, w)_K,\qquad\langle\mu,\lambda\rangle_{\partial\Omega_h} = \sum_{K\in\Omega_h}\langle\mu,\lambda\rangle_{\partial K}.\]

1.2. Reaching the formulation

1.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.

1.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

\[\begin{align*} (\Lambda\mathbf u_h,\mathbf v)_K - (p_h, \nabla\cdot\mathbf v)_K + \langle\hat p_h, \mathbf v\cdot\mathbf n\rangle_{\partial K} &= 0, & &(11)\\ -(\mathbf u_h, \nabla w)_K + \langle\hat{\mathbf u}_h\cdot\mathbf n,w\rangle_{\partial K} + (dp_h, w)_K &= (f,w)_K, & &(12) \end{align*}\]

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\).

1.2.3. The global problem

From \((8)-(10)\), We determine \(\hat p_h\in M_h\) by requiring that

\[\begin{align*} \langle\mu,[[ \hat{\mathbf u}_h ]]\rangle_F &= 0 & &\forall\,\mu\in M(F) & &\text{if }F\in\mathcal E_h^o, & &(13)\\ \langle\mu,[[ \hat{\mathbf u}_h ]]\rangle_F &= \langle\mu,h_N\rangle_F & &\forall\,\mu\in M(F) & &\text{if }F\in\Gamma_N, & &(14)\\ \langle\mu,\hat p_h\rangle_F &= \langle\mu,h_D\rangle_F & &\forall\,\mu\in M(F) & &\text{if }F\in\Gamma_D. & &(15) \end{align*}\]

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)\)).

1.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

\[\begin{align*} (\Lambda\mathbf u_h,\mathbf v)_{\Omega_h} - (p_h, \nabla\cdot\mathbf v)_{\Omega_h} + \langle\hat p_h, \mathbf v\cdot\mathbf n\rangle_{\partial\Omega_h} &= 0 & &\forall\mathbf v\in \mathbf V_h, & &(16)\\ -(\mathbf u_h, \nabla w)_{\Omega_h} + \langle\hat{\mathbf u}_h\cdot\mathbf n,w\rangle_{\partial\Omega_h} + (d p_h, w)_{\Omega_h} &= (f,w)_{\Omega_h} & &\forall w\in W_h, & &(17)\\ \langle\mu,\hat{\mathbf u}_h\cdot\mathbf n\rangle_{\partial\Omega_h\setminus\Gamma} &= 0 & &\forall \mu\in M^o_h, & &(18)\\ \langle\mu,\hat{\mathbf u}_h\cdot\mathbf n\rangle_{\Gamma_N} &= \langle\mu,h_N\rangle_{\Gamma_N} & &\forall\mu\in M^N_h, & &(19)\\ \langle\mu,\hat p_h\rangle_{\Gamma_D} &= \langle\mu,h_D\rangle_{\Gamma_D} & &\forall\mu\in M^D_h, & &(20) \end{align*}\]

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.

2. 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\).

2.1. The RT-H method

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

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

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

  3. \(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 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

\[[[ \hat{\mathbf u}_h]] = 0\quad\text{on }\mathcal E_h^o,\]

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

2.2. The BDM-H method

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:

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

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

  3. Same \(M_h\) of the RT-H method.

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

2.3. The HDG 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:

  1. 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\).

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

  3. Same \(M_h\) of the RT-H 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].

2.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\)

\begin{align} (\nabla p_h^*, \nabla w)_K &= -(\Lambda\mathbf u_h, \nabla w)_K & &\forall\,w\in P_{k+1}(K),\\ (p^*_h, 1)_K &= (p_h, 1)_K, & & \end{align}

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

\begin{align} (\boldsymbol\eta_h,\mathbf v)_K &= 0 & &\forall\,\mathbf v\in[P_{k-1}(K)]^n,\\ \langle\boldsymbol\eta_h\cdot\mathbf n, \mu\rangle_F &= \langle(\hat{\mathbf u}_h-\mathbf u_h)\cdot\mathbf n,\mu\rangle_F & &\forall\,F\in P_k(F),\quad\forall\,F\in\partial K. \end{align}

2.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 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

\[\begin{align*} &(\Lambda\mathbf u_h,\mathbf v)_{\Omega_h} & &- (p_h, \nabla\cdot\mathbf v)_{\Omega_h} & &+ \langle\hat p_h, \mathbf v\cdot\mathbf n\rangle_{\partial\Omega_h} & &= 0, & &(21)\\ &(\nabla\cdot\mathbf u_h, w)_{\Omega_h} & &+ \langle\tau p_h,w\rangle_{\partial\Omega_h} + (d p_h, w)_{\Omega_h} & &- \langle\tau \hat p_h,w\rangle_{\partial\Omega_h} & &= (f,w)_{\Omega_h}, & &(22)\\ &\langle\mathbf u_h\cdot\mathbf n,\mu_1\rangle_{\partial\Omega_h\setminus\Gamma} & &+ \langle\tau p_h,\mu_1\rangle_{\partial\Omega_h\setminus\Gamma} & &- \langle\tau \hat p_h,\mu_1\rangle_{\partial\Omega_h\setminus\Gamma} & &= 0, & &(23)\\ &\langle\mathbf u_h\cdot\mathbf n,\mu_2\rangle_{\Gamma_N} & &+ \langle\tau p_h,\mu_2\rangle_{\Gamma_N} & &- \langle\tau \hat p_h,\mu_2\rangle_{\Gamma_N} & &= \langle h_N,\mu_2\rangle_{\Gamma_N}, & &(24)\\ & & & & &\langle\hat p_h,\mu_3\rangle_{\Gamma_D} & &= \langle h_D,\mu_3\rangle_{\Gamma_D}, & &(25) \end{align*}\]

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\).

2.4.1. Local solvers

Introduce the matrices related to the local bilinear forms

\[\begin{align} A_{11}^K &\leftrightarrow (\Lambda\mathbf u_h,\mathbf v)_K, & &A_{12}^K\leftrightarrow- (p_h, \nabla\cdot\mathbf v)_K, & &A_{13}^K\leftrightarrow\langle\hat p_h, \mathbf v\cdot\mathbf n\rangle_{\partial K},\\ A_{21}^K &\leftrightarrow(\nabla\cdot\mathbf u_h, w)_K, & &A_{22}^K\leftrightarrow\langle\tau p_h,w\rangle_{\partial K} + (d p_h, w)_K, & &A_{23}^K\leftrightarrow\langle\tau \hat p_h,w\rangle_{\partial K},\\ A_{31}^K &\leftrightarrow\langle\mathbf u_h\cdot\mathbf n, \mu\rangle_{\partial K}, & &A_{32}^K\leftrightarrow\langle\tau p_h,\mu\rangle_{\partial K}, & &A_{33}^K\leftrightarrow\langle\tau \hat p_h,\mu\rangle_{\partial K},\\ & & &A_f^K\leftrightarrow (f,w)_K \end{align}\]

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 element-by-element. 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

\[\begin{align*} A^K &= \begin{bmatrix} A_{11}^K & A_{12}^K\\ A_{21}^K & A_{22}^K\\ \end{bmatrix}, & B^K &= \begin{bmatrix} A_{13}^K\\ A_{23}^K \end{bmatrix}, & F^K &= \begin{bmatrix} \mathbf 0\\ A_f^K \end{bmatrix}. \end{align*}\]

Then, the matrix representation of the local solutions is

\[\begin{align} &\begin{bmatrix} \mathbf u_K\\ \mathbf p_K \end{bmatrix} = -(A^K)^{-1}B^K \mathbf p_{\partial K} + (A^K)^{-1}F^K. & &(26) \end{align}\]

Let us define

\[C^K = \begin{bmatrix} A_{31}^K & A_{32}^K \end{bmatrix}.\]

The flux prescribed by the HDG method

\[\mathbf u_h\cdot\mathbf n + \tau(p_h-\hat p_h)\colon\partial K\to\mathbb R\]

creates a bilinear form

\[\mu\in M(\partial K)\to\langle\mathbf u_h\cdot\mathbf n + \tau (p_h-\hat p_h), \mu\rangle_{\partial K} = \langle\mathbf u_h\cdot\mathbf n + \tau p_h, \mu\rangle_{\partial K} - \langle \tau\hat p_h, \mu\rangle\]

whose matrix representation is (using \((26)\))

\[\begin{split} C^K\begin{bmatrix} \mathbf u_K\\ \mathbf p_K \end{bmatrix} - A_{33}^K\mathbf p_{\partial K} &= -C^K(A^K)^{-1}B^K \mathbf p_{\partial K} + C^K(A^K)^{-1}F^K - A_{33}^K\mathbf p_{\partial K}\qquad(27)\\ &= D_f^K - D^K\mathbf p_{\partial K}, \end{split}\]


\[\begin{align*} D_f^K &= C^K(A^K)^{-1}F^K, & D^K &= C^K(A^K)^{-1}B^K + A_{33}^K. \end{align*}\]

2.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

\begin{equation*} \mathbb H\,\mathbf{\hat p} = \mathbf F + \mathbf G_N,\qquad(28) \end{equation*}

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.

2.5. Orders of accuracy for RT-H, BDM-H, 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.



\(\mathbf u_h\)


\(\overline p_h\)







\(\geq 0\)






\(\geq 2\)






\(\geq 1\)






\(\geq 1\)












\(\geq 1\)

2.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:

  1. For any simplex \(K\in\Omega_h\), we take

\begin{equation*} \hat{\mathbf u}_h = \mathbf u_h + \tau_K(p_h - \hat p_h)\mathbf n\quad\text{on }\partial K, \end{equation*}

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\).

    1. The space of approximate traces is

\begin{equation*} M_h := \{\mu\in L^2(\mathcal E_h):\mu|_F\in P_{k(F)}\quad\forall\,F\in\mathcal E_h\}. \end{equation*}

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

\[\tau_K|_F\in(0,\infty)\quad\text{if }F\text{ is not a face of }K.\qquad(16)\]

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 RT-H, BDM-H, and HDM methods considered above use a single local solver in each of the elements \(K\) of the conforming triangulation \(\Omega_h\). A variable-degree 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 RT-H 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>.

3. References

  • [] Superconvergent discontinuous Galerkin methods for second-order elliptic problems B Cockburn, J Guzmán, H Wang - Mathematics of Computation, 2009

  • [[[ cockburnGS2010]]] The hybridizable discontinuous Galerkin methods, B Cockburn - Proceedings of the International Congress of …, 2010

  • [] Mixed finite element methods and applications, D Boffi, F Brezzi, M Fortin - 2013

  • [] Algorithm 949: MATLAB tools for HDG in three dimensions Z Fu, LF Gatica, FJ Sayas - ACM Transactions on Mathematical Software (TOMS), 2015

  • [] From Raviart-Thomas to HDG: a personal voyage, FJ Sayas - arXiv preprint arXiv:1307.2491, 2013