% --------------------------------
\subsection{Slope limiters}
% --------------------------------

% Voir~\citet[p.~108]{PieErn-2012} et~\citet{HesWar-2008}, pages~145 et~224.
% On a juste a visiter les voisins par face~:
% pas plus de communications que DG en distribue.
%
% Pas evident de faire un schema implicite BDF(p)~:
% le limiteur est non-lineaire et non-differentiable 
% et le schema ne conserve pas les quantites entre $[0,1]$ 
% Pour que Newton converge bien, il faut etre pseudo-differentiable,
% et meme, le schema pourrait sortir de $[0,1]$ a cause de BDF(p).
%
% En attendant d'avoir la biblio complete, il vaut mieux essayer
% le limiteur classique, en schema explicite.
%
Slope limiters are required when the solution develops
discontinuities: this is a very classical feature of most
solutions of nonlinear hyperbolic equations.
A preliminary version of the slope limiter 
proposed by \citet[p.~208]{Coc-1998}
% See also~\citep{CocHouShu-1990,HotAckMosErhPhi-2004,KriXinRemCheFla-2004,Kri-2007,LuoBauLoh-2007,Kuz-2010} or~\citep[p.~108]{PieErn-2012}
% for the implementation of limiters.
is implemented in \Rheolef: this preliminary version
only supports the $d=1$ dimension and $k=1$ polynomial degree.
Recall that the $k=0$ case do not need any limiter.
More general implementation will support the $d=2,3$ and $k\geq 2$ 
cases in the future.
The details of the limiter implementation is presented in this section:
the impatient reader, who is interested by applications,
can jump to the next section, devoted to the Burgers equation.

\begin{figure}[htb]
     \begin{center}
          \includegraphics[height=6cm]{limiter-fig.pdf}
     \end{center}
     \caption{Limiter: the neighbors elements and the middle edge points.}
   \label{fig-limiter}
\end{figure}
Fig.~\ref{fig-limiter} shows the $d+1$ neighbor elements
$K_i$, $i=0\ldots d$ around an element $d$.
Let $S_i=\partial K\cap \partial K_i$, $i=0\ldots d$ be the $i$-th side of $K$.
We denote by $\boldsymbol{x}_K$, $\boldsymbol{x}_{K_i}$ ad $\boldsymbol{x}_{S_i}$
the barycenters of these elements and sides, $i=0\ldots d$.
When $d=2$, the barycenter $\boldsymbol{x}_{S_i}$ of the edge
belongs to the interior of a triangle 
$(\boldsymbol{x}_{K},\boldsymbol{x}_{K_i},\boldsymbol{x}_{K_{J_{i,1}}})$
for exactly one of the two possible ${J_{i,1}}\neq i$ and $0\leq {J_{i,1}}\leq d$. 
When $d=3$, the barycenter $\boldsymbol{x}_{S_i}$ of the face
belongs to the interior of a tetrahedron
$(\boldsymbol{x}_{K},\boldsymbol{x}_{K_i},\boldsymbol{x}_{K_{J_{i,1}}},\boldsymbol{x}_{K_{J_{i,2}}})$
for exactly one pair $(J_{i,1},J_{i,2})$, up to a permutation,
of the three possible pairs
$J_{i,1},J_{i,2}\neq i$ and $0\leq J_{i,1},J_{i,2}\leq d$.
Let us denote $J_{i,0}=i$.
Then, the vector
$\overrightarrow{\boldsymbol{x}_K \boldsymbol{x}_{S_i}}$
decompose on the basis
$(\overrightarrow{\boldsymbol{x}_K \boldsymbol{x}_{K_{J_{i,k}}}})_{0\leq k\leq d-1}$
as
\begin{eqnarray}
  \overrightarrow{\boldsymbol{x}_K \boldsymbol{x}_{S_i}}
  =
  \sum_{k=0}^{d-1}
  \alpha_{i,k}
  \overrightarrow{\boldsymbol{x}_K \boldsymbol{x}_{K_{J_{i,k}}}}
  \label{eq-hyp-barycenter}
\end{eqnarray}
where $\alpha_{i,k} \geq 0$, $k=0\ldots d-1$.
Let us consider now the patch 
$\omega_K$
composed of $K$
and its $d$ neighbors:
\[
  \omega_K
  =
  K
  \cup
  K_0
  \cup
  \ldots
  \cup
  K_{d}
\]
For any affine function $\xi\in P_1(\omega_K)$ over this patch,
let us denote
\begin{eqnarray*}
  \delta_{K,i}(\xi)
  &=&
  \sum_{k=0}^{d-1}
  \alpha_{i,k}
  \left(
    \xi( \boldsymbol{x}_{K_{J_{i,k}}} )
    -
    \xi( \boldsymbol{x}_K )
  \right)
  , \ \ i=0\ldots d-1
  \\
  &=&
  \xi( \boldsymbol{x}_{S_i} )
  -
  \xi( \boldsymbol{x}_K )
  \ \ \mbox{ from~\eqref{eq-hyp-barycenter}}
\end{eqnarray*}
In other terms, 
  $\delta_{K,i}(\xi)$
represents the departure of the value of $\xi$
at $\boldsymbol{x}_{S_i}$ from its
average $\xi( \boldsymbol{x}_K$ on the element $K$.

Let now $(\varphi_i)_{0\leq i\leq d-1}$ denote the Lagrangian
basis in $K$ associated to the set of nodes 
$(\boldsymbol{x}_{S_i})_{0\leq i\leq d-1}$:
\begin{eqnarray*}
	\varphi_i (\boldsymbol{x}_{S_j})
	= \delta_{i,j}
	, \ \ \ 0\leq i,j \leq d-1
	\\
	\sum_{i=0}^{d-1} \varphi_i (\boldsymbol{x}) = 1
        ,\ \ \forall \boldsymbol{x} \in K
\end{eqnarray*}
The affine function $\xi\in P_1(\omega_K)$ expresses on this basis as
\begin{eqnarray*}
	\xi (\boldsymbol{x})
	&=&
	\xi (\boldsymbol{x}_K)
	+
	\sum_{i=0}^{d-1} \delta_{K,i}(\xi) \, \varphi_i (\boldsymbol{x})
        ,\ \ \forall \boldsymbol{x} \in K
\end{eqnarray*}
Let now $u_h \in \mathbb{P}_{1d}(\mathscr{T}_h)$.
On any element $K\in\mathscr{T}_h$,
let us introduce its average value:
\begin{eqnarray*}
        \bar{u}_{K}
	&=& 
	\Frac{1}{{\rm meas(K)}}
	\int_K
            u_{h} (\boldsymbol{x})
	    \,{\rm d}x
\end{eqnarray*}
and its departure from its average value:
\begin{eqnarray*}
	\tilde{u}_K (\boldsymbol{x})
	&=& 
        u_{h|K} (\boldsymbol{x})
	-
        \bar{u}_{K}
        ,\ \ \forall \boldsymbol{x} \in K
	\\
\end{eqnarray*}
Note that $u_h\not\in P_1(\omega_K)$.
Let us extends $\delta_{K,i}$ to $u_h$ as
\begin{eqnarray*}
  \delta_{K,i}(u_h)
  &=&
  \sum_{k=0}^{d-1}
  \alpha_{i,k}
  \left(
    \bar{u}_{K_{J_{i,k}}}
    -
    \bar{u}_{K}
  \right)
  , \ \ i=0\ldots d-1
\end{eqnarray*}
Since $u_h\not\in P_1(\omega_K)$,
we have
\mbox{$
  \tilde{u}_K (\boldsymbol{x}_{K_{J_{i,k}}})
  \neq
  \delta_{K,i}(u_h)
$}
in general.
The idea is then to capture oscillations by controlling 
the departure of the values
$\tilde{u}_K (\boldsymbol{x}_{K_{J_{i,k}}})$
from the values 
$\delta_{K,i}(u_h)$.
Thus, associate to $u_h \in \mathbb{P}_{1d}(\mathscr{T}_h)$
the quantities
\begin{eqnarray*}
  \Delta_{K,i}(u_h)
  &=&
  {\rm minmod}_{\rm T VB}
  \left(
    \tilde{u}_K (\boldsymbol{x}_{K_{J_{i,k}}})
    ,\ 
    \theta \delta_{K,i}(u_h)
  \right)
\end{eqnarray*}
for all $i=0\ldots d-1$
and where $\theta \geq 1$ is a parameter of the limiter 
and
\begin{eqnarray*}
  {\rm minmod}_{\rm TVB} (a,b)
  &=&
  \left\{ \begin{array}{ll}
	a
	& \mbox{ when }
	  |a| \leq Mh^2
	  \\
        {\rm minmod} (a,b)
	& \mbox{ otherwise }
  \end{array} \right.
\end{eqnarray*}
where $M>0$ is a tunable parameter
which can be evaluated
from the curvature of the initial
datum at its extrema by setting
\begin{eqnarray}
   M
   =
   \sup_{
     \boldsymbol{x} \in\Omega
     , \nabla u_0 (\boldsymbol{x}) = 0
   }
   |\nabla\otimes\nabla u_0|
   \label{eq-dg-limiter-M}
\end{eqnarray}
Introduced by \citet{Shu-1987}, the basic idea is to deactivate the limiter 
when space derivatives are of order $h^2$.
This improves the limiter behavior near smooth local extrema.
The minmod function is defined by
\begin{eqnarray*}
  {\rm minmod} (a,b)
  &=&
  \left\{ \begin{array}{ll}
    	{\rm sgn}(a)\, {\rm min}(|a|,|b|)
	& \mbox{ when }
    	  {\rm sgn}(a) = {\rm sgn}(b)
	\\
	0
	& \mbox{ otherwise }
  \end{array} \right.
\end{eqnarray*}
Then, for all $i=0\ldots d-1$ we define
\begin{eqnarray*}
  r_K(u_h)
  &=&
  \Frac{ 
    {\displaystyle
  	\sum_{j=0}^{d-1}
	\max (0, -\Delta_{K,j}(u_h))
    }
  }{
    {\displaystyle
  	\sum_{j=0}^{d-1}
	\max (0, \Delta_{K,j}(u_h))
    }
  }
  \ \geq 0
  \\
  \hat{\Delta}_{K,i}(u_h)
  &=&
  \phantom{+}
    \min(1,r_K(u_h))
  \,\max (0, \Delta_{K,i}(u_h))
  \\
  && 
  -
     \min(1,1/r_K(u_h))
  \, \max (0, -\Delta_{K,i}(u_h))
  ,\ \ i=0\ldots d-1
    \mbox{ when } r_K(u_h) \neq 0
\end{eqnarray*}
Finally, the limited function $\Lambda_h(u_h)$ is defined 
element by element for all element $K\in\mathscr{T}_h$
for all $\boldsymbol{x} \in K$ by
\begin{eqnarray*}
  \Lambda_h(u_h)_{|K}
  (\boldsymbol{x})
  &=&
  \left\{ \begin{array}{ll}
    {\displaystyle
      \bar{u}_K
      + 
      \sum_{i=1}^{d-1}
        \Delta_{K,i}(u_h)
        \ \varphi_i (\boldsymbol{x})
    }
    & \mbox{ when } r_K(u_h) = 0
    \\
    {\displaystyle
      \bar{u}_K
      + 
      \sum_{i=1}^{d-1}
        \hat{\Delta}_{K,i}(u_h)
        \ \varphi_i (\boldsymbol{x})
    }
    & \mbox{ otherwise }
  \end{array} \right.
\end{eqnarray*}
Note that there are two types
of computations involved in the limiter:
one part is independent of $u_h$ and depends
only upon the mesh:
 $J_{i,k}$ and $\alpha_{i,k}$ on each element.
It can be computed one time for all.
The other part depends upon the values of $u_h$.

Note that the limiter preserves the average value
of $u_h$ on each element $K$ and also the functions that
are globally affine on the patch $\omega_K$.
Also we have, inside each element $K$
and for all side index $i=0\ldots d-1$:
\begin{eqnarray*}
  \left|
    \Lambda_h(u_h)_{|K} (\boldsymbol{x}_{S_i})
    -
    \bar{u}_{K}
  \right|
  \ \leq \ 
  \max
  \left(
     |\Delta_{K,i}(u_h)|
     ,\,
     |\hat{\Delta}_{K,i}(u_h)|
  \right)
  \ \leq \ 
  |\Delta_{K,i}(u_h)|
  \ \leq \ 
  \left|
    u_{h|K} (\boldsymbol{x}_{S_i})
    -
    \bar{u}_{K}
  \right|
\end{eqnarray*}
It means that, inside each element, the gradient of
the $P_1$ limited function is no larger than that
of the original one.

The limiter on an element close to the boundary 
should takes into account the inflow condition,
see \citet{CocHouShu-1990}.