The aim of this chapter is to introduce to
the Raviart-Thomas element~\cite{RavTho-1977}
for building an approximation of the~$H(\mathrm{div},\Omega)$ space.

There is a subtle issue.
The \Rheolef\  implementation choice for this element
bases on internal interpolation nodes instead of moments.
This choice leads to more efficient computation of degrees
of freedom, but the standard Lagrange interpolation~$\pi_h$
no more satisfies the comutation diagram and
optimal error in the~$H(\mathrm{div},\Omega)$ norm.
Instead of the Lagrange interpolation~$\pi_h$, we propose
a projection operator, that satisfies these desired properties.
We start building this projection operator
for a piecewise discontinuous version of
this element: it allows one to build a projection
that requires only local operations and converges optimaly.
\index{approximation!discontinuous!hybrid}%
Moreover, the piecewise discontinuous Raviart-Thomas approximation
is used during the post-processing stage of
the hybrid discontinuous Galekin method,
that will be developped in a forthcoming chapter.

\red{TODO: it remains to merge degreees of freedom along sides
for obtaining a projector for the continuous Raviart-Thomas approximation,
and check that the obtained projection still satisfies
the commutation diagram and converges optimaly in~$H(\mathrm{div},\Omega)$.
Indeed, even with a $C^1(\bar{\Omega})$ function, it is not clear
that the obtained projection have degrees of freedom that matches
along internal sides.
}%red

Let
\begin{eqnarray*}
  V_h
  &=&
  \left\{
    \boldsymbol{v}_h \in \left(L^2(\Omega)\right)^d ; \ \ 
    \boldsymbol{v}_{h|K} \in RT_k(K), \ \forall K \in \mathcal{T}_h
  \right\}
  \\
  Q_h
  &=&
  \left\{ 
    q_h \in L^2(\Omega) ~; \ \ 
    q_{h|K}\in P_k, \ \forall K \in \mathcal{T}_h
  \right\}
\end{eqnarray*}
\apindex{discontinuous}%
\apindex{Raviart-Thomas}%
Here, $V_h$ represents the space of discontinuous and piecewise
$k$-th order Raviart-Thomas~$RT_k$ vector-valued functions
while~$Q_h$ is the space of piecewise discontinuous polynomials.

\cindex{projection!in {$H(\mathrm{div})$}}%
Let~$\pi_{Q_h}$ denote the~$L_2$ projection
from~$L^2(\Omega)$ into~$Q_h$.
For all~\mbox{$p\in L^2(\Omega)$},
it is defined as 
\mbox{$\pi_{Q_h}(p)=p_h\in Q_h$},
where~$q_h$ is the solution of the following quadratic minimization problem: 
\begin{eqnarray*}
  p_h
  &=&
  \arginf_{q_h\in Q_h}
  \ 
  \int_\Omega
    (p - q_h)^2
    \mathrm{d}x
\end{eqnarray*}
Its solution is characterized as the unique solution
of the following linear system, expressed in variational form:

\ \ $(P_1)$: \emph{find $p_h\in Q_h$ such that}
\begin{eqnarray*}
  \int_\Omega 
    p_h\,q_h
    \,\mathrm{d}x
  &=&
  \int_\Omega 
    p\,q_h
    \,\mathrm{d}x
  ,\ \ \forall q_h\in Q_h
\end{eqnarray*}
Following Roberts and Thomas~\cite[p.~551-552]{CiaLio-vol02-part04-1991},
our aim is to define~$\pi_{V_h}$ as
the~$L_2$ projection from~$H(\mathrm{div},\Omega)$
into $V_h$ and satisfying the following \emph{commuting} property:
\cindex{projection!commuting}%
\begin{center}
    \begin{tikzcd}
        {H(\mathrm{div},\Omega)}
          \arrow[r, "\mathrm{div}"]
          \arrow[d, "\pi_{V_h}"]
        & 
        {L^2(\Omega)}
          \arrow[d, "\pi_{Q_h}"]
        \\
        V_h
          \arrow[r, "\mathrm{div}"]
        & 
        Q_h
    \end{tikzcd}
\end{center}
i.e.
\begin{subequations}
\begin{eqnarray}
  \mathrm{div} (\pi_{V_h}(\boldsymbol{u}))
  &=&
  \pi_{Q_h}(\mathrm{div}\,\boldsymbol{u})
  ,\ \ \forall \boldsymbol{u}\in H(\mathrm{div},\Omega)
  \label{eq-rt-commute}
\end{eqnarray}
In~\cite[p.~553]{CiaLio-vol02-part04-1991}, theorem~6.3,
this projection operator then satisfies an optimal error bound
in the~$H(\mathrm{div},\Omega)$ norm, i.e.:
\begin{eqnarray}
  \| u-\pi_{V_h}(\boldsymbol{u}) \|_{0,2,\Omega}
  +
  \| \mathrm{div} (u-\pi_{V_h}(\boldsymbol{u})) \|_{0,2,\Omega}
  &=&
  \mathcal{O}(h^{k+1}) 
  \label{eq-rt-proj-bound}
\end{eqnarray}
\end{subequations}
Remark that the Lagrange interpolation 
operator~$\pi_h$ from~$H(\mathrm{div},\Omega)$ to~$V_h$
do not necessarily satisfy the commuting property~\eqref{eq-rt-commute}.
Indeed, this depends upon the way the Raviart-Thomas internal
degrees of freedom are chosen and implemented.
When \emph{internal degrees of freedom} are chosen as
integrals over polynomials of degree~$\ell \leq k-1$,
e.g. as in~\cite[p.~551]{CiaLio-vol02-part04-1991}, eqn~(6.8),
then the Lagrange interpolation~$\pi_h$ satisfies
both~\eqref{eq-rt-commute} and~\eqref{eq-rt-proj-bound},
as shown~\cite{CiaLio-vol02-part04-1991}, theorem~6.1 and~6.3.

In practice, it is more efficient to choose
for all the internal degrees of freedom of the Raviart-Thomas
some values of the function on a set of internal points:
this implementation choice has been chosen in~\Rheolef.
In that case, the Lagrange interpolation~$\pi_h$ neither
satisfies the commutation~\eqref{eq-rt-commute}
nor the bound~\eqref{eq-rt-proj-bound}.
More precisely, the Lagrange interpolation error is optimal 
in~$L^2$ norm only while its divergence converges
sub-optimally. 
Thus, with the present choice of the internal degrees of freedom,
there is a need to explicitly compute the projection~$\pi_{V_h}$
that satisfies both~\eqref{eq-rt-commute} and~\eqref{eq-rt-proj-bound}.

For all~\mbox{$\boldsymbol{u}\in H(\mathrm{div},\Omega)$},
this projection is defined as the~$L^2$ projection of~$\boldsymbol{u}$
under the constraint~\eqref{eq-rt-commute}:
\begin{eqnarray*}
  \pi_{V_h}(\boldsymbol{u})
  &=&
  \arginf_{\boldsymbol{v}_h\in V_h}
      \| \boldsymbol{u} - \boldsymbol{v}_h \|^2_{0,2,\Omega}
  \\
  &&
  \mbox{subject to} \ \ 
    \mathrm{div} (\boldsymbol{v}_v)
    \ =\ 
    \pi_{Q_h}(\mathrm{div}\,\boldsymbol{u})
\end{eqnarray*}
Then, $\boldsymbol{u}_h=\pi_{V_h}(\boldsymbol{u})\in V_h$
is equivalently characterized as
the solution of the following saddle-point problem: 
\begin{eqnarray*}
  (\boldsymbol{u}_h,p_h)
  &=&
  \arginf_{\boldsymbol{v}_h\in V_h}
  \ 
  \argsup_{q_h\in Q_h}
  \ 
  L(\boldsymbol{v}_h,q_h)
\end{eqnarray*}
where the Lagrangian~$L$ is defined,
for all $(\boldsymbol{v}_h,q_h)\in V_h\times Q_h$, 
by
\begin{eqnarray*}
  L(\boldsymbol{v}_h,q_h)
  &=&
  \int_\Omega
    \left(
      |
        \boldsymbol{u}
        -
        \boldsymbol{v}_h
      |^2
      +
      q_h
      \mathrm{div}
      (
        \boldsymbol{u}
        -
        \boldsymbol{v}_h
      )
    \right)
    \,\mathrm{d}x
\end{eqnarray*}
The saddle-point of~$L$
is characterized as the unique solution of
a linear system, expressed in variational form.
Moreover, since both~$V_h$ and~$Q_h$ are
spaces of piecewise discontinuous functions,
the linear system writes as a collection
of local linear systems on each element.

\ \ $(P)$: \emph{find $(\boldsymbol{u}_h,p_h)\in V_h\times Q_h$ such that,
   on each element $K\in\mathcal{T}_h$, we have}
\begin{subequations}
\begin{eqnarray}
  \int_K
    (
      \boldsymbol{u}_h.
      \boldsymbol{v}_h
      +
      p_h \mathrm{div}\,\boldsymbol{v}_h
    )
    \mathrm{d}x
  &=&
  \int_K
    \boldsymbol{u}.
    \boldsymbol{v}_h
    \,\mathrm{d}x
  \label{eq-rt-proj-fv-u}
  \\
  \int_K
    q_h \mathrm{div}\,\boldsymbol{u}_h
    \,\mathrm{d}x
  &=&
  \int_K
    q_h \mathrm{div}\,\boldsymbol{u}
    \,\mathrm{d}x
  \label{eq-rt-proj-fv-lambda}
\end{eqnarray}
\end{subequations}
for all $(\boldsymbol{v}_h,q_h)\in V_h\times Q_h$.
\cindex{Lagrange!multiplier}%
Observe that~$q_h$ represents the Lagrange multiplier associated
to the commutation constraint~\eqref{eq-rt-proj-fv-lambda},
which is equivalent to~\eqref{eq-rt-commute}.

Let us introduce the following forms:
\begin{eqnarray*}
  a(\boldsymbol{u}_h,p_h;\ \boldsymbol{v}_h,q_h)
  &=&
  \int_\Omega 
    \left(
      \boldsymbol{u}_h.
      \boldsymbol{v}_h
      +
      p_h \mathrm{div}\,\boldsymbol{v}_h
      +
      q_h \mathrm{div}\,\boldsymbol{u}_h
    \right)
    \mathrm{d}x
  \\
  \ell(\boldsymbol{v}_h,q_h)
  &=&
  \int_\Omega 
    \left(
      \boldsymbol{u}.
      \boldsymbol{v}_h
      +
      q_h \mathrm{div}\,\boldsymbol{u}
    \right)
    \mathrm{d}x
\end{eqnarray*}
The previous problem writes equivalently in abstract form:

\ \ $(P)$: \emph{find $(\boldsymbol{u}_h,p_h)\in V_h\times Q_h$ such that}
\begin{eqnarray*}
  a(\boldsymbol{u}_h,p_h;\ \boldsymbol{v}_h,q_h)
  &=&
  \ell(\boldsymbol{v}_h,q_h)
  ,\ \ \forall (\boldsymbol{v}_h,q_h)\in V_h\times Q_h
\end{eqnarray*}
Note that the matrix associated to the bilinear form~$a$
is symmetric and block-diagonal: it can thus be efficiently inverted
on the fly at the element level during the assembly process.
The following code implement this efficient approach.

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

\myexamplenoinput{commute_rtd_error.cc}

\myexamplenoinput{cosinus_vector.h}
% -------------------------------------

\begin{figure}[htb]
  \begin{center}
  \begin{tabular}{ll}
    \includegraphics[width=0.42\textwidth]{commute-rtd-t-p-err-l2.pdf} &
    \includegraphics[width=0.42\textwidth]{commute-rtd-t-p-err-div-l2.pdf}  \\
    \includegraphics[width=0.42\textwidth]{commute-rtd-t-i-err-l2.pdf} &
    \includegraphics[width=0.42\textwidth]{commute-rtd-t-i-err-div-l2.pdf} 
  \end{tabular}
  \end{center}
  \caption{Raviart-Thomas approximation:
	$\pi_{V_h}$~projection (top) and
	$\pi_h$~interpolation (bottom) errors
	in $L^2$ norm for the approximation and its divergence.
  }
  \label{commute-rtd}
\end{figure}

\subsection*{How to run the program ?}
%-------------------------------------
\begin{verbatim}
  make commute_rtd commute_rtd
  mkgeo_grid -t 10 > square.geo ./commute_rtd square.geo | ./commute_rtd_error 
\end{verbatim}
The program \file{commute_rt.cc} compute both the projection~$\pi_{V_h}(u)$
and the standard Lagrange interpolation~$\pi_h(u)$,
while \file{commute_rtd_error.cc} performs the computation of the 
corresponding errors.
The file \file{cosinus_vector.h} furnishes the function used for the
present test:
\begin{equation*}
  \boldsymbol{u}(x)
  =
  \left( \begin{array}{c}
    \cos(x_0+2x_1) \\
    \sin(x_0-2x_1)
  \end{array} \right)
\end{equation*}
These two last files are not listed here but are available in the \Rheolef\  example directory.
Observe on Fig.~\ref{commute-rtd} that the error for the projection~$\pi_{V_h}(u)$
and its divergence behave as~$\mathcal{O}(h^{k+1})$, which is optimal.
Conversely, the error for the Lagrange interpolation~$\pi_h(u)$
is sub-optimal for the divergence.

In conclusion, the projection~$\pi_{V_h}$ should be used instead of
the interpolation~$\pi_h$
when we want to build an optimal Raviart-Thomas approximation.