\subsection{The linear elasticity problem}
\label{sec-elasticity}
\pbindex{elasticity}%
\cindex{boundary condition!Dirichlet}%
\cindex{boundary condition!Neumann}%
\cindex{benchmark!embankment}%

\subsubsection*{Formulation}

\cindex{tensor!Cauchy stress}%
  The total Cauchy stress tensor expresses:
  \begin{equation}
    \sigma({\bf u}) = \lambda \, {\rm div}({\bf u}).I + 2 \mu D({\bf u})
    \label{eq-elasticity-cauchy}
  \end{equation}
\cindex{Lam\'e coefficients}
  where $\lambda$ and $\mu$ are the Lam\'e coefficients.
\cindex{operator!gradient}%
\cindex{operator!divergence}%
\cindex{operator!gradient!symmetric part}%
  Here, $D({\bf u})$ denotes the symmetric part of the gradient operator
  and ${\rm div}$ is the divergence operator.
  Let us consider the elasticity problem for the 
  {\em embankment}, in $\Omega = ]0,1[^d$, $d = 2,3$.
  The problem writes:

  \ \ \ $(P)$: {\it find ${\bf u}=(u_0,\ldots,u_{d-1})$, defined in $\Omega$, such that:}
  \begin{equation}
    \begin{array}{rcl}
      -\ {\bf div}\,\sigma({\bf u}) &=& {\bf f} \ {\rm in}\ \Omega, \\
      \Frac{\partial {\bf u}}{\partial {\bf n}} &=& 0 \ {\rm on}\ \Gamma_{\rm top} \cup \Gamma_{\rm right} \\
      {\bf u} &=& 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm bottom}, \\
      {\bf u} &=& 0 \ {\rm on}\ \Gamma_{\rm front} \cup \Gamma_{\rm back}, \mbox{ when } d=3
    \end{array}
    \label{eq-elasticity}
  \end{equation}
  where ${\bf f} = (0,-1)$ when $d=2$ and ${\bf f}=(0,0,-1)$ when $d=3$.
  The Lam\'e coefficients are assumed to satisfy $\mu>0$ and
  $\lambda+\mu>0$.
  Since the problem is linear, we can suppose that $\mu=1$ without
  any loss of generality.
 \begin{figure}[htb]
    \centerline{\includegraphics{square-cube-fig.pdf}}
    \caption{The boundary domains for the square and the cube.}
    \label{fig-square-cube}
 \end{figure}
  recall that, in order to avoid mistakes with the \code{C++} style indexes,
  we denote by $(x_0,\ldots,x_{d-1})$ the Cartesian coordinate
  system in $\mathbb{R}^d$.

  For $d=2$ we define the boundaries:
  \[
     \begin{array}{llllll}
	\Gamma_{\rm left}   &=& \{0\} \times ]0,1[, & \Gamma_{\rm right} &=& \{1\}\times ]0,1[  \\
	\Gamma_{\rm bottom} &=& ]0,1[ \times \{0\}, & \Gamma_{\rm top}   &=& ]0,1[ \times \{1\} \\
     \end{array}
  \]
  and for $d=3$:
  \[
     \begin{array}{llllll}
	\Gamma_{\rm back} &=& \{0\} \times ]0,1[^2,       & \Gamma_{\rm front}   &=& \{1\} \times ]0,1[^2 \\
	\Gamma_{\rm left}   &=& ]0,1[ \times \{0\} \times ]0,1[,     & \Gamma_{\rm right} &=& ]0,1[ \times \{1\}\times ]0,1[ \\
	\Gamma_{\rm bottom} &=& ]0,1[^2 \times \{0\}, & \Gamma_{\rm top}   &=& ]0,1[^2 \times \{1\}
     \end{array}
  \]
  These boundaries are represented on Fig.~\ref{fig-square-cube}.

  The variational formulation of this problem expresses:
 
  \ \ \ $(VF)$: {\it find ${\bf u}\in {\bf V}$ such that:}
  \begin{equation}
      a({\bf u},{\bf v}) = l({\bf v}), \ \forall {\bf v}\in {\bf V},
    \label{eq-elasticity-fv}
  \end{equation}
  where
  \begin{eqnarray*}
       a({\bf u},{\bf v}) &=& \int_\Omega \left( \lambda {\rm div}\,{\bf u} \,{\rm div}\,{\bf v} 
				+ 2D({\bf u}) : D({\bf v}) \right) \, dx, \\
       l({\bf v}) &=& \int_\Omega {\bf f} . {\bf v} \, dx, \\
       {\bf V} &=& \{ {\bf v} \in (H^1(\Omega))^2; \  
      		{\bf v} = 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm bottom} \}, \mbox{ when } d=2 \\
       {\bf V} &=& \{ {\bf v} \in (H^1(\Omega))^3; \  
      		{\bf v} = 0 \ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm bottom}
      		\cup \Gamma_{\rm right} \cup \Gamma_{\rm back} \}, \mbox{ when } d=3
  \end{eqnarray*}

\subsubsection*{Approximation}
  We introduce a mesh ${\cal T}_h$ of $\Omega$
  and for $k\geq 1$, the following finite dimensional spaces:
  \begin{eqnarray*}
      {\bf X}_h &=& \{ {\bf v}_h \in (H^1(\Omega))^d; \
          {\bf v}_{h/K} \in (P_k)^d, \
          \forall K \in {\cal T}_h \}, \\
      {\bf V}_h &=& {\bf X}_h \cap {\bf V}
  \end{eqnarray*}
  The approximate problem writes:

  \ \ \ $(VF)_h$: {\it find ${\bf u}_h \in {\bf V}_h$ such that:}
  \[
      a({\bf u}_h,{\bf v}_h) = l({\bf v}_h), \ \forall {\bf v}_h\in {\bf V}_h
  \]

%--------------------------
\myexamplelicense{embankment.cc}
\myexamplelicense{embankment.icc}
%--------------------------

\subsubsection*{Comments}
\cindex{form!{$2D({\bf u}):D({\bf v})+\lambda{\rm div}\,{\bf u}\,{\rm div}\,{\bf v}$}}%
\apindex{P1}%
\apindex{P2}%
The space is defined in a separate file \reffile{embankment.icc},
since it will be reused in others examples along this chapter:
\begin{lstlisting}[numbers=none,frame=none]
  space Vh (omega, "P2", "vector");
\end{lstlisting}
Note here the multi-component specification \code{"vector"} as a supplementary
argument to the \code{space} constructor.
The boundary condition contain a special cases for bi- and tridimensional cases.
The right-hand-side ${\bf f}_h$ represents the dimensionless gravity forces,
oriented on the vertical axis: the last component of ${\bf f}_h$ is set
to $-1$ as:
\begin{lstlisting}[numbers=none,frame=none]
    fh [d-1] = -1;
\end{lstlisting}
The code for the bilinear form $a(.,.)$ and the linear one $l(.)$ are closed
to their mathematical definitions: 
\begin{lstlisting}[numbers=none,frame=none]
  form  a  = integrate (lambda*div(u)*div(v) + 2*ddot(D(u),D(v)));
  field lh = integrate (dot(f,v));
\end{lstlisting}

\findex{catchmark}%
Finally, the $1/\lambda$ parameter and the
multi-field result are printed, using
mark labels, thanks to the \code{catchmark} stream manipulator.
Labels are convenient for post-processing purpose,
as we will see in the next paragraph.

\subsubsection*{How to run the program}
\label{sec-howtorun-elasticity}
  \begin{figure}[htb]
     \mbox{}\hspace{-2cm}
     %\begin{center}
       \begin{tabular}{cccc}
	    \includegraphics[height=7cm]{embankment-2d-fig.pdf} &
	    \hspace{0.2cm} &
	    \includegraphics[height=7cm]{embankment-2d-fill-fig.pdf} &
	    \\
	    \includegraphics[width=8cm]{embankment-3d-fig.png} &
	    \hspace{0.2cm} &
	    \includegraphics[width=8.5cm]{embankment-3d-fill-fig.png} &
	    \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{The linear elasticity for $\lambda=1$ and $d=2$ and $d=3$:
	both wireframe and filled surfaces ; stereoscopic anaglyph mode for 3D solutions.
     }
     \label{fig-embankment-deformation}
  \end{figure}
Compile the program as usual (see page~\pageref{makefile}):
\begin{verbatim}
  make embankment
\end{verbatim}
and enter the commands:
\begin{verbatim}
  mkgeo_grid -t 10 > square.geo
  geo square.geo
\end{verbatim}
The triangular mesh has four boundary domains, named \code{left},
\code{right}, \code{top} and \code{bottom}.
Then, enter:
\begin{verbatim}
  ./embankment square.geo P1 > square-P1.field
\end{verbatim}
\fiindex{\filesuffix{.field} multi-component field}%
The previous command solves the problem for the corresponding mesh
and writes the multi-component solution in the \filesuffix{.field} file format.
\pindexoptopt{visualization}{mesh}{deformed}%
Run the deformation vector field visualization:
\pindexopt{field}{-nofill}%
\begin{verbatim}
  field square-P1.field
  field square-P1.field -nofill
\end{verbatim}
It bases on the default \code{paraview} render for 2D and 3D geometries.
The view is shown on Fig.~\ref{fig-embankment-deformation}. 
\pindex{gnuplot}%
If you are in trouble with \code{paraview}, you can switch to
the simpler \code{gnuplot} render:
\begin{verbatim}
  field square-P1.field -nofill -gnuplot
\end{verbatim}
The unix manual for the \code{field} command is available as:
\clindex{reference manual}%
\clindex{man}%
\begin{verbatim}
  man field
\end{verbatim}
A specific field component can be also selected for a scalar visualization:
\pindexopt{field}{-comp}%
\begin{verbatim}
  field -comp 0 square-P1.field
  field -comp 1 square-P1.field
\end{verbatim}
Next, perform a $P_2$ approximation of the solution:
\begin{verbatim}
  ./embankment square.geo P2 > square-P2.field
  field square-P2.field -nofill
\end{verbatim}
Finally, let us consider the three dimensional case
\pindexopt{field}{-fill}%
\pindexopt{field}{-stereo}%
\begin{verbatim}
  mkgeo_grid -T 10 > cube.geo
  ./embankment cube.geo P1 > cube-P1.field
  field cube-P1.field -stereo
  field cube-P1.field -stereo -fill
\end{verbatim}
\cindex{visualization!stereoscopic anaglyph}%
The two last commands show the solution in 3D stereoscopic anaglyph mode.
The graphic is represented on Fig.~\ref{fig-embankment-deformation}.
The $P_2$ approximation writes:
\begin{verbatim}
  ./embankment cube.geo P2 > cube-P2.field
  field cube-P2.field
\end{verbatim}

% ===================================
\subsection{Computing the stress tensor}
% ===================================

\subsubsection*{Formulation and approximation}
  The following code computes the total Cauchy stress tensor,
  reading the Lam\'e coefficient $\lambda$ and the deformation
  field ${\bf u}_h$ from a \filesuffix{.field} file.
  Let us introduce:
  \[
     T_h = \{ \tau_h \in (L^2(\Omega))^{d\times d}; \ 
		\tau_h = \tau_h^T \mbox{ and }
		\tau_{h;ij/K} \in P_{k-1}, \ \forall K\in {\cal T}_h,\ 
		1\leq i,j \leq d \}
  \]
  This computation expresses:

  \ \ \ {\it find $\sigma_h$ such that:}
  \[
	m(\sigma_h, \tau) = b(\tau, {\bf u}_h), \forall \tau \in T_h
  \]
  where
  \begin{eqnarray*}
	m(\sigma, \tau) &=& \int_\Omega \sigma : \tau \, dx
	  ,
	  \\
        b(\tau, {\bf u})
	&=&
	\int_\Omega 
	  \left(
		2 D({\bf u}) : \tau \, dx
        	+ 
		\lambda {\rm div}({\bf u})\,{\rm tr}(\tau)
	  \right)
	  \, dx
	  ,
  \end{eqnarray*}
where ${\rm tr}(\tau) = \sum_{i=1}^d \tau_{ii}$ is the trace
of the tensor $\tau$.

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

\subsubsection*{Comments}
  In order to our code \code{stress.cc} to apply also for the forthcoming incompressible case
  $\lambda=+\infty$, the Lam\'e coefficient is introduced as $1/\lambda$.
  Its value is zero in the incompressible case.
  By this way, the previous code applies for any deformation
  field, and is not restricted to our {\em embankment} problem.
\apindex{discontinuous}%
  The stress tensor is obtained by a direct interpolation of the $u_h$ first derivatives.
  As $u_h$ is continuous and piecewise polynomial $P_k$,
  its derivatives are also piecewise polynomials with degree $k-1$,
  but \emph{discontinuous} at inter-element boundaries~: this approximation is denoted as $P_{k-1,d}$.
  Thus, the stress tensor belongs to the space $T_h$ with the $P_{k-1,d}$ element.

\subsubsection*{How to run the program}
\cindex{tensor!visualization as ellipsoid}
  \begin{figure}[htb]
     %\begin{center}
       \begin{tabular}{cccc}
	  \includegraphics[scale=0.30]{embankment-2d-tensor-fig.png} &
	  \hspace{0.1cm} &
	  \includegraphics[scale=0.30]{embankment-3d-tensor-fig.png} &
	  \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{The stress tensor visualization (linear elasticity $\lambda=1$).}
     \label{fig-embankment-tensor}
  \end{figure}
  
  First, compile the program:
\begin{verbatim}
  make stress
\end{verbatim}
\cindex{tensor!field}%
\pindexopt{field}{-proj}%
The visualization for the stress tensor as ellipses writes:
\pindex{paraview}%
\begin{verbatim}
  ./stress < square-P1.field > square-stress-P1.field
  ./stress < square-P2.field > square-stress-P2.field
  field square-stress-P1.field -proj
  field square-stress-P2.field -proj
\end{verbatim}
This visualization based on \code{paraview} requires
the \code{TensorGlyph} feature, available since \code{paraview-5.0}.
\apindex{discontinuous}%
\cindex{projection!in {$L^2$} norm}%
\apindex{P0}%
\apindex{P1d}%
\apindex{P1}%
\apindex{P2}%
\pindexopt{field}{-stereo}%
  Recall that the stress, as a derivative of the deformation, is 
  P0 (resp. P1d) and discontinuous when the deformation is P1 (resp. P2) and continuous.
  The approximate stress tensor field is projected as
  a continuous piecewise linear space, using the \code{-proj} option.
  Conversely, the 3D visualization bases on ellipsoids:
\begin{verbatim}
  ./stress < cube-P1.field > cube-stress-P1.field
  field cube-stress-P1.field -proj -stereo
\end{verbatim}
  \begin{figure}[htb]
     \begin{center}
       \begin{tabular}{cccc}
	    \includegraphics[width=6cm]{embankment-2d-s01-P0-fig.pdf} &
	    \hspace{0.0cm} &
	    \includegraphics[width=8cm]{embankment-2d-s01-P1d-elevation-fig.png} &
		\\
	    \includegraphics[width=7cm]{embankment-3d-s01-P0-fig.png} &
	    \hspace{0.0cm} &
	    \includegraphics[width=7cm]{embankment-3d-s01-P1d-fig.png} &
	    \stereoglasses
       \end{tabular}
     \end{center}
     \caption{The $\sigma_{01}$ stress component (linear elasticity $\lambda=1$):
        $d=2$ (top) and $d=3$ (bottom) ;
	$P_0$ (left) and $P_1$ discontinuous approximation (right).
     }
     \label{fig-embankment-s01}
  \end{figure}
\apindex{discontinuous}%
\pindexopt{field}{-elevation}%
You can observe a discontinuous constant or piecewise linear representation
of the approximate stress component $\sigma_{01}$
(see Fig.~\ref{fig-embankment-s01}):
\begin{verbatim}
  field square-stress-P1.field -comp 01
  field square-stress-P2.field -comp 01 -elevation
  field square-stress-P2.field -comp 01 -elevation -stereo
\end{verbatim}
\pindexopt{field}{-stereo}%
\pindex{paraview}%
Note that the \code{-stereo} implies the \code{paraview} render:
this feature available with \code{paraview}.
The approximate stress field can be also projected on
a continuous piecewise space:
\begin{verbatim}
  field square-stress-P2.field -comp 01 -elevation -proj
\end{verbatim}
  The tridimensional case writes simply (see Fig.~\ref{fig-embankment-s01}):
\pindexopt{field}{-comp}%
\begin{verbatim}
  ./stress < cube-P1.field > cube-stress-P1.field
  ./stress < cube-P2.field > cube-stress-P2.field
  field cube-stress-P1.field -comp 01 -stereo
  field cube-stress-P2.field -comp 01 -stereo 
\end{verbatim}
and also the P1-projected versions write:
\begin{verbatim}
  field cube-stress-P1.field -comp 01 -stereo -proj
  field cube-stress-P2.field -comp 01 -stereo -proj
\end{verbatim}
  These operations can be repeated for each $\sigma_{ij}$ components
  and for both \code{P1} and \code{P2} approximation of the
  deformation field.

% =======================
\subsection{Mesh adaptation}
% =======================
\cindex{mesh!adaptation!anisotropic}%
\findex{adapt}%
The main principle of the auto-adaptive mesh
writes~\citep{Hec-2006-bamg,BorGeoHecLauMohSal-1995,CasHecMohPir-1997,Val-1990,RoqMicSar-2000}:
\begin{verbatim}
      din >> omega;
      uh = solve(omega);
      for (unsigned int i = 0; i < n; i++) {
          ch = criterion(uh);
          omega = adapt(ch);
          uh = solve(omega);
      }
\end{verbatim}
The initial mesh is used to compute a first solution.
The adaptive loop compute an {\em adaptive criterion},
denoted by \code{ch},
that depends upon the problem under consideration and the
polynomial approximation used.
Then, a new mesh is generated, based on this criterion.
A second solution on an adapted mesh can be constructed.
The adaptation loop converges generally in roughly 5 to 20 iterations.

Let us apply this principle to the elasticity problem.

% -----------------------------------
\myexamplelicense{embankment_adapt.cc}
\myexamplelicense{elasticity_solve.icc}
\myexamplelicense{elasticity_criterion.icc}
\clindex{odiststream}
% -----------------------------------

\subsubsection*{Comments}

The criterion is here:
\[
  c_h = \left\{
    \begin{array}{ll}
	|{\bf u}_h|                            & \mbox{ when using } P_1 \\ 
	(\sigma({\bf u}_h):D({\bf u}_h))^{1/2} & \mbox{ when using } P_2
    \end{array}
  \right.
\]
\findex{interpolate}%
The \code{elasticity_criterion} function compute it as
\begin{lstlisting}[numbers=none,frame=none]
  return interpolate (Xh, norm(uh));
\end{lstlisting}
\findex{sqr}%
\findex{norm2}%
when using $P_1$, and as
\begin{lstlisting}[numbers=none,frame=none]
  return interpolate (T0h, sqrt(2*norm2(D(uh)) + lambda*sqr(div(uh))));
\end{lstlisting}
when using $P_2$.
The \code{sqr} function returns the square of a scalar.
Conversely, the \code{norm2} function returns the square of the norm.
%
\findex{adapt}%
In the min program, the result of the \code{elasticity_criterion} function is send to
the \code{adapt} function:
\begin{lstlisting}[numbers=none,frame=none]
  field ch = elasticity_criterion (lambda, uh);
  omega = adapt (ch, options);
\end{lstlisting}


\clindex{adapt_option}%
The \code{adapt_option} declaration is used by \Rheolef\  to
send options to the mesh generator.
The \code{err} parameter controls the error via the edge length of the mesh:
the smaller it is, the smaller the edges of the mesh are.
In our example, is set by default to one.
Conversely, the \code{hmin} parameter controls minimal edge length.

\subsubsection*{How to run the program}

  \begin{figure}[htb]
     \begin{center}
       \begin{tabular}{ccc}
	  $P_1$: 6771 elements, 3663 vertices
	  &&
	  $P_2$: 1636 elements, 920 vertices
	  \\
	  \includegraphics[height=7cm]{embankment-adapt-2d-P1-fig.png}
	  &
	  \hspace{0.3cm}
	  &
	  \includegraphics[height=7cm]{embankment-adapt-2d-P2-fig.png}
       \end{tabular}
     \end{center}
     \caption{Adapted meshes: the deformation visualization for $P_1$ and $P_2$ approximations.}
     \label{fig-embankment-adapt-deformation}
  \end{figure}
The compilation command writes:
\begin{verbatim}
    make embankment_adapt
\end{verbatim}
The mesh loop adaptation is initiated from a \code{bamg} mesh
(see also appendix~\ref{sec-bamg}).
\pindex{bamg}%
\exindex{square.bamgcad}%
\exindex{square.dmn}%
\fiindex{\filesuffix{.bamg} bamg mesh}%
\fiindex{\filesuffix{.bamgcad} bamg geometry}%
\begin{verbatim}
    bamg -g square.bamgcad -o square.bamg
    bamg2geo square.bamg square.dmn > square.geo
    ./embankment_adapt square P1 2e-2
\end{verbatim}
The last command line argument is the target error.
The code performs a loop of five mesh adaptations: the corresponding meshes
are stored in files, from \code{square-001.geo.gz} to \code{square-005.geo.gz},
and the associated solutions in files, from
\code{square-001.field.gz} to \code{square-005.field.gz}.
\pindex{gzip}%
\fiindex{\filesuffix{.gz} gzip compressed file}%
The additional \filesuffix{.gz} suffix expresses that the files are
compressed using \code{gzip}.
\pindexopt{field}{-nofill}%
\begin{verbatim}
    geo square-005.geo
    field square-005.field -nofill
\end{verbatim}
Note that the \filesuffix{.gz} suffix is automatically
assumed by the \code{geo} and the \code{field} commands.

For a piecewise quadratic approximation:
\begin{verbatim}
    ./embankment_adapt square P2 5e-3
\end{verbatim}
Then, the visualization writes:
\begin{verbatim}
    geo square-005.geo
    field square-005.field -nofill
\end{verbatim}
A one-dimensional mesh adaptive procedure is also possible:
\exindex{line.mshcad}%
\pindex{gmsh}%
\fiindex{\filesuffix{.msh} gmsh mesh}
\fiindex{\filesuffix{.mshcad} gmsh geometry}
\begin{verbatim}
    gmsh -1 line.mshcad -format msh2 -o line.msh
    msh2geo line.msh > line.geo
    geo line.geo
    ./embankment_adapt line P2
    geo line-005.geo
    field line-005.field -comp 0 -elevation
\end{verbatim}
The three-dimensional extension of this mesh adaptive procedure is in development.