\subsection{The incompressible elasticity problem}
\label{sec-incompressible-elasticity}
\cindex{problem!elasticity!incompressible}%
\cindex{boundary condition!Dirichlet}%
\cindex{boundary condition!Neumann}%
\cindex{boundary condition!mixed}%
\cindex{stabilization}%

\subsubsection*{Formulation}

  Let us go back to the linear elasticity problem.

  When $\lambda$ becomes large, this problem is related
  to the \emph{incompressible elasticity} and
  cannot be solved as it was previously done.
  To overcome this difficulty, the pressure is introduced~:
  \[
	p = - \lambda {\rm div}\, {\bf u}
  \]
  and the problem becomes:

  \ \ \ $(E)$ {\it find ${\bf u}$ and $p$ defined in $\Omega$ such that:}
  \[
    \begin{array}{ccccl}
      -\ {\bf div}(2D({\bf u})) &+& \bnabla p &=& {\bf f} \ {\rm in}\ \Omega, \\
      -\ {\rm div}\,{\bf u}     &-& \Frac{1}{\lambda} p       &=& 0 \ {\rm in}\ \Omega, \\
	+ B. C.
    \end{array}
  \]
\cindex{form!{$2D({\bf u}):D({\bf v})$}}%
\cindex{form!{${\rm div}({\bf u})\,q$}}%
  The variational formulation of this problem expresses:
 
  \ \ \ $(VFE)$ {\it find ${\bf u}\in V(1)$ and $p \in L^2(\Omega)$ such that:}
  \[
     \begin{array}{lcccl}
      a({\bf u},{\bf v}) &+& b({\bf v}, p) &=& m({\bf f},{\bf v}), \ \forall {\bf v}\in V(0), \\
      b({\bf u},q) &-& c(p,q)&=& 0, \ \forall q \in L^2_0(\Omega),
     \end{array}
  \]
  where
  \begin{eqnarray*}
       m({\bf u},{\bf v}) &=& \int_\Omega {\bf u} . {\bf v} \, dx, \\
       a({\bf u},{\bf v}) &=& \int_\Omega D({\bf u}) : D({\bf v}) \, dx, \\
       b({\bf v}, q) &=& - \int_\Omega {\rm div}({\bf v}) \, q \, dx.  \\
       c(p, q) &=& \Frac{1}{\lambda} \int_\Omega p \, q \, dx.  \\
	V &=& \{ {\bf v} \in (H^1(\Omega))^2; \  
      		{\bf v} = 0 \ {\rm on}\ \Gamma_{left} \cup \Gamma_{bottom} \}
  \end{eqnarray*}
  When $\lambda$ becomes large, we obtain the incompressible
  elasticity problem, that coincides with the Stokes problem.

\subsubsection*{Approximation}
\apindex{P1}%
\apindex{P2}%
  As for the Stokes problem,
  the \citet*{hood-taylor-73} finite element approximation
  is considered.
  We introduce a mesh ${\cal T}_h$ of $\Omega$
  and the following finite dimensional spaces:
  \begin{eqnarray*}
      X_h &=& \{ {\bf v} \in (H^1(\Omega)); \
          {\bf v}_{/K} \in (P_2)^2, \
          \forall K \in {\cal T}_h \}, \\
      V_h(\alpha) &=& X_h \cap V, \\
      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*}
  The approximate problem writes:

  \ \ \ $(VFE)_h$ {\it find ${\bf u}_h \in V_h(1)$ and $p \in Q_h$ such that:}
  \[
    \begin{array}{lcccl}
      a({\bf u}_h,{\bf v}) &+& b({\bf v}, p_h) &=& 0, \ \forall {\bf v}\in V_h(0), \\
      b({\bf u}_h,q) &-& c(p,q) &=& 0, \ \forall q \in Q_h.
    \end{array}
  \]

% --------------------------------------
\myexamplelicense{incompressible-elasticity.cc}
% --------------------------------------

% ======================
\subsubsection*{Comments}
% ======================
\cindex{method!conjugate gradient algorithm}%
\cindex{preconditioner!for nearly incompressible elasticity}%
  The problem admits the following matrix form:
  \[
     \left( \begin{array}{cc} 
	{\tt a} & {\tt trans(b)} \\ 
        {\tt b} & -{\tt c} 
     \end{array} \right)
     \left( \begin{array}{c} 
        {\tt uh} \\ 
        {\tt ph} 
     \end{array} \right)
	=
     \left( \begin{array}{c} 
	{\tt lh} \\
        {\tt 0}
     \end{array} \right)
  \]
  The problem is similar to the Stokes one (see page~\pageref{ref-pcg-abtb}).
  This system is solved by:
\clindex{problem_mixed}%
\begin{lstlisting}[numbers=none,frame=none]
  problem_mixed elasticity (a, b, c);
  elasticity.solve (lh, field(Qh,0), uh, ph);
\end{lstlisting}
  For two-dimensional problems, a direct solver is used by default.
  In the three-dimensional case, an iterative algorithm is the default:
  the preconditioned conjugate gradient.
  The preconditioner is here the mass matrix \code{mp} for the
  pressure. 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
  and is uniformly bounded when $\lambda$ becomes small,
  i.e. in the incompressible case.
% ======================
\subsubsection*{How to run the program}
% ======================
  \begin{figure}[htb]
     %\begin{center}
       \mbox{}\hspace{-0.5cm}
       \begin{tabular}{ccc}
	  \includegraphics[height=7cm]{incompressible-elasticity-square-fig.pdf} &
	  \includegraphics[height=7.5cm]{incompressible-elasticity-cube-fig.png} &
	  \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{The incompressible linear elasticity ($\lambda=+\infty$) for $N=2$ and $N=3$.}
     \label{fig-incompressible-elasticity-deformation}
  \end{figure}
  We assume that the previous code is contained in
  the file \reffile{incompressible-elasticity.cc}.
  Compile the program as usual (see page~\pageref{makefile}):
\begin{verbatim}
  make incompressible-elasticity
\end{verbatim}
and enter the commands:
\pindexopt{field}{-scale}%
\begin{verbatim}
  mkgeo_grid -t 10 > square.geo
  ./incompressible-elasticity square.geo 0 > square.field
  field square.field -nofill

  mkgeo_grid -T 10 > cube.geo
  ./incompressible-elasticity cube.geo 0 > cube.field
  field cube.field -fill -scale 2
\end{verbatim}
  The visualization is performed as usual: see section~\ref{sec-howtorun-elasticity},
  page~\pageref{sec-howtorun-elasticity}.
  Compare the results on {\sc Fig}.~\ref{fig-incompressible-elasticity-deformation},
  obtained for $\lambda=+\infty$ with those of
  {\sc Fig}.~\ref{fig-embankment-deformation}, page~\pageref{fig-embankment-deformation},
  obtained for $\lambda=1$.

  Finally, the stress computation and the mesh adaptation loop is left as an exercise to the reader.
