\subsection{The Stokes problem}
\label{sec-stokes}
\cindex{problem!Stokes}%
\cindex{boundary condition!Dirichlet}%
\cindex{benchmark!driven cavity flow}%

\subsubsection*{Formulation}
  Let us consider the Stokes problem for the driven cavity
  in $\Omega = ]0,1[^d$, $d = 2,3$.
  The problem writes:

  \ \ \ $(S)$ {\it find ${\bf u}=(u_0,\ldots,u_{d-1})$ and $p$ defined in $\Omega$ such that:}
  \[
    \begin{array}{ccccl}
      -\ {\bf div}(2D({\bf u})) &+& \bnabla p &=& 0 \ {\rm in}\ \Omega, \\
      -\ {\rm div}\,{\bf u}     &&            &=& 0 \ {\rm in}\ \Omega, \\
      {\bf u} && &=& (1,0) \ {\rm on}\ \Gamma_{\rm top}, \\
      {\bf u} && &=& 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm right} \cup \Gamma_{\rm bottom}, \\
      \Frac{\partial u_0}{\partial {\bf n}} &=&
      \Frac{\partial u_1}{\partial {\bf n}} &=&
      u_2 = 0 \ {\rm on}\ \Gamma_{\rm back} \cup \Gamma_{\rm front} \ \mbox{ when } d=3,
    \end{array}
  \]
  where $D({\bf u}) = (\bnabla {\bf u} + \bnabla {\bf u}^T)/2$.
  The boundaries are represented on 
  Fig.~\ref{fig-square-cube}, page~\pageref{fig-square-cube}.

  The variational formulation of this problem expresses:
 
  \ \ \ $(VFS)$ {\it find ${\bf u}\in {\bf V}(1)$ and $p \in L^2_0(\Omega)$ such that:}
  \[
     \begin{array}{lcccl}
      a({\bf u},{\bf v}) &+& b({\bf v}, p) &=& 0, \ \forall {\bf v}\in {\bf V}(0), \\
      b({\bf u},q) && &=& 0, \ \forall q \in L^2_0(\Omega),
     \end{array}
  \]
  where
  \begin{eqnarray*}
       a({\bf u},{\bf v}) &=& \int_\Omega 2 D({\bf u}) : D({\bf v}) \, dx, \\
       b({\bf v}, q) &=& -\int_\Omega {\rm div}({\bf v}) \, q \, dx.  \\
	{\bf V}(\alpha) &=& \{ {\bf v} \in (H^1(\Omega))^2; \  
      		{\bf v} = 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm right} \cup \Gamma_{\rm bottom}
		\mbox{ and }
      		{\bf v} = (\alpha,0) \ {\rm on}\ \Gamma_{\rm top} \}, 
		\ \mbox{ when } d=2,
	    \\
	{\bf V}(\alpha) &=& \{ {\bf v} \in (H^1(\Omega))^3; \  
      		{\bf v} = 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm right} \cup \Gamma_{\rm bottom},
            \\
          && \hspace{2cm}
      		{\bf v} = (\alpha,0,0) \ {\rm on}\ \Gamma_{\rm top} 
		\mbox{ and }
      		v_2 = 0 \ {\rm on}\ \Gamma_{\rm back} \cup \Gamma_{\rm front} \}, 
		\ \mbox{ when } d=3,
	    \\
       L^2_0(\Omega) &=& \{ q \in L^2(\Omega); \ \int_\Omega q \, dx = 0 \}.
  \end{eqnarray*}
\subsubsection*{Approximation}
\apindex{mixed}%
\apindex{P2-P1, Taylor-Hood}%
  The \citet*{hood-taylor-73} finite element approximation
  of the Stokes problem is considered.
  We introduce a mesh ${\cal T}_h$ of $\Omega$
  and the following finite dimensional spaces:
  \begin{eqnarray*}
      {\bf X}_h &=& \{ {\bf v} \in (H^1(\Omega))^d; \
          {\bf v}_{/K} \in (P_2)^d, \
          \forall K \in {\cal T}_h \}, \\
      {\bf V}_h(\alpha) &=& {\bf X}_h \cap {\bf V}(\alpha), \\
      Q_h &=& \{ q \in L^2(\Omega))\cap C^0(\bar{\Omega}); \
          q_{/K} \in P_1, \
          \forall K \in {\cal T}_h \},
  \end{eqnarray*}
\cindex{form!{$2D({\bf u}):D({\bf v})$}}%
\cindex{form!{${\rm div}({\bf u})\,q$}}%
  The approximate problem writes:

  \ \ \ $(VFS)_h$ {\it find ${\bf u}_h \in {\bf V}_h(1)$ and $p \in Q_h$ such that:}
  \begin{equation}
    \begin{array}{lcccl}
      a({\bf u}_h,{\bf v}) &+& b({\bf v}, p_h) &=& 0, \ \forall {\bf v}\in {\bf V}_h(0), \\
      b({\bf u}_h,q) && &=& 0, \ \forall q \in Q_h.
    \end{array}
    \label{eq-fvh-stokes}
  \end{equation}

% ----------------------------
\myexamplelicense{cavity.h}

\myexamplelicense{stokes_cavity.cc}
% ----------------------------
\subsubsection*{Comments}
  The spaces and boundary conditions and grouped
  in specific functions, defined in file \reffile{cavity.h}.
  This file is suitable for a future re-usage.
  Next, forms are defined as usual, in file \reffile{stokes_cavity.cc}.

\cindex{method!conjugate gradient algorithm}%
\cindex{preconditioner}%
\cindex{preconditioner!for Stokes problem}%
  The problem admits the following matrix form:
  \[
     \left( \begin{array}{cc} 
	{\tt a} & {\tt trans(b)} \\ 
        {\tt b} & 0 
     \end{array} \right)
     \left( \begin{array}{c} 
        {\tt uh} \\ 
        {\tt ph} 
     \end{array} \right)
	=
     \left( \begin{array}{c} 
	0 \\
        0
     \end{array} \right)
  \]
  An initial value for the pressure field is provided:
\begin{lstlisting}[numbers=none,frame=none]
  field ph (Qh, 0);
\end{lstlisting}
\label{ref-pcg-abtb}%
\clindex{problem_mixed}%
The solve call to the Stokes problem writes:
\begin{lstlisting}[numbers=none,frame=none]
  problem_mixed stokes (a, b);
  stokes.solve (field(Xh,0), field(Qh,0), uh, ph);
\end{lstlisting}
The two first arguments of the \code{solve} member function
represents the zero right-hand-side of the problem.
For two-dimensional geometries ($d=2$), this system is solved by 
a direct method (see \citealp[p.~41]{Sar-2016-cfma}).
Conversely, for tridimensional geometries ($d=3$), it is solved
by the preconditioned conjugate gradient algorithm
(see \citealp[p.~56]{Sar-2016-cfma}).
In that case, the preconditioner is by default
the mass matrix \code{mp} for the pressure space:
as showed by \citet{Kla-1998}, the number of iterations need by 
the conjugate gradient algorithm to reach a given precision
is then independent of the mesh size.
For more details, see the Rheolef reference manual related to mixed solvers,
available on the web site and via the unix command:
\clindex{reference manual}%
\clindex{man}%
\begin{verbatim}
  man problem_mixed
\end{verbatim}
%% \clindex{solver}%
%%   When $d=2$, it is interesting to turn to direct methods and factorize
%%   the whole matrix of the linear system. As the pressure is defined up to a constant,
%% \cindex{Lagrange!multiplier}%
%%   the whole matrix is singular. By adding a Lagrange multiplier that impose 
%%   a null average pressure value, the system becomes regular and the modified
%%   matrix can be inverted. Such a technique has already been presented 
%%   in section~\ref{sec-neumann-laplace} for the Neumann-Laplace problem.
%%   Finally, he choice between iterative and direct algorithm for the Stokes solver is
%%   automatically done, regarding the geometry dimension.

\subsubsection*{How to run the program}
  \begin{figure}[htb]
     %\begin{center}
       \mbox{} \hspace{-0.2cm}
       \begin{tabular}{ccc}
	  \includegraphics[height=7cm]{stokes-cavity-2d-fig.png} &
          \includegraphics[height=7cm]{stokes-cavity-3d-fig.png} &
	  \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{The velocity visualization for $d=2$ and $d=3$ with stereo anaglyph.}
     \label{fig-stokes-cavity-velocity}
  \end{figure}
  We assume that the previous code is contained in
  the file \reffile{stokes_cavity.cc}.
  Then, compile the program as usual (see page~\pageref{makefile}):
\begin{verbatim}
  make stokes_cavity
\end{verbatim}
  and enter the commands:
\begin{verbatim}
  mkgeo_grid -t 10 > square.geo
  ./stokes_cavity square > square.field
\end{verbatim}
  The previous command solves the problem for the corresponding mesh
  and writes the solution in a \filesuffix{.field} file.
  Run the velocity vector visualization :
\pindexopt{field}{-velocity}
\begin{verbatim}
  field square.field -velocity
\end{verbatim}
  Run also some scalar visualizations:
\pindexopt{field}{-mark}
\begin{verbatim}
  field square.field -comp 0
  field square.field -comp 1
  field square.field -mark p
\end{verbatim}
\findex{catchmark}%
  Note the \code{-mark} option to the \code{field} command:
  the file reader jumps to the label and then starts to read the selected field.
  Next, perform another computation on a finer mesh:
\begin{verbatim}
  mkgeo_grid -t 20 > square-20.geo
  ./stokes_cavity square-20.geo > square-20.field
\end{verbatim}
  and observe the convergence.

  Finally, let us consider the three dimensional case:
\begin{verbatim}
  mkgeo_grid -T 5 > cube.geo
  ./stokes_cavity cube.geo > cube.field
\end{verbatim}
  and the corresponding visualization:
\begin{verbatim}
  field cube.field -velocity -scale 3
  field cube.field -comp 0
  field cube.field -comp 1
  field cube.field -comp 2
  field cube.field -mark p
\end{verbatim}
 
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Computing the vorticity}
%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Formulation and approximation}
\cindex{vorticity}%
\cindex{operator!curl}%

\cindex{form!{${\rm curl}({\bf u}).\boldsymbol{\xi}$}}%
When $d=2$, we define~\citep[page~30]{GirRav-1986} for
any distributions $\phi$ and ${\bf v}$:
\begin{eqnarray*}
    {\bf curl}\,\phi &=& \left(
			    \Frac{\partial \phi}{\partial x_1},\,
			   -\Frac{\partial \phi}{\partial x_0} 
			 \right), \\
    {\rm curl}\,{\bf v} &=& \Frac{\partial v_1}{\partial x_0} - \Frac{\partial v_0}{\partial x_1},
\end{eqnarray*}
and when $d=3$:
\[
    {\bf curl}\,{\bf v} = \left(
			    \Frac{\partial v_2}{\partial x_1} -\Frac{\partial v_1}{\partial x_2}, \,
			    \Frac{\partial v_0}{\partial x_2} -\Frac{\partial v_2}{\partial x_0}, \,
			    \Frac{\partial v_1}{\partial x_0} -\Frac{\partial v_0}{\partial x_1}
			    \right)
\]
Let ${\bf u}$ be the solution of the Stokes problem~$(S)$.
The vorticity is defined by:
\[
    \begin{array}{cccc}
	\omega  &=& {\rm curl} \, {\bf u} & \mbox{ when } d=2, \\
	\bomega &=& {\bf curl} \, {\bf u} & \mbox{ when } d=3.
    \end{array}
\]
Since the approximation of the velocity is piecewise quadratic,
we are looking for a discontinuous piecewise linear vorticity
field that belongs to:
\[
    \begin{array}{cccc}
      Y_h &=& \{ \xi \in L^2(\Omega); \, \xi_{/K} \in P_1, \ \forall K \in {\cal T}_h \},
	  & \mbox{ when } d=2 \\
      {\bf Y}_h &=& \{ \bxi \in (L^2(\Omega))^3; \, \xi_{i/K} \in P_1, \ \forall K \in {\cal T}_h \},
	  & \mbox{ when } d=3
    \end{array}
\]
The approximate variational formulation writes:
\[
    \begin{array}{cccc}
	\omega_h \in Y_h, \ \ 
	  \int_\Omega \omega_h \, \xi \, dx 
	  &=& \int_\Omega {\rm curl} \, {\bf u}_h \, \xi \, dx, \ \forall \xi \in Y_h
	  & \mbox{ when } d=2, \\
	&&& \\
	\bomega \in {\bf Y}_h, \ \ 
	  \int_\Omega \bomega_h . \bxi \, dx 
	  &=& \int_\Omega {\bf curl} \, {\bf u}_h . \bxi \, dx, \ \forall \bxi \in {\bf Y}_h
          & \mbox{ when } d=3.
    \end{array}
\]
% -------------------------
\myexamplelicense{vorticity.cc}
% -------------------------

\subsubsection*{Comments}
\apindex{discontinuous}%
  As for the stress tensor (see \code{stress.cc}, page~\pageref{stress.cc}),
  the vorticity is obtained by a direct interpolation of the $u_h$ first derivatives
  and its approximation is \emph{discontinuous} at inter-element boundaries.

\subsubsection*{How to run the program}
  \begin{figure}[htb]
     %\begin{center}
       \mbox{}\hspace{-0.5cm}
       \begin{tabular}{ccc}
	  \includegraphics[height=7cm]{vorticity-cavity-2d-fig.png} &
	  \includegraphics[height=7cm]{vorticity-cavity-3d-fig.png} &
	  \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{The vorticity: elevation view for $d=2$ and 
     	vector representation for $d=3$ (with anaglyph).}
     \label{fig-stokes-cavity-vorticity}
  \end{figure}
For $d=2$, just enter:
\begin{verbatim}
  make vorticity
  ./vorticity < square.field | field -elevation -stereo -
\end{verbatim}
and you observe a discontinuous piecewise linear representation
of the approximate vorticity.
\cindex{singular solution}%
Also, the vorticity presents two sharp peaks at the upper corners
of the driven cavity: the vorticity is unbounded and the peaks will
increase with mesh refinements. This singularity of the solution
is due to the boundary condition for the first component of the
velocity $u_0$ that jumps from zero to one at the corners.
\cindex{projection!in {$L^2$} norm}%
The approximate vorticity field can also be projected on
a continuous piecewise linear space, using the \code{-proj}
option (See Fig.~\ref{fig-stokes-cavity-vorticity} left):
\begin{verbatim}
  ./vorticity < square.field | field -elevation -stereo -nofill -
  ./vorticity < square.field | field -elevation -stereo -proj -
\end{verbatim}
For $d=3$, the whole vorticity vector can also be visualized
(See Fig.~\ref{fig-stokes-cavity-vorticity} right):
\begin{verbatim}
  ./vorticity < cube.field | field -proj -velocity -stereo -
\end{verbatim}
In the previous command, the \code{-proj} option has been
used: since the 3D render has no support for discontinuous 
piecewise linear fields, the P1-discontinuous field is transformed into a
P1-continuous one, thanks to a $L^2$ projection.
P1
The following command shows the second component of the vorticity vector,
roughly similar to the bidimensional case.
\begin{verbatim}
  ./vorticity < cube.field | field -comp 1 - 
  ./vorticity < cube.field | field -comp 1 -proj -
\end{verbatim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Computing the stream function}
\label{sec-streamf}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection*{Formulation and approximation}
\cindex{stream function}
When $d=3$, the stream function is a vector-valued field
$\bpsi$ that satisfies~\citep[page~90]{GirRav-1986}:
${\bf curl}\,\bpsi = {\bf u}$ and ${\rm div}\,\bpsi=0$.
From the identity:
\[
   {\bf curl}\,{\bf curl}\,\bpsi = -\Delta\bpsi + \nabla({\rm div}\,\bpsi)
\]
we obtain the following characterization of $\bpsi$~:
\[
    \begin{array}{rccl}
	-\Delta \,\bpsi &=& {\bf curl} \, {\bf u} & \mbox{ in } \Omega, \\
	\bpsi              &=& 0                   & \mbox{ on } 
      		\Gamma_{\rm back} \cup \Gamma_{\rm front} \cup \Gamma_{\rm top} \cup \Gamma_{\rm bottom}, \\
	\Frac{\partial \bpsi}{\partial n} &=& 0    &  \mbox{ on } 
      		\Gamma_{\rm left} \cup \Gamma_{\rm right}.
    \end{array}
\]
When $d=2$, the stream function $\psi$ is a scalar-valued field the
solution of the following problem~\citep[page~88]{GirRav-1986}:
\[
    \begin{array}{rccl}
	-\Delta \,\psi   &=& {\rm curl} \, {\bf u} &\mbox{ in } \Omega, \\
	\psi             &=& 0                      & \mbox{ on } \partial\Omega.
    \end{array}
\]
% --------------------------------
\myexamplelicense{streamf_cavity.cc}
% --------------------------------

\subsubsection*{How to run the program}

  \begin{figure}[htb]
     %\begin{center}
       \mbox{}\hspace{-0.3cm}
       \begin{tabular}{ccc}
	  \includegraphics[width=7cm]{streamf-cavity-2d-fig.pdf} &
          \includegraphics[width=8.5cm]{streamf-cavity-3d-fig.png} &
	  \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{The stream function visualization:
     	isolines for $d=2$, and combined vectors and isonorm surface for $d=3$.}
     \label{fig-stokes-cavity-streamf}
  \end{figure}
  For $d=2$, just enter (see Fig.~\ref{fig-stokes-cavity-streamf} left):
\begin{verbatim}
  make streamf_cavity
  ./streamf_cavity  < square.field | field -bw -
\end{verbatim}
  For $d=3$, the whole stream function vector can be visualized:
\begin{verbatim}
  ./streamf_cavity < cube.field | field -velocity -scale 10 -
\end{verbatim}
  The second component of the stream function is showed by:
\begin{verbatim}
  ./streamf_cavity < cube.field | field -comp 1 -
\end{verbatim}
The combined representation of Fig.~\ref{fig-stokes-cavity-streamf}.right
has been obtained in two steps.
First, enter:
\pindexopt{field}{-noclean}%
\pindexopt{field}{-noexecute}%
\begin{verbatim}
  ./streamf_cavity < cube.field | field -comp 1 -noclean -noexecute -
  mv output.vtk psi1.vtk
  ./streamf_cavity < cube.field | field -velocity -
\end{verbatim}
\pindex{paraview}
\fiindex{\filesuffix{.vtk} vtk file}
The \code{-noclean -noexecute} options cause the creation of the \filesuffix{.vtk} file
for the second component, without running the \code{paraview} visualization.
Next, in the \code{paraview} window associated to the whole stream function, 
select the \code{File->Open} menu and load \file{psi1.vtk} and click on the
green button \code{Apply}.
Finally, select the \code{Filters/Common/Contours} menu:
the isosurface appears.
Observe that the 3D stream function is mainly represented by its second component.