The aim of this chapter is to introduce to
hybridization of discontinuous Galerkin methods
within the \Rheolef\  environment.
For a review of hybridizable discontinuous Galerkin methods,
see~\cite{NguPerCoc-2011}.
The hybridization technique allows an efficient finite element
implementation of many problems of importance,
such as the Navier-Stokes one.
Let us start by some model problems.

% -----------------------------------------------------------------------
\subsection{The Poisson problem}
% -----------------------------------------------------------------------
Let us consider the Poisson problem with mixed 
Dirichlet and Neumann boundary conditions:

\ \ $(P)$: \emph{find $u$, defined in $\Omega$, such that}
\begin{eqnarray*}
  -\Delta u &=& f 
	\ \mbox{ in } \ \Omega
	\\
  u &=& g_d
	\ \mbox{ on } \ \Gamma_d
	\\
  \Frac{\partial u}{\partial n} &=& g_n
	\ \mbox{ on } \ \Gamma_n
\end{eqnarray*}
where $\partial \Omega=\Gamma_d\cup\Gamma_n$,
and the interior of the two boundary domains 
$\Gamma_d$ and $\Gamma_n$ are disjoints.
The data $f$, $g_d$ and $g_n$ are given.
Let us introduce the gradient $\boldsymbol{\sigma}=\nabla u$ as an independent variable.
The problem writes equivalently (see e.g.~\cite{NguPerCoc-2011}):

\ \ $(P)$: \emph{find $\boldsymbol{\sigma}$ and $u$, defined in $\Omega$, such that}
\begin{subequations}
\begin{empheq}[left={\empheqlbrace}]{align}
  \boldsymbol{\sigma} - \nabla u &= 0 
	\ \mbox{ in } \ \Omega
	\label{eq-dirichlet-hdg-sigma}
	\\
  {\rm div}(\boldsymbol{\sigma}) &= -f 
	\ \mbox{ in } \ \Omega
	\label{eq-dirichlet-hdg-u}
	\\
  u &= g_d
	\ \mbox{ on } \ \Gamma_d
	\label{eq-dirichlet-hdg-bc-d}
	\\
  \boldsymbol{\sigma}.\boldsymbol{n} &= g_n
	\ \mbox{ on } \ \Gamma_n
	\label{eq-dirichlet-hdg-bc-n}
\end{empheq}
\end{subequations}
Let us multiply \eqref{eq-dirichlet-hdg-sigma} by a test function $\boldsymbol{\tau}$
and integrate by parts on any element $K$:
\begin{eqnarray*}
  \int_K
    \left(
      \boldsymbol{\sigma}.\boldsymbol{\tau}
      +
      {\rm div}(\boldsymbol{\tau})
      \, u
    \right)
    \,{\rm d}x
  -
  \int_{\partial K}
    (\boldsymbol{\tau}.\boldsymbol{n})
    \, u
    \,{\rm d}s
  &=&
  0 
\end{eqnarray*}
In the discontinuous Galerkin method, the trace of $u$ on $\partial K$ will
be discontinuous across the inter-element boundaries.
The hybridization replaces this trace by a new independent variable,
denoted as $\lambda$,
that is only definied on sides of the mesh:
\begin{subequations}
\begin{eqnarray}
  \int_K
    \left(
      \boldsymbol{\sigma}.\boldsymbol{\tau}
      +
      {\rm div}(\boldsymbol{\tau})
      \, u
    \right)
    \,{\rm d}x
  -
  \int_{\partial K}
    (\boldsymbol{\tau}.\boldsymbol{n})
    \, \lambda
    \,{\rm d}s
  &=&
  0 
  \label{eq-dirichlet-hdg-fv-sigma}
\end{eqnarray}
The $\lambda$ variable will be uni-valued on boundary inter-elements
and accounts strongly for the boundary condition:
$\lambda=g_d$ on $\Gamma_d$.
%
Then, let us multiply \eqref{eq-dirichlet-hdg-u} by a test function $v$
with a zero average value on~$K$:
\begin{eqnarray*}
  \int_K
    {\rm div}(\boldsymbol{\sigma})
    \,v
    \,{\rm d}x
  &=&
  -
  \int_K
    f 
    \,v
    \,{\rm d}x
\end{eqnarray*}
In order to weakly impose the condition $u=\lambda$ on $\partial K$,
we add a penalization term:
\begin{eqnarray}
  \int_K
    {\rm div}(\boldsymbol{\sigma})
    \,v
    \,{\rm d}x
  -
  \beta
  h^n
  \int_{\partial K}
    h^n
    (u-\lambda)
    \,v
    \,{\rm d}s
  &=&
  -
  \int_K
    f 
    \,v
    \,{\rm d}x
  \label{eq-dirichlet-hdg-fv-u}
\end{eqnarray}
\end{subequations}
Here, $h$ denotes the local mesh size in the element~$K$.
The two constants~$\beta > 0$ and~$n\in\mathbb{R}$
are respectively a penalization coefficient
and power index.
Following~\cite{CocGuzWan-2009}, we consider the three cases~$n=0$,
$1$ and~$-1$.
We have three unknowns $\boldsymbol{\sigma}$, $u$ and $\lambda$
and only the two
equations~\eqref{eq-dirichlet-hdg-fv-sigma}
      and~\eqref{eq-dirichlet-hdg-fv-u}.
For the problem to be complete, we add
an equation for~$\lambda$.
Observe that~\eqref{eq-dirichlet-hdg-fv-u} writes equivalently,
after a second integration by part on~$K$:
\begin{eqnarray*}
  -
  \int_K
    \boldsymbol{\sigma}.\nabla v
    \,{\rm d}x
  +
  \int_{\partial K}
    \left(
      \boldsymbol{\sigma}.\boldsymbol{n}
      -
      \beta 
      h^n
      (u-\lambda)
    \right)
    v
    \,{\rm d}s
  &=&
  -
  \int_K
    f 
    \,v
    \,{\rm d}x
\end{eqnarray*}
\cindex{numerical flux}%
The quantity involved in the integral on~$\partial K$
is denoted as
\mbox{$
  \widehat{\boldsymbol{\sigma}}
  =
  \boldsymbol{\sigma}
  -
  \beta h^n (u-\lambda)
  \,\boldsymbol{n}
$}
and refered as the \emph{numerical flux}
across inter-element boundaries.
As an additional equation, we impose the continuity
of the normal component of this numerical flux.
On the $\Gamma_n$ boundary domain, the normal component
of the numerical flux is imposed to be 
the prescribed flux $g_n$, while the Dirichlet condition 
$\lambda=g_d$ is precribed on the $\Gamma_d$ boundary domain:
\begin{subequations}
\begin{eqnarray}
  \jump{
    \boldsymbol{\sigma}
    -
    \beta h^n (u-\lambda)
    \,\boldsymbol{n}
  }
  .\boldsymbol{n}
  &=&
  0
  \ \mbox{ on }\ S
  ,\ \forall S\in \mathscr{S}_h^{(i)}
  \label{eq-dirichlet-hdg-flux-internal}
  \\
  \boldsymbol{\sigma}
  .\boldsymbol{n}
  -
  \beta h^n (u-\lambda)
  &=&
  g_n
  \ \mbox{ on }\ \Gamma_n
  \\
  \lambda
  \label{eq-dirichlet-hdg-flux-boundary-neumann}
  &=&
  g_d
  \ \mbox{ on }\ \Gamma_d
  \label{eq-dirichlet-hdg-flux-boundary-dirichlet}
\end{eqnarray}
\end{subequations}
%TODO: cross ref with main volume
%Let us consider the notation of
%    Fig.~\ref{fig-dg-neighbor},
%page~\pageref{fig-dg-neighbor}.
For an internal side
\mbox{$
  S = \partial K_+ \cap \partial K_- \in\mathscr{S}_h^{(i)}
$}
between two elements 
\mbox{$
  K_1, K_2 \in \mathscr{T}_h
$},
relation~\eqref{eq-dirichlet-hdg-flux-internal} writes:
\begin{eqnarray*}
  \boldsymbol{\sigma}_+
  .
  \boldsymbol{n}_+
  -
  \beta h^n (u_+-\lambda)
  +
  \boldsymbol{\sigma}_-
  .
  \boldsymbol{n}_-
  -
  \beta h^n (u_--\lambda)
  &=&
  0
  \ \mbox{ on }\ S
\end{eqnarray*}
where 
$\boldsymbol{\sigma}_\pm$ and
$u_\pm$ are the trace on~$S$ of
$\boldsymbol{\sigma}$ and~$u$
in~$K_\pm$
and
$\boldsymbol{n}_\pm$ are the
outer normal of~$K_\pm$ on~$S$.
Since~$S$ is oriented, let us choose without loss of generality
\mbox{$
 \boldsymbol{n}
 =
 \boldsymbol{n}_-
 =
 -
 \boldsymbol{n}_+
$}.
\cindex{operator!\code{average}, across sides}%
\cindex{operator!\code{jump}, across sides}%
Then, the previous relation writes:
\begin{eqnarray*}
  &&
  \left(
    \boldsymbol{\sigma}_-
    -
    \boldsymbol{\sigma}_+
  \right)
  .
  \boldsymbol{n}
  -
  \beta h^n (u_- + u_+)
  +
  2\beta  h^n \lambda
  \ =\ 
  0
  \ \mbox{ on }\ S
  \\
  \Longleftrightarrow
  &&
  \jump{\boldsymbol{\sigma}}.\boldsymbol{n}
  -
  2\beta h^n (\average{u}-\lambda)
  \ =\ 
  0
  \ \mbox{ on }\ S
\end{eqnarray*}
where we have used the jump and average across the internal side~$S$.
The previous relation, together
with~\eqref{eq-dirichlet-hdg-flux-boundary-neumann}
 and~\eqref{eq-dirichlet-hdg-flux-boundary-dirichlet},
writes in variational form:
\begin{eqnarray}
  \int_{\mathscr{S}_h^{(i)}}
    \left(
      \jump{\boldsymbol{\sigma}}.\boldsymbol{n}
      -
      2\beta h^n
      (\average{u}-\lambda)
    \right)
    \, \mu
    \,{\rm d}s
  +
  \int_{\partial\Omega}
    \left(
      \boldsymbol{\sigma}.\boldsymbol{n}
      -
      \beta
      h^n
      (u-\lambda)
    \right)
    \, \mu
    \,{\rm d}s
  &=&
  \int_{\Gamma_n}
    g_n
    \, \mu
    \,{\rm d}s
  \label{eq-dirichlet-hdg-fv-lambda}
\end{eqnarray}
for all test function $\mu$, defined on all
internal sides of $\mathscr{S}_h^{(i)}$
and on all sides of the boundary domain $\Gamma_n$
and that vanishes on all sides of the boundary domain $\Gamma_d$.

Grouping~\eqref{eq-dirichlet-hdg-fv-sigma},
         \eqref{eq-dirichlet-hdg-fv-u}
     and~\eqref{eq-dirichlet-hdg-fv-lambda},
we obtain the discrete variational formulation:

\ \ $(FV)_h$: \emph{find $(\boldsymbol{\sigma}_h,u_h,\lambda_h)\in T_h\times X_h\times \Lambda_h(g_d)$ such that}
\begin{eqnarray*}
  \int_\Omega
     \boldsymbol{\sigma}
    .\boldsymbol{\tau}
    \,{\rm d}x
  +
  \int_\Omega
    {\rm div}_h(\boldsymbol{\tau})
    \, u
    \,{\rm d}x
  -
  \int_{\mathscr{S}_h^{(i)}}
    (\jump{\boldsymbol{\tau}}.\boldsymbol{n})
    \, \lambda
    \,{\rm d}s
  -
  \int_{\partial\Omega}
    (\boldsymbol{\tau}.\boldsymbol{n})
    \, \lambda
    \,{\rm d}s
  &=&
  0 
  \\
  \int_\Omega
    {\rm div}_h(\boldsymbol{\sigma})
    \,v
    \,{\rm d}x
  -
  \beta
  \sum_{K\in\mathscr{T}_h}
  \int_{\partial K}
    h^n
    u\,v
    \,{\rm d}s
  +
  \int_{\mathscr{S}_h^{(i)}}
    2\beta
    h^n
    \average{v}
    \,\lambda
    \,{\rm d}s
  +
  \int_{\partial\Omega}
    \beta
    h^n
    v
    \,\lambda
    \,{\rm d}s
  &=&
  -
  \int_\Omega
    f 
    \,v
    \,{\rm d}x
  \\
  \int_{\mathscr{S}_h^{(i)}}
    \left(
      2\beta
      h^n
      \average{u}
      -
      \jump{\boldsymbol{\sigma}}.\boldsymbol{n}
    \right)
    \mu
    \,{\rm d}s
  +
  \int_{\partial\Omega}
    \left(
      \beta
      h^n
      u
      -
      \boldsymbol{\sigma}.\boldsymbol{n}
    \right)
    \mu
    \,{\rm d}s
  -
  \int_{\mathscr{S}_h^{(i)}}
    2\beta
    h^n
    \lambda
    \, \mu
    \,{\rm d}s
  -
  \int_{\partial\Omega}
    \beta
    h^n
    \lambda
    \, \mu
    \,{\rm d}s
  &=&
  -
  \int_{\Gamma_n}
    g_n
    \, \mu
    \,{\rm d}s
\end{eqnarray*}
for all $(\boldsymbol{\tau}_h,v_h,\tau_h)\in T_h\times X_h\times \Lambda_h(0)$,
where
\begin{eqnarray*}
  T_h
  &=&
  \left\{
    \boldsymbol{\tau} \in \left(L^2(\Omega)\right)^{d}
    \ ;\ 
    \boldsymbol{\tau}_{|K} \in (P_k)^d
    ,\ \forall K\in \mathscr{T}_h
  \right\}
  \\
  X_h
  &=&
  \left\{
    v\in L^2(\Omega)
    \ ;\ 
    v_{|K} \in P_k
    ,\ \forall K\in \mathscr{T}_h
  \right\}
  \\
  M_h
  &=&
  \left\{
    \mu\in L^2(\mathscr{S}_h)
    \ ;\ 
    \mu_{|S} \in P_k
    ,\ \forall S\in \mathscr{S}_h
  \right\}
  \\
  \Lambda_h(g_d)
  &=&
  \left\{
    \mu\in M_h
    \ ;\ 
    \mu_{|S} = \pi_h(g_d)
    ,\ \forall S\subset \Gamma_d
  \right\}
\end{eqnarray*}
and $k\geq 0$ is the polynomial degree.
This is a symmetric system with a mixed structure.
Let:
\begin{eqnarray*}
  a_h(
    (\boldsymbol{\sigma},u),
    (\boldsymbol{\tau},v))
  &=&
  \int_\Omega
     \left(
       \boldsymbol{\sigma}
      .\boldsymbol{\tau}
      +
      {\rm div}_h(\boldsymbol{\tau})
      \, u
      +
      {\rm div}_h(\boldsymbol{\sigma})
      \,v
    \right)
    \,{\rm d}x
  -
  \beta
  \sum_{K\in\mathscr{T}_h}
  \int_{\partial K}
    h^n
    u\,v
    \,{\rm d}s
  \\
  b_h((\boldsymbol{\tau},v),\mu)
  &=&
  \int_{\mathscr{S}_h^{(i)}}
    \left(
      -
      \jump{\boldsymbol{\tau}}.\boldsymbol{n}
      +
      2\beta
      h^n
      \average{v}
    \right)
    \mu
    \,{\rm d}s
  +
  \int_{\partial\Omega}
    \left(
      -
      \boldsymbol{\tau}.\boldsymbol{n}
      +
      \beta
      h^n
      \average{v}
    \right)
    \mu
    \,{\rm d}s
  \\
  c_h(\lambda,\mu)
  &=&
  \int_{\mathscr{S}_h^{(i)}}
    2\beta
    h^n
    \lambda
    \, \mu
    \,{\rm d}s
  +
  \int_{\partial\Omega}
    \beta
    h^n
    \lambda
    \, \mu
    \,{\rm d}s
  \\
  \ell_h(\boldsymbol{\tau},v)
  &=&
  -
  \int_\Omega
    f 
    \,v
    \,{\rm d}x
  \\
  k_h(\mu)
  &=&
  -
  \int_{\Gamma_n}
    g_n
    \, \mu
    \,{\rm d}s
\end{eqnarray*}
Then, the variational formulation writes equivalently:

\ \ $(FV)_h$: \emph{find $(\boldsymbol{\sigma}_h,u_h,\lambda_h)\in T_h\times X_h\times \Lambda_h(g_d)$ such that}
\begin{eqnarray*}
  a_h(
    (\boldsymbol{\sigma},u),
    (\boldsymbol{\tau},v))
  +
  b_h((\boldsymbol{\tau},v),\lambda)
  &=&
  \ell_h(\boldsymbol{\tau},v)
  \\
  b_h((\boldsymbol{\sigma},u),\mu)
  -
  c_h(\lambda,\mu)
  &=&
  k_h(\mu)
\end{eqnarray*}
for all $(\boldsymbol{\tau}_h,v_h,\tau_h)\in T_h\times X_h\times \Lambda_h(0)$.
%
Let $\chi_h=(\boldsymbol{\sigma}_h,u_h)$.
\cindex{matrix!Schur complement}%
This linear symmetric system admits the following matrix structure:
\begin{eqnarray*}
  &&
  \left( \begin{array}{cc}
    A & B^T \\
    B & -C
  \end{array} \right)
  \left( \begin{array}{c}
    \boldsymbol{\chi}_h \\
    \lambda_h
  \end{array} \right) 
  \ =\ 
  \left( \begin{array}{c}
    F \\
    G
  \end{array} \right) 
\end{eqnarray*}
A carreful study shows that the $A$ matrix
has a block-diagonal structure at the element level
and can be easily inverted on the fly during the assembly process.
The matrix structure writes equivalently:
\begin{eqnarray*}
  &\Longleftrightarrow&
  \left\{ \begin{array}{rcl}
	A\boldsymbol{\chi}_h + B^T\lambda_h &=& F \\
	B\boldsymbol{\chi}_h -   C\lambda_h &=& G
  \end{array} \right.
  \\
  &\Longleftrightarrow&
  \left\{ \begin{array}{l}
	\boldsymbol{\chi}_h = A^{-1}(F-B^T\lambda_h) \\
	(C+BA^{-1}B^T)\lambda_h = BA^{-1}F - G
  \end{array} \right.
\end{eqnarray*}
The second equation is solved first:
the only remaining unknown is the Lagrange multiplier $\lambda_h$,
in a linear system involving the Schur complement matrix
\mbox{$
  S
  =
  C+BA^{-1}B^T
$}.
Then, the two variables $\boldsymbol{\chi}_h=(\boldsymbol{\sigma}_h,u_h)$
are simply obtained by a direct computation.

% -------------------------------------
\myexamplelicense{dirichlet_hdg.cc}

\myexamplenoinput{sinusprod_dirichlet.h}

\myexamplenoinput{sinusprod_error_hdg.cc}
% -------------------------------------

\begin{figure}[htb]
  \begin{tabular}{lll}
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-p1-error-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-00-error-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-m1-error-l2-t.pdf} \\
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-p1-error-sigma-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-00-error-sigma-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-m1-error-sigma-l2-t.pdf} 
% TODO: lambda_h - proj_l2(lambda) : super-cv k+2 : cf CocGuzWan-2009 table 3
%   \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-00-error-lambda-l2-t.pdf} &
  \end{tabular}
  \caption{Hybrid discontinuous Galerkin method: 
	Poisson problem with Dirichlet boundary conditions in 2D geometry.
        Convergence vs $h$ and $k$ for the approximation of the solution~$u_h$
          and of its gradient~$\boldsymbol{\sigma}_h$.
        (left) $n=1$~; (center) $n=0$~; (right) $n=-1$.
  }
  \label{fig-dirichlet-hdg-error} 
\end{figure}
The right-hand-side~$f$ and the Dirichlet boundary condition~$g$
has been chosen such that the exact solution is given by:
\begin{equation*}
  u(x) = \prod_{i=0}^{d-1} \sin(\pi x_i)
\end{equation*}
The files~\file{sinusprod_dirichlet.h}, that defines the data~$f$ and~$g$,
and~\file{sinusprod_error_hdg.cc}, that compute the error
in various norms, are not listed here but are available in the \Rheolef\  example directory.

\subsubsection*{How to run the program}
The compilation and run write:
\begin{verbatim}
  make dirichlet_hdg
  mkgeo_grid -t 10 > square.geo
  ./dirichlet_hdg square.geo P1d > square.field 
  field square.field -elevation
  field square.field -elevation -mark lambda
  field square.field -velocity  -mark sigma
  make sinusprod_error_hdg
  ./sinusprod_error_hdg < square.field 
\end{verbatim}
Fig.~\ref{fig-dirichlet-hdg-error} plots the errors vs $h$ and $k$
for a two dimensional geometry.
Observe that the approximation~$u_h$
converges optimaly, as $h^{k+1}$, in the~$L^2$ norm
for any $k\geq 0$ when the power index~$n=1$.
When~$n=1$, the convergence is suboptimal, as~$h^k$ only and
note that the lowest order approximation~$k=0$ is not convergent.
When~$n=-1$, the convergence is optimal, as~$h^{k+1}$,
only when~$k\geq 1$,
while the lowest order approximation~$k=0$ is not convergent.
The approximation~$\boldsymbol{\sigma}_h$ of the gradient
converges optimaly, as~$h^{k+1}$, for both~$n=0$ and~$1$,
while it is suboptimal, as~$h^k$, when~$n=-1$.
All these results are consistent with the approxiamtion theory
of the HDG method, see~\cite{CocGuzWan-2009} and tables~2, 3, 6.

% TODO: lambda_h - proj_l2(lambda) : super-cv k+2 : cf CocGuzWan-2009 table 3
% The Lagrange multiplier~$\lambda_h$ converges also to~$u$
% on the boundaries of the elements as~$h^{k+1}$
% in both the $L^2$ and $L^\infty$ norms.

% -----------------------------------------------------------------------
\subsection{Superconvergence of the Lagrange multiplier}
% -----------------------------------------------------------------------
\cindex{convergence!error!superconvergence}%

\begin{figure}[htb]
  %\begin{center}
  \begin{tabular}{lll}
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-p1-error-lambda-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-00-error-lambda-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-m1-error-lambda-l2-t.pdf}
  \end{tabular}
  %\end{center}
  \caption{Hybrid discontinuous Galerkin method:
	Poisson problem with Dirichlet boundary conditions in 2D geometry.
        Super-convergence vs $h$ and $k$ for the Lagrange multiplier~$\lambda_h$
        to the $L^2$ projection on~$M_h$ of the exact solution.
        (left) $n=1$~; (center) $n=0$~; (right) $n=-1$.
  }
  \label{fig-dirichlet-hdg-multiplier-error}
\end{figure}
Let~$\pi_{M_h}$ denote the~$L_2$ projection
from~$L^2(\mathscr{S}_h)$ into~$M_h$.
We consider the special case of computing
the projection~$\pi_{M_h}(u)\in M_h$ of the restriction
to~$\mathscr{S}_h$ of an element~$u\in L^2(\mathscr{S}_h)$.
It is defined as 
\begin{eqnarray*}
  \pi_{M_h}(u)
  &=&
  \arginf_{\mu_h\in M_h}
  \ 
  \Frac{1}{2}
  \sum_{S\in\mathscr{S}_h}
  \int_{S}
    (u - \mu_h)^2
    \mathrm{d}s
\end{eqnarray*}
The following bilinear forms are introduced:
\begin{eqnarray*}
  m_s(\lambda,\mu)
  &=&
  \sum_{S\in\mathscr{S}_h}
  \int_{S}
    \lambda \, \mu
    \,\mathrm{d}s
  \\
  k_h(\mu)
  &=&
  \sum_{S\in\mathscr{S}_h}
  \int_{S}
    u \, \mu
    \,\mathrm{d}s
\end{eqnarray*}
Then, the projection reduces to:

\ \ \emph{find $\bar{\lambda}_h=\pi_{M_h}(u)\in M_h$ such that}
\begin{eqnarray*}
  m_s(\bar{\lambda}_h,\mu_h)
  &=&
  k_h(\mu_h)
  ,\ \ \forall \mu_h\in M_h
\end{eqnarray*}
The result of projection operator can be obtained
by a resolution of a linear system.
\cindex{norm!mesh-dependent on {$\mathscr{S}_h$}}%
Following~\cite[p.~1896]{CocDonGuz-2008-b},
in order to measure the error on the set of sides~$\mathscr{S}_h$ of the mesh,
we introduce the mesh-dependent norm~$\|.\|_{0,2,\mathscr{S}_h}$, defined
for all~$\mu\in L^2(\mathscr{S}_h)$ by
\[
  \|\mu\|_{0,2,\mathscr{S}_h}^2
  =
  \sum_{S\in\mathscr{S}_h}
  \int_{S}
    h_S
    \mu^2
    \,\mathrm{d}s
\]
where~$h_S$ is a characteristic length on the side $S\in\mathscr{S}_h$.
%
Fig.~\ref{fig-dirichlet-hdg-multiplier-error} plots the difference
\mbox{$
  \pi_{M_h}(u) - \lambda_h
$}
in this mesh-dependent norm.
In agrement with the theoretical results~\cite[p.~1896]{CocDonGuz-2008-b},
we observe the superconvergence of the multiplier~$\lambda_h$
to the $L^2$ projection~$\pi_{M_h}(u)$.
More precisely, when~$n=1$ the order of convergence is~$k+2$
for any~$k\geq 0$.
When~$n=0$, the superconvergence occurs only when~$k\geq 1$
and when~$n=-1$ there is no superconvergence.
This error is also computed by the code~\file{sinusprod_error_hdg.cc}.

% -----------------------------------------------------------------------
\subsection{Superconvergence of the piecewise averaged solution}
% -----------------------------------------------------------------------
\cindex{convergence!error!superconvergence}%

% -----------------------------------------------------------------------
\myexamplelicense{dirichlet_hdg_average.icc}

\myexamplelicense{dirichlet_hdg_average.cc}
% -----------------------------------------------------------------------


\begin{figure}[htb]
  %\begin{center}
  \begin{tabular}{lll}
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-p1-average-error-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-00-average-error-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-m1-average-error-l2-t.pdf}
  \end{tabular}
  %\end{center}
  \caption{Hybrid discontinuous Galerkin method:
	Poisson problem with Dirichlet boundary conditions in 2D geometry.
        Super-convergence vs $h$ and $k$ for the piecewise average~$\bar{u}_h$ of the approximate solution.
        (left) $n=1$~; (center) $n=0$~; (right) $n=-1$.
  }
  \label{fig-dirichlet-hdg-average-error}
\end{figure}

% -------------------------------------
\myexamplenoinput{sinusprod_error_hdg_average.cc}
% -------------------------------------
The superconvergence of the piecewise average values of the
solution obtained by the hybrid discontinuous Galerkin method
was first observed in~\cite{CocDonGuz-2008-b,CocGuzWan-2009}
and then exploited for many postprocessing application (see e.g.~\cite{NguPerCoc-2011}).
%
The file \file{dirichlet_hdg_average.icc}
compute~$\bar{u}_h=\bar{\pi}_h(u_h,\lambda_h)$,
the averaged solution, defined by (see~\cite{CocDonGuz-2008-b}, eqn~(2.9b)):
\[
  \bar{\pi}_h(u_h,\lambda_h)
  =
  \left\{ \begin{array}{ll}
     \dfrac{1}{d}
     \displaystyle
     \sum_{S\subset\partial K} \lambda_h 
     & \textrm{when } k=0
     \\
     &\\
     \dfrac{1}{\mathrm{meas}(K)} 
     \displaystyle
     \int_K u_h\,\mathrm{d}x
     & \textrm{when } k\geq 1
  \end{array} \right.
\]
The file~\file{sinusprod_error_hdg_average.cc} that compute the error
\mbox{$\bar{u}_h-\bar{\pi}_h(u)$}
is not listed here but is available in the \Rheolef\  example directory.
The computation of the error is obtained by:
\begin{verbatim}
  make dirichlet_hdg dirichlet_hdg_average sinusprod_error_hdg_average
  mkgeo_grid -t 10 > square.geo
  ./dirichlet_hdg square.geo P1d | dirichlet_hdg_average | \
      ./sinusprod_error_hdg_average
\end{verbatim}
Observe on Fig.~\ref{fig-dirichlet-hdg-average-error}
that, when~$n=1$ and for any~$k\geq 0$,
we obtain~\mbox{$\bar{u}_h-\bar{\pi}_h(u)$} of order~$k+2$.
This is a remarkable result since~$u_h$ is piecewise polynomial
of order~$k$ and is expected to converge at order~$k+1$.
When~$n=0$, for any~$k\geq 1$ we also observe this superconvergence
while, when~$k=0$, the average value converges only
at first order.
Finaly, when~$n=-1$, for any~$k\geq 1$ we do no more observe
These observations are consistent with those of
\cite[p.~16]{CocGuzWan-2009}
(see also~\cite[p.~1613]{CocDonGuz-2008-b}).

%
% -----------------------------------------------------------------------
\subsection{Improving the solution by local postprocessing}
% -----------------------------------------------------------------------
\label{subsec-hdg-postprocessing}%
\cindex{convergence!postprocessing}%
By combining the approximation~$\boldsymbol{\sigma}_h$ of the gradient~$\nabla u$,
that converges at rate~$k+1$
  (see Fig.~\ref{fig-dirichlet-hdg-error}),
with the average~$\bar{u}_h$ that super-converges at rate~$k+2$
  (see Fig.~\ref{fig-dirichlet-hdg-average-error}),
it is then possible, by a local integration inside each element,
to obtain a new approximation~$u_h^*$,
that is piecewise discontinuous polynomial of order~$k+1$,
and that converges at rate~$k+2$.

The postprocess step is nothing than the resolution of
the following local Neumann problem in any element~$K\in\mathscr{T}_h$
(see~\cite[p.~1360]{CocGopSay-2010}, eqn~(5.1)):

\ \ $(P^*)$: \emph{find $u^*$, defined in $K$, such that}
\begin{empheq}[left={\empheqlbrace}]{align*}
  -\Delta u^*
	&= f
	\ \mbox{ in } K
	\\
  \Frac{\partial u^*}{\partial n}	
	&= \boldsymbol{\sigma}_h.\boldsymbol{n}
	\ \mbox{ on } \partial K
	\\
  \int_K u^* \,\mathrm{d}x
  	&=
  	\int_K \bar{u}_h \,\mathrm{d}x
\end{empheq}
where $\boldsymbol{\sigma}_h$ and~$\bar{u}_h$ are given
by the previous resolution.

For any $\beta\in\mathbb{R}$,
let us introduce the following functional space:
\begin{eqnarray*}
  X^*_K(\beta)
  &=&
  \left\{
    v \in H^1(K)
    \ ;\ 
    \int_K v \,\mathrm{d}x = \beta
  \right\}
\end{eqnarray*}
Then, the variational formulation of the problem writes:

%ICI
\ \ $(FV^*)$: \emph{find
  $u^* \in X_K\left( \int_K \bar{u}_h \,\mathrm{d}x\right)$
  such that}
\begin{eqnarray*}
  \int_K
    \nabla u^* .\nabla v^*
    \,{\rm d}x
  &=&
  \int_K
    f \, v^*
    \,{\rm d}x
  +
  \int_{\partial K}
    \boldsymbol{\sigma}_h.\boldsymbol{n}
    \, v^*
    \,{\rm d}s
  ,\ \ \forall v^*\in X_K(0)
\end{eqnarray*}
%
The linear constraint for the imposition of the
average value is not easy to impose in~$X_K(\beta)$.
Indeed, $X_K(\beta)$ is not a vectorial space when~$\beta\neq 0$.
Following the methodology previously introduced
%TODO: cross ref with main volume
% in section~\ref{sec-neumann-laplace},
%   page~\pageref{sec-neumann-laplace},
for the Neumann boundary conditions for the Laplace operator,
we introduce a Lagrange multiplier
denoted here as~$\zeta$,
which is constant inside each element~$K$.

Finally, the postprocessing step compute the
approximation~$u^*_h$,
when~$\boldsymbol{\sigma}_h$ and~$u_h$ are known,as:

\ \ $(FV^*)_h$: \emph{find
  $(u^*_h,\zeta_h)
  \in
  X^*_h \times Z_h$
  such that,
  on each element $K\in\mathscr{T}_h$, we have}
\begin{empheq}[left={\empheqlbrace}]{align*}
  \int_K
    \nabla u_h^*.\nabla v_h^*
    \, {\rm d}x
  +
  \int_K
    v_h^*
    \,\zeta_h
    \, {\rm d}x
  &=
  \int_K
    f \, v_h^*
    \, {\rm d}x
  +
  \int_{\partial K}
    \boldsymbol{\sigma}_h.\boldsymbol{n}
    \, v_h^*
    \,{\rm d}s
  ,\ \ \forall v_h^*\in X_h^*
  \\
  \int_K
    u^*_h
    \,\xi_h
    \, {\rm d}x
  \phantom{aaaaaaaaaaaaaa}
  &=
  \int_K
    \bar{u}_h
    \,\xi_h
    \, {\rm d}x
  ,\ \ \forall \xi_h\in X_h^*
\end{empheq}
where
\begin{eqnarray*}
  X^*_h
  &=&
  \left\{
    v\in L^2(\Omega)
    \ ;\ 
    v_{|K} \in P_{k+1}
    ,\ \forall K\in \mathscr{T}_h
  \right\}
  \\
  Z_h
  &=&
  \left\{
    \xi\in L^2(\Omega)
    \ ;\ 
    \xi_{|K} \in P_{0}
    ,\ \forall K\in \mathscr{T}_h
  \right\}
\end{eqnarray*}
Let
\begin{eqnarray*}
  a^*_h 
  (u,\zeta;\,
   v,\xi)
  &=&
  \int_\Omega
    \left(
      \nabla_h u.\nabla_h v
      +
      v \,\zeta
      +
      u \,\xi
    \right)
    \, {\rm d}x
  \\
  \ell^*_h (v,\xi)
  &=&
  \int_\Omega
    \left(
      f\, v
      +
      \bar{u}_h
      \,\xi
    \right)
    \, {\rm d}x
  +
  \sum_{K\in \mathscr{T}_h}
  \int_{\partial K}
    (\boldsymbol{\sigma}_h.\boldsymbol{n})
    \, v
    \, {\rm d}s
\end{eqnarray*}
The second step writes in this abstract setting:

\ \ $(FV^*)_{h}$: \emph{find
\mbox{$
  (u^*_h,\zeta^*_h)
  \in
  X^*_h \times Z_h
$}
such that}
\begin{eqnarray*}
  a^*_h 
  (u^*_h,\zeta_h;\,
   v_h^*,\xi_h)
  &=&
  \ell^*_h (v_h^*,\xi_h)
  ,\ \ \forall (v_h^*,\xi_h) \in X^*_h \times Z_h
\end{eqnarray*}
Note that the matrix associated to the bilinear form~$a^*_h$
is symmetric and block-diagonal:
it can thus be easily be inverted on the fly at the element level.
The present postprocessing stage was first 
introduced in~\cite[p.~1360]{CocGopSay-2010}, eqn~(5.1)
as a replacement and a simplification of
those previously introduced in~\cite{CocDonGuz-2008-b,CocGuzWan-2009}.
The following code implement this postprocessing.

\myexamplelicense{dirichlet_hdg_post.cc}

\begin{figure}[htb]
  \begin{tabular}{lll}
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-p1-post-error-us-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-00-post-error-us-l2-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-m1-post-error-us-l2-t.pdf}
  \end{tabular}
  \caption{Solution post-processing of the hybrid discontinuous Galerkin method:
	Poisson problem with Dirichlet boundary conditions in 2D geometry.
        Convergence vs $h$ and $k$ for the post-treated approximate solution.
        (left) $n=1$~; (center) $n=0$~; (right) $n=-1$.
  }
  \label{fig-dirichlet-hdg-post-error} 
\end{figure}

\subsubsection*{How to run the program}
The compilation and run write:
\begin{verbatim}
  make dirichlet_hdg dirichlet_hdg_post ./sinusprod_error_hdg
  mkgeo_grid -t 10 > square.geo
  ./dirichlet_hdg square.geo P1d > square.field 
  ./dirichlet_hdg_post < square.field > square-post.field
  field square-post.field -elevation
  ./sinusprod_error_hdg < square.field
  ./sinusprod_error_hdg < square-post.field
\end{verbatim}
The results are shown on Fig.~\ref{fig-dirichlet-hdg-post-error}.
When~$n=1$, observe that, for any~$k\geq 0$,
the error for the post-treated solution~$u^*_h$ converges to zero at rate~$k+2$
in~$L^2$ norm, which is optimal,
since~$u^*_h$ is a piecewise $k+1$ degree polynomial.
When~$n=0$, this result is obtained only for~$k\geq 1$
while, when~$n=-1$, the convergence is suboptimal.

% -----------------------------------------------------------------------
\subsection{Improving the gradient with the Raviart-Thomas element}
% -----------------------------------------------------------------------
\label{subsec-hdg-post-raviart-thomas}%
\apindex{Raviart-Thomas}%
Let $(\boldsymbol{\sigma}_h,u_h,\lambda_h)$ be the solution
of problem~$(FV)_h$.
This solution is considered here as known.
%
Observe that the normal component of the numerical flux
\mbox{$
  \widehat{\boldsymbol{\sigma}}_h
  .\boldsymbol{n}
  =
  \boldsymbol{\sigma}_h
  .\boldsymbol{n}
  -
  \beta h^n (u_h-\lambda_h)
$}
is continuous across any inter-element boundaries.
Indeed, relation~\eqref{eq-dirichlet-hdg-flux-internal}
is imposed weakely in the variational formulation~$(FV)_h$,
with a Lagrange multiplier~$\mu_h\in M_h$
which is piecewise polynomial of degree~$k$.
Since~$\widehat{\boldsymbol{\sigma}}_h$ is also
piecewise polynomial of degree~$k$,
the discrete version of~\eqref{eq-dirichlet-hdg-flux-internal}
is also satisfied strongly:
\begin{eqnarray*}
  \widehat{\boldsymbol{\sigma}}_h
  .\boldsymbol{n}
  \ = \ 
  \jump{
    \boldsymbol{\sigma}_h
    -
    \beta h^n (u_h-\lambda_h)
    \,\boldsymbol{n}
  }
  .\boldsymbol{n}
  &=&
  0
  \ \mbox{ on }\ S
  ,\ \forall S\in \mathscr{S}_h^{(i)}
\end{eqnarray*}
Since any element of $H({\rm div},\Omega)$ presents the
continuity of its normal component 
across any inter-element boundaries,
is possible to define a new $H({\rm div},\Omega)$ conform
approximation of $\boldsymbol{\sigma}$,
denoted as $\tilde{\boldsymbol{\sigma}}_h$,
which satisfies
\mbox{$
  \widetilde{\boldsymbol{\sigma}}_h
  .\boldsymbol{n}
  \ = \ 
  \widehat{\boldsymbol{\sigma}}_h
  .\boldsymbol{n}
$}
on any internal sides.
Similarly to~\eqref{eq-dirichlet-hdg-flux-internal},
$\tilde{\boldsymbol{\sigma}}_h$ also should
satisfy a discrete version
of~\eqref{eq-dirichlet-hdg-flux-boundary-neumann}-\eqref{eq-dirichlet-hdg-flux-boundary-dirichlet} i.e.
\mbox{$
  \widehat{\boldsymbol{\sigma}}_h
  .\boldsymbol{n}
$}
is equals to~$g_n$ on~$\Gamma_n$ and to
\mbox{$
  \boldsymbol{\sigma}_h
  .\boldsymbol{n}
  -
  \beta  h^n (u_h-g_d)
$}
on~$\Gamma_d$.
It is characterized
as the unique element
\mbox{$
  \tilde{\boldsymbol{\sigma}}_h \in \tilde{T}_h
$}
satisfying in any element $K\in\mathscr{T}_h$
the following local variational formulation~\cite[p.~1360]{CocGopSay-2010}:
\begin{eqnarray*}
  \int_K
    \tilde{\boldsymbol{\sigma}}_h
    .\tilde{\boldsymbol{\gamma}}_h
    \, {\rm d}x
  &=&
  \int_K
    \boldsymbol{\sigma}_h
    .\tilde{\boldsymbol{\gamma}}_h
    \, {\rm d}x
  ,\ \ \forall \tilde{\boldsymbol{\gamma}}_h \in \tilde{G}_h
  \\
  \int_{\partial K}
    \left(
      \tilde{\boldsymbol{\sigma}}_h.\boldsymbol{n}
    \right)
    \tilde{\mu}_h
    \, {\rm d}s
  &=&
  \int_{\partial K}
    \left(
      \boldsymbol{\sigma}_h.\boldsymbol{n}
      -
      \beta  h^n (u_h - \lambda_h)
    \right)
    \tilde{\mu}_h
    \, {\rm d}s
  ,\ \ \forall \tilde{\mu}_h \in \tilde{M}_h
\end{eqnarray*}
where we have introduced the finite dimensional spaces
\begin{eqnarray*}
  \tilde{T}_h
  &=&
  \left\{
    \tilde{\boldsymbol{\tau}}_h \in \left( L^2(\Omega)\right)^d ;\ 
    \tilde{\boldsymbol{\tau}}_{h|K} \in RT_{k}(K)
    ,\ \forall K \in \mathscr{T}_h
  \right\}
  \\
  \tilde{W}_h
  &=&
  \left\{
    \tilde{\boldsymbol{\gamma}}_h \in \left( L^2(\Omega)\right)^d ;\ 
    \tilde{\boldsymbol{\gamma}}_{h|K} \in P_{k-1}(K)
    ,\ \forall K \in \mathscr{T}_h
  \right\}
  \\
  \tilde{M}_h
  &=&
  \displaystyle
  \prod_{K \in \mathscr{T}_h}
  \left\{
    v_{|\partial K} \ ;\ 
    v\in P_{k}(K)
  \right\}
\end{eqnarray*}
\cindex{approximation!Raviart-Thomas}%
The space~$\tilde{T}_h$
contains discontinuous and piecewise
$k$-th order Raviart-Thomas~$RT_k$ polynomial vector-valued functions.
\apindex{discontinuous!trace}%
The space~$\tilde{M}_h$ contains, inside each element,
the normal trace of functions of~$\tilde{T}_h$
while~$\tilde{W}_h$ contains piecewise discontinuous
polynomial of degree~$k\!-\!1$.
By convention, when~$k=0$, then the space~$\tilde{W}_h$ is empty.
Observe that
\mbox{$
  \mathrm{dim}(\tilde{T}_h)
  =
  \mathrm{dim}(\tilde{W}_h)
  +
  \mathrm{dim}(\tilde{M}_h)
$}
and that the previous local variational problem is well-posed.
Let
\begin{eqnarray*}
  \tilde{a}_h (\tilde{\boldsymbol{\sigma}}_h ;\ 
               [\tilde{\boldsymbol{\gamma}}_h, \tilde{\mu}_h])
  &=&
  \int_\Omega
    \tilde{\boldsymbol{\sigma}}_h
    .\tilde{\boldsymbol{\gamma}}_h
    \, {\rm d}x
  +
  \sum_{K\in\mathscr{T}_h}
  \sum_{S\subset \partial K}
  \int_{S}
    (\tilde{\boldsymbol{\sigma}}_h.\boldsymbol{n})
    \,\tilde{\mu}_h
    \, {\rm d}s
  \\
  \tilde{\ell}_h ([\tilde{\boldsymbol{\gamma}}_h, \tilde{\mu}_h])
  &=&
  \int_\Omega
    \boldsymbol{\sigma}_h
    .\tilde{\boldsymbol{\gamma}}_h
    \, {\rm d}x
  +
  \sum_{K\in\mathscr{T}_h}
  \int_{\partial K}
    \left(
      \boldsymbol{\sigma}_h.\boldsymbol{n}
      -
      \beta  h^n (u_h - \lambda_h)
    \right)
    \,\tilde{\mu}_h
    \, {\rm d}s
\end{eqnarray*}
For known
$(\boldsymbol{\sigma}_h,u_h,\lambda_h)$,
the postprocessing of the gradient writes:

\ \ $(\widetilde{FV})_h$: \emph{find $\tilde{\boldsymbol{\sigma}}_h\in \tilde{T}_h$ such that}
\mbox{$
  \tilde{a}_h (\tilde{\boldsymbol{\sigma}}_h,
               [\tilde{\boldsymbol{\gamma}}_h,\tilde{\mu}_h])
  =
  \tilde{\ell}_h ([\tilde{\boldsymbol{\gamma}}_h,\tilde{\mu}_h])
$}
for all $(\tilde{\boldsymbol{\gamma}}_h,\tilde{\mu}_h) \in \tilde{W}_h\times \tilde{M}_h$.

A carreful study shows that the matrix
associated with the~$\tilde{a}_h$ bilinear form
has a block-diagonal structure at the element level
and can be efficiently inverted on the fly during the assembly process.
This property is due to the usage of discontinuous
Raviart-Thomas approximation space~$\tilde{T}_h$.
Moreover, since the normal component
of $\tilde{\boldsymbol{\sigma}}_h$ is
equal to the numerical flux 
\mbox{$
  \widehat{\boldsymbol{\sigma}}_h.\boldsymbol{n}
$}
which is continuous across inter-element boundaries
this is also the case for
$\tilde{\boldsymbol{\sigma}}_h.\boldsymbol{n}$
and finally
$\tilde{\boldsymbol{\sigma}}_h \in H({\rm div},\Omega)$.

% -----------------------------------------
\myexamplelicense{dirichlet_hdg_post_rt.cc}

\myexamplenoinput{sinusprod_error_hdg_post_rt.cc}
% -----------------------------------------

\begin{figure}[htb]
  \begin{center}
  \begin{tabular}{lll}
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-p1-post-error-sigmat-hdiv-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-00-post-error-sigmat-hdiv-t.pdf} &
    \includegraphics[width=0.32\textwidth]{dirichlet-hdg-n-m1-post-error-sigmat-hdiv-t.pdf}
  \end{tabular}
  \end{center}
  \caption{Solution post-processing of the hybrid discontinuous Galerkin method:
	Poisson problem with Dirichlet boundary conditions in 2D geometry.
        Convergence vs $h$ and $k$ for the post-treated gradient~$\tilde{\boldsymbol{\sigma}}_h$ approximate solution.
        (left) $n=1$~; (center) $n=0$~; (right) $n=-1$.
  }
  \label{fig-dirichlet-hdg-post-rt-error} 
\end{figure}

\subsubsection*{How to run the program}
The compilation and run write:
\begin{verbatim}
  make dirichlet_hdg dirichlet_hdg_post_rt
  mkgeo_grid -t 10 > square.geo
  ./dirichlet_hdg square.geo P1d > square.field 
  ./dirichlet_hdg_post_rt < square.field > square-rt.field
  field square-rt.field -velocity -mark sigmat -proj P1
  make sinusprod_error_hdg_post_rt
  ./sinusprod_error_hdg_post_rt < square-rt.field
\end{verbatim}
Fig.~\ref{fig-dirichlet-hdg-post-rt-error} plots the errors vs $h$ and $k$
for a two dimensional geometry.
\cindex{norm!in $H(\mathrm{div})$}%
Observe that, when~$n=1$, the approximation~$\tilde{\boldsymbol{\sigma}}_h$
of the gradient converges to~$\nabla u$
with a~$k+1$ rate in~$H(\mathrm{div}$ norm, with:
\[
    \|\boldsymbol{\tau}\|_{\mathrm{div},2,\Omega}^2
    =
    \|\boldsymbol{\tau}\|_{0,2,\Omega}^2
    +
    \|\mathrm{div}(\boldsymbol{\tau})\|_{0,2,\Omega}^2
\]
for any~$\boldsymbol{\tau}\in H(\mathrm{div},\Omega)$.
It means that~$\mathrm{div}(\tilde{\boldsymbol{\sigma}}_h)$ converges to~$\Delta u$
with a~$k+1$ rate in both~$L^2$ norm.
This is optimal with respect to the classical interpolation results
for the $k$-th order Raviart-Thomas element
(see section~\ref{sec-raviart-thomas},
    page~\pageref{sec-raviart-thomas}).
When~$n=0$, the convergence is also optimal while, when~$n=-1$ it is suboptimal.

\begin{table}[htb]
  \begin{equation*}
    \begin{array}{|c||c|c||c|c|}
      \hline
      n & 
      u_h &
      \boldsymbol{\sigma}_h &
      u_h^* &
      \tilde{\boldsymbol{\sigma}}_h 
      \\ \hline \hline
      1 &
      k &
      k+1 &
      \red{k+2} &
      \red{k+1}
      \\ \hline
      0 &
      k+1 &
      k+1 &
      \left\{ \begin{array}{cl}
        1, & k=0 \\
        k+2, & k\geq 0
      \end{array} \right.
      &
      k+1
      \\ \hline
      -1 &
      \left\{ \begin{array}{cl}
        0, & k=0 \\
        k+1, & k\geq 0
      \end{array} \right.
      &
      k &
      \left\{ \begin{array}{cl}
        0, & k=0 \\
        k+1, & k\geq 0
      \end{array} \right.
      &
      k
      \\ \hline
    \end{array}
  \end{equation*}
  \caption{Hybrid discontinuous Galerkin method: 
        convergence order versus~$k$
        for~$n\in\{1,0,-1\}$.
  }
  \label{tab-dirichlet-hdg-summary} 
\end{table}
Fig.~\ref{tab-dirichlet-hdg-summary} summarizes the convergence orders versus
the mesh size~$h$ in terms of~$k$ and~$n$.
Observe that~$n=1$ allows one to obtain an optimal convergence
of both the the post-processed solution~$u_h^*$ and its gradient~$\tilde{\boldsymbol{\sigma}}_h$
for any polynomial degree~$k\geq 0$.
When~$n=0$, it is optimal only for~$k\geq 1$
while when~$n=-1$ it is suboptimal.
Thus the HDG method and its postprocessed stage could be considered as a whole.
It allows, by solving a linear system for the $k$-th piecewise polynomial
approximation of the Lagrange multiplier only, to obtain a~$k\!+\!2$ order approximation
of the solution in~$L^2$ and a~$k\!+\!1$ approximation of its gradient in~$H(\mathrm{div})$.