\documentclass[reqno]{amsart}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
2004-Fez conference on Differential Equations and Mechanics \newline
{\em Electronic Journal of Differential Equations},
Conference 11, 2004, pp. 33--40.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu  (login: ftp)}
\thanks{\copyright 2004 Texas State University - San Marcos.}
\vspace{9mm}}
\setcounter{page}{33}

\begin{document}

\title[\hfilneg EJDE/Conf/11 \hfil mixed finite elements for crack problems]
{Mixed finite element discretization of some variational
inequalities arising in elasticity problems in domains with
cracks}

\author[Z. Belhachmi, S. Tahir\hfil EJDE/Conf/11 \hfilneg]
{Zakaria Belhachmi, Souad Tahir}  % in alphabetical order

\address{Zakaria Belhachmi \hfill\break
Laboratoire de Math\'ematiques (LMAM),
Universit\'e de Metz, ISGMP,
Batiment A, Ile du Saulcy, 57045 Metz, France.}
\email{belhach@poncelet.sciences.univ-metz.fr}

\address{Souad Tahir \hfill\break
Laboratoire de Math\'ematiques (LMAM),
Universit\'e de Metz, ISGMP,
Batiment A, Ile du Saulcy, 57045 Metz, France}
\email{tahir@poncelet.sciences.univ-metz.fr}

\date{}
\thanks{Published October 15, 2004.}
\subjclass[2000]{65N30, 74M15, 35J85}
\keywords{Crack problems; variational inequalities; smooth domain
method; \hfill\break\indent
augmented Lagrangian formulation; mixed finite elements;
a priori estimates}


\begin{abstract}
  We consider some mixed variational formulations
  of elasticity system in domains with cracks.
  Inequality type conditions are prescribed at the crack faces
  which results in a model of unilateral contact.
  Relying in a new variational formulation of these problems
  in the smooth domain, we study and implement various mixed
  finite elements methods. We derive convergence rates and
  optimal error estimates.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}


\section{Introduction}

The numerical approximation of solutions of partial differential
equations in non smooth domains is often a very difficult and
challenging task. Among them, unilateral contact crack problems in
elasticity, arises some specific difficulties, both in theoretical
and approximation grounds. Such problems are characterized by
inequality type and nonlinear boundary conditions prescribed on
non smooth part of the boundary \cite{KK}, describing the mutual
non penetration between the crack faces.

The variational formulations associated to these problems lead to
solving variational inequalities arising in contact mechanics. We
refer the reader to \cite{DuL,KS,KO} for mathematical foundations
and \cite{KK} for crack problems. A new approach to crack theory
for linear elastic bodies proposed in \cite{SK} allows us to solve
the problem in the entire domain (including the cracks), reducing
the difficulties which occur when dealing with the numerical
simulation of such problems. In the framework of this formulation,
the discretization with finite elements method leads to similar
results, from the accuracy point of view and for developing
efficient algorithms, as for classical variational inequalities in
unilateral contact problems \cite{BSS2}. For the discretization of
variational inequalities, we refer the reader to
 \cite{BHR,GL,HHN1,HHN,KO, Z}, and for more recent
developments to \cite{BBR,CHLS,BLS,BBZ}.

In \cite{BSS1,BSS2}, the approach called smooth domain method is
considered in the case of an elastic membrane. Using a
regularization technique, the discretization by various mixed
finite element methods is considered, the numerical analysis of
the method is carried out and numerical simulations are performed.

We extend in this work the formulation in the smooth domain to the
case of the elasticity system. However, our discretization by
finite elements method is based on an augmented Lagrangian
formulation \cite{GLT} since this approach yields many efficient
numerical algorithms.

\section{The continuous problem}

\subsection*{Problem formulation}

Let $\Omega$ be a bounded domain in ${\mathbb R}^2$ with smooth boundary
$\Gamma$, and $\Gamma_c\subset\Omega$ be a smooth curve without
selfintersections.
We assume that $\Gamma_c$ can be extended to
a closed smooth curve $\Sigma\subset\Omega$, with $\Sigma$ of class $C^{1,1}$,  and
$\Omega=\Omega^1\cup\Sigma\cup\Omega^2$ divided into two sub-domains $\Omega^1$, $\Omega^2$. In this
case, $\Sigma=\partial\Omega^1 $ is the boundary of $\Omega^1$ and $\Sigma\cup\Gamma=\partial\Omega^2$
is the boundary of $\Omega^2$. Let $\Omega_c$ be the domain $\Omega\setminus\overline{\Gamma_c}$,
then $\Gamma_c$ is called a crack in the elastic body of the
reference configuration $\Omega_c$ (see Figure \ref{fg:fig1-1}).

The equilibrium problem for a linear elastic body occupying the
domain $\Omega_c$ with the interior crack $\Gamma_c$ can be formulated
as follows: find $\textbf{u}=(u_1,u_2)$, and
$\sigma=(\sigma_{ij})$, $i,j=1,2$, such that
\begin{gather}
  -\textbf{div}\,\sigma  = \textbf{f}   \quad\text{in  }\Omega_c, \label{eq:2-1}  \\
  C\sigma-\epsilon (\textbf{u}) = 0  \quad \text{in  }\Omega_c,  \label{eq:2-2} \\
  \textbf{u} = 0 \quad \text{on  }\Gamma,   \label{eq:2-3} \\
  \left[\textbf{u}\right]\nu\geq 0, \quad
  [\sigma\,\nu]  = 0,
  \quad \sigma_\nu [\textbf{u}] =0  \quad    \text{on }\Gamma_C, \label{eq:2-4} \\
  \sigma_\nu\leq 0, \quad \sigma_t =0
 \quad \text{on }\Gamma_C^{\pm}, \label{eq:2-5}.
  \end{gather}
Here $[\textbf{u}] = \textbf{u}^{+}-\textbf{u}^{-}$ denotes the
jump of the displacement field across $\Gamma_c$, and the signs
$\pm$ indicate the positive and negative directions of the normal
$\nu$. $\textbf{f}=(f_1,f_2)\in L^2(\Omega)^2$ is a given external
force acting on the body. We have used the following standard
notation:
\begin{gather*}
\sigma_\nu=\sigma_{ij}\nu_j\nu_i,\ \sigma_t=\sigma\,\nu-
\sigma_\nu\,\nu=\left\{\sigma_t^i\right\}_{i=1}^2,\
\sigma\nu=\left\{\sigma_{ij}\nu_j\right\}_{i=1}^2,\\
\epsilon_{ij}(\textbf{u})=\frac{1}{2}(u_{i,j}+u_{j,i}),\
i,j=1,2,\ \epsilon (\textbf{u})=(\epsilon_{ij})_{i,j=1}^2,\\
\left\{ C\sigma\right\}_{ij}=c_{ijk\ell}\sigma_{k\ell},\
c_{ijk\ell}=c_{jik\ell}=c_{k\ell ij},\
c_{ijk\ell}\in L^\infty(\Omega).
\end{gather*}
The tensor $C$ satisfies the ellipticity condition
\begin{equation}\label{eq:ellip1}
c_{ijk\ell}\xi_{ji}\xi_{k\ell}\geq c_0\vert\xi\vert^2,\quad
\forall \xi_{ji}=\xi_{ij},\; c_0 > 0.
\end{equation}
We use the summation convention over repeated indices.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.5\textwidth]{fig1}
\caption{A crack in the reference configuration. \label{fg:fig1-1}}
\end{center}
\end{figure}

\subsection*{The smooth domain method}

The smooth domain method is based on a mixed variational
formulation of problem \eqref{eq:2-1}-\eqref{eq:2-5}.
We consider the space
\begin{equation*}
\mathbf{X}=\left\{\sigma=\left\{\sigma_{ij}\right\}\in (L^2(\Omega_c))^4 :
\textbf{div}\,\sigma\in (L^2(\Omega_c))^2\right\},
\end{equation*}
equipped with the norm
\begin{equation*}
\Vert \sigma\Vert_{\mathbf{X}}=\left(\Vert
\sigma\Vert_{(L^2(\Omega_c))^4}^2+ \Vert\textbf{div}\,
\sigma\Vert_{(L^2(\Omega_c))^2}^2\right)^{1/2},
\end{equation*}
and we define the closed convex set
\[
\mathbf{K}=\left\{\sigma\in\mathbf{X} : [ \sigma.\nu]=0,\
\ \text{on}\,\Gamma_c,\ \sigma_\nu\leq 0,\,\sigma_t =0
\ \text{on}\,\Gamma_c^{\pm}\right\}
\]
where $[ \cdot ]$ denotes the jump across $\Gamma_c$.

\begin{remark} \rm
The contact condition on $\Gamma_c$ is to be understood in the
weak sense of traces by using fractional Sobolev spaces
$H_{00}^{1/2}(\Gamma_c)$ and its dual space \cite{LM}.
\end{remark}

The mixed formulation
for problem \eqref{eq:2-1}-\eqref{eq:2-5} reads: Find
$\textbf{u}=(u_1,u_2)\in (L^2(\Omega_c))^2$,
$\sigma=\left\{\sigma_{ij}\right\}\in\mathbf{K}$, such that
\begin{equation}\label{eq:Fhyb}
\begin{gathered}
(C\sigma,\tau-\sigma)_{\Omega_c} + (\textbf{u},
\textbf{div}\,\tau-
\textbf{div}\sigma)_{\Omega_c}\geq 0,\quad\forall
\overline{\sigma}\in\mathbf{K}, \\
 \textbf{div}\,\sigma = -\textbf{f}, \quad\text{in}\ \Omega_c.
\end{gathered}
\end{equation}
The well-posedness of this problem  follows from the
Brezzi-Babuska theory \cite{BF}. A proof based on a regularization
argument is given in \cite{SK} (see also \cite{T} and \cite{BLS}); we
state only the result.

\begin{proposition}\label{pro:wp1}
There exists a unique solution to problem \eqref{eq:Fhyb}.
\end{proposition}

The formulation in the smooth domain consists of extending
$\sigma$ and $\mathbf{u}$ to the whole domain
$\Omega=\Omega_c\cup\overline\Gamma_c$, and it results in the (close) formulation
obtained by replacing $\Omega_c$ with $\Omega$ in \eqref{eq:Fhyb}, with
the obvious modification of the spaces that we denote also $\mathbf{X}
(\Omega)$ and $\mathbf{K} (\Omega)$ (see. \cite{SK} for details). This
formulation allows us to solve the problem in the smooth geometric
domain while the constraint at the crack faces is included in the
functional space.

The extended problem is also well-posed and it is readily checked
that the restriction of the solution $(\mathbf{u},\sigma)$ to $\Omega_c$ is
the solution of problem \eqref{eq:Fhyb}. The converse  is
(trivially) true if additional  smoothness holds for solutions of
problem \eqref{eq:Fhyb}.

The Discretization of problem \eqref{eq:Fhyb} by various finite
elements of low order based on affine finite elements and
Raviart-Thomas elements, is considered in \cite{T}. In particular,
the study of approximation properties as well as numerical
simulations of the discrete problem is performed. In what follows,
we consider a new formulation based on an augmented Lagrangian
technique.

\section{Augmented Langrangian formulation}

\subsection*{Mixed variational formulation}

The main advantage of the augmented Lagrangian method, compared to
the method used in \cite{BSS2}, is that it converges without
imposing to the penalization parameter to be too small, so the
numerical resolution is more stable and efficient. The new
formulation is described by the following problem: For $r >0$,
find $\textbf{u}=(u_1,u_2)\in (L^2(\Omega))^2$,
$\sigma=\{\sigma_{ij}\}\in\mathbf{K}$, such that for all
$\tau\in\mathbf{K}$,
\begin{equation}\label{eq:Fhybal}
\begin{gathered}
(C\sigma,\tau-\sigma) + (\textbf{u},
\textbf{div}\,\tau - \textbf{div}\sigma) +  r
(\textbf{div}\,\sigma, \textbf{div}\,\tau- \textbf{div}\sigma)
\geq r\,(-\textbf{f},\textbf{div}\,\tau-\textbf{div}\sigma), \\
 \textbf{div}\,\sigma = -\textbf{f},  \text{in}\ \Omega.
\end{gathered}
\end{equation}
Clearly, this problem is equivalent to the minimization
over $L^2(\Omega)^2\times\mathbf{K}$ of the functional
\begin{equation*}
{\mathcal L} (\mathbf{v},\tau) = \frac{1}{2} (C\tau,\tau) + (\mathbf{v},\textbf{div}\,\tau+\textbf{f})
+\frac{r}{2}\Vert\textbf{div}\,\tau+\textbf{f}\Vert_{L^2(\Omega)^2}^2.
\end{equation*}
The bilinear form
\[
a_r(\sigma,\tau)=(C\sigma,\tau)+r(\textbf{div}\,\sigma,
\textbf{div}\,\tau),
\]
is elliptic on $\mathbf{K}$ and the bilinear form $b(.,.)$ satisfies the
usual Brezzi-Babuska condition (\cite{BF}). The well-posedness of
problem \eqref{eq:Fhybal} and the convergence of the solution
$(\sigma(r),\mathbf{u}(r))$ to the solution $(\sigma,\mathbf{u})$ of problem
\eqref{eq:Fhyb} (see \cite{T}), follows from standard elliptic variational
inequalities \cite{KS,GL} by applying Stamppacchia's
theorem on the cone $\mathbf{K}\times L^2(\Omega)$. Thus, we have the
following result.

\begin{proposition}\label{pro:wp2}
There exists a unique solution $(\sigma(r),\mathbf{u}(r))$ to problem \eqref{eq:Fhybal}.
Moreover, $(\sigma(r),\mathbf{u}(r))$ converges to
 $(\sigma,\mathbf{u})$, the solution of \eqref{eq:Fhyb} when $r$ goes to zero.
\end{proposition}

\subsection*{Mixed hybrid formulation}

The algorithms used to solve the discrete problem corresponding to
\eqref{eq:Fhybal} are the steepest descent methods. Another
strategy for solving the problem consists in solving a hybrid
formulation where the constraint in $\mathbf{K}$ is taken into account
with a Lagrange multiplier. The resulting algorithm yields to
solve a quadratic programming problem of small size to compute the
Lagrange multiplier and then a Large linear system to compute the
other unknowns. This approach is based on the following mixed
hybrid formulation. For the sake of brevity, we will denote by
$\mathbf{M}$ the space $H^{1/2}(\Gamma_c)^2$. We introduce also
the following closed convex cone
\[
M_+ = \left\{ \mu\in H^{1/2}(\Gamma_c);\ \mu\geq 0, \text{a.e.}\right\}.
\]
and we define the new functional $\tilde{\mathcal L}$ over
$\mathbf{V}\times\mathbf{X}\times \mathbf{M}\times M_+$
\[
\tilde{{\mathcal L}}(\mathbf{v},\tau,\mu_t,\mu_n) = {\mathcal L}(\mathbf{v},\tau)
 + \langle\langle \tau_t,\mu_t\rangle\rangle _{1/2 , \Gamma_c}
 + \langle\tau_\nu,\mu_n\rangle_{1/2 ,  \Gamma_c}.
\]
where $\langle \cdot,\cdot\rangle_{1/2 ,\Gamma_c}$ denotes the duality product
 between $(H^{1/2}(\Gamma_c))$ and its dual space.
$\langle\langle \cdot,\cdot \rangle\rangle _{1/2 ,\Gamma_c}$ is the duality
product defined by
\[
\langle\langle \tau_t,\phi\rangle\rangle _{1/2 ,\Gamma_c}
=\langle \tau_{t1},\phi_1\rangle _{1/2 ,\Gamma_c}
+\langle \tau_{t2},\phi_2\rangle_{1/2 ,\Gamma_c},\quad
\forall \phi=(\phi_1,\phi_2)\in\mathbf{M},\ \phi_i\nu_i=0.
\]
Therefore, problem \eqref{eq:Fhybal} can be written as:
Find
$(\mathbf{u},\sigma,\lambda_t, \lambda_n)\in \mathbf{V}\times\mathbf{X}\times \mathbf{M}\times
M_+$, $\lambda_{ti}\nu_i=0$, such that
\begin{equation}\label{eq:Fhybaldm}
\begin{gathered}
a_r(\sigma,\tau)+b(\mathbf{u},\tau)
+ \langle\langle \tau_t,\lambda_t\rangle\rangle _{1/2 ,\Gamma_c}
+ \langle\tau_\nu,\lambda_n\rangle_{1/2 ,\Gamma_c} =
-r (\textbf{f},\textbf{div}\,\tau),\quad\forall \tau\in\mathbf{X}, \\
b(\mathbf{v},\sigma) = -(\textbf{f},\mathbf{v}), \quad\forall\mathbf{v}\in\mathbf{V}, \\
\langle\langle \sigma_t,\mu_t\rangle\rangle _{1/2 ,\Gamma_c}= 0,
 \quad\forall\mu_t=(\mu_1,\mu_2)\in \mathbf{M},\ \mu_i\nu_i=0, \\
\langle \sigma_\nu,\mu_n-\lambda_n\rangle _{1/2 ,\Gamma_c}\leq 0,
 \quad\forall\mu\in M_+.
\end{gathered}
\end{equation}
It is readily checked that solution
$(\mathbf{u},\sigma,\lambda_t,\lambda_n)$ of problem \eqref{eq:Fhybaldm}
is a saddle-point of the Lagrangian $\tilde{\mathcal L}$.

\section{Discrete mixed hybrid problem}

We denote by $T_h$ a triangulation of $\Omega$ made of elements which
are triangles (or quadrilateres) with a maximum size $h$
satisfying the usual admissibility assumption, i.e. the
intersection of two different elements is either empty,  a vertex,
or a whole edge. In addition, $\mathcal{T}_{h}$ is assumed regular, i.e.
the ratio of the diameter  of any element $T\in \mathcal{T}_{h}$ to the
diameter of its largest inscribed ball is bounded by a constant
$\sigma$ independent of $T$ and $h$. We will assume that the
endpoints of $\Gamma_c$ are vertices of the triangulation. The
set of nodes on $\Gamma_c$ are denoted by
$c_1=x_0,x_1,\dots,x_{I-1},x_I=c_2$ and we set
$t_i=]x_{i-1},x_i[$.

In \cite{T} we study various discretizations based on the
finite elements method for both problem \eqref{eq:Fhybal} and problem
\eqref{eq:Fhybaldm}.
In particular, we perform the complete analysis and
the numerical simulations with
the (lower order) PEERS element and the BDMS (Brezzi-Douglas-Marini and Steinberg)
elements.
The use of the above elements is based on a formulation obtained by modifying
the Hellinger-Reissner principle. In this approach, the symmetry of the stress
tensor $\sigma$ is relaxed and only imposed by means of the Lagrange multiplier.

We present here another example of finite elements also studied in \cite{T} due to
 Arnold and Winther \cite{AW}, based on the unaltered Hellinger-Reissner principle.

For any $T\in \mathcal{T}_h$, and for any Banach space $E$, we denote by $P_k(T,E)$ the
space of polynomial functions over $T$ with values in $E$,  of degree less than $k$.
As usual we denote by $P_k(T)$ the space $P_k(T,{\mathbb R})$. We define
the associated local space $\Sigma_T$ as
\begin{align*}
\Sigma_T &=P_2(T,{\mathbb R}_{\text{sym}}^{2\times 2})+\left\{\tau\in P_3(T,{\mathbb R}_{\text{sym}}^{2\times 2});\ \text{div}\,\tau=0\right\} \\
& = \left\{\tau\in P_3(T,{\mathbb R}_{\text{sym}}^{2\times 2});
\ \text{div}\,\tau\in P_1(T)^2\right\}\,.
\end{align*}
We introduce the
following discrete spaces: For $h>0$,
\begin{gather*}
\mathbf{X}_{h}=\left\{ \sigma_h\in\mathbf{X}; \sigma_{h\vert T}\in
\Sigma_T,\ \forall T\in\mathcal{T}_h\right\}, \\
\mathbf{V}_h=\bigl\{ \mathbf{v}_h\in (\mathcal{C}(\overline\Omega))^2;\ \mathbf{v}_{h\vert T}\in
(P_1(T))^2,\ \forall T\in\mathcal{T}_h\bigr\}.
\end{gather*}
Concerning the approximation of Lagrange multipliers, we introduce
the space
\begin{equation*}
W_h^1(\Gamma_c)=\left\{ \mu_h\in C(\overline{\Gamma}_c),\
\mu_{h\vert t_i}\in P_1(t_i),\ 0\leq i\leq I-1\right\}.
\end{equation*}
We also denote by $\mathbf{W}_h^1(\Gamma_c)$ the space $(W_h^1(\Gamma_c))^2$.

Therefore, we can define various version of the discrete space
$\mathbf{M}_h$ and the discrete convex cone $M_{h+}$. An example of
such choices consists to take \cite{T}
\begin{align*}
\mathbf{M}_h&=\mathbf{W}_h^1(\Gamma_c) \\
M_{h+}&=\big\{\mu_h\in W_h^1(\Gamma_c),\
\int_{\Gamma_c}\mu_h\psi_h\,d\Gamma\geq 0,\ \forall\psi_h\in
W_h^1, \psi_h\geq 0\big\}.
\end{align*}
We will denote for brevity ${\mathcal M}_h$ the convex cone
$\mathbf{M}_h\times M_{h+}$ and
$\lambda_h=(\lambda_{ht}=(\lambda_{ht1},\lambda_{ht2}),\lambda_{hn})$,
$\lambda_{hti}\nu_i=0$. The discrete problem is
now the same as problem \eqref{eq:Fhybaldm} when replacing the
unknowns and spaces by the finite dimensional analogues.

The following result is proved in \cite{T}.

\begin{proposition}\label{pro:cv}
Assume that the set $\mathbf{M}_h$ and $M_{h+}$ are given as above.
Then, the discrete problem corresponding to \eqref{eq:Fhybaldm}
admits a unique solution.
\end{proposition}

The error analysis and the study of convergence rates is based on
uniform inf-sup condition (with respect to $h$) (see \cite{BSS2})
we only states the main error estimate when the multiplier spaces
are chosen as in Proposition \ref{pro:cv} (see \cite{T} for details). Let
us define
\begin{equation*}
\mathbf{Z}=\left\{\sigma=\left\{\sigma_{ij}\right\}\in (H^1(\Omega))^4 :
\textbf{div}\,\sigma\in (H^1(\Omega))^2\right\},
\end{equation*}
and denotes by $\mathbf{Z}(\Omega^\ell)$, $\ell=1,2$ the space of the restrictions
of functions of $\mathbf{Z}$ to $\Omega^\ell$.
We have the following result.

\begin{theorem}\label{thm:approx}
Let $(\mathbf{u},\sigma,\lambda=(\lambda_{t1},\lambda_{t2},\lambda_n))$ be the solution of problem
\eqref{eq:Fhybaldm}. Suppose that $\mathbf{u}_{\vert\Omega^1}\in H^2(\Omega^1)$,
$\mathbf{u}_{\vert\Omega^2}\in H^2(\Omega^2)$ and also $\sigma_{\vert\Omega^1}\in
\mathbf{Z}(\Omega^1)$, $\sigma_{\vert\Omega^2}\in \mathbf{Z}(\Omega^2)$. Let
$(\mathbf{U}_h,\lambda_h=(\lambda_{ht1},\lambda_{ht2},\lambda_{hn}))$ be the solution of the discrete problem
associated to \eqref{eq:Fhybaldm} with the specified choice of
$\mathbf{M}_h$ and $M_{h+}$ given above. Then the following estimate
holds
\begin{equation}\label{eq:errMh0}
\begin{aligned}
&\Vert \mathbf{u}-\mathbf{u}_h\Vert_V +\Vert\sigma-\sigma_h\Vert_{\mathbf{X}} +
\Vert\lambda_t- \lambda_{ht}\Vert_{(H^{1/2}(\Gamma_c))^2}
+ \Vert\lambda_n- \lambda_{hn}\Vert_{H^{1/2}(\Gamma_c)}\\
&\leq C(r,u,\sigma,\lambda)h^{3/4}\, .
\end{aligned}
\end{equation}
The constant $C(r,\mathbf{u},\sigma,\lambda)$ depends linearly on $\Vert
\mathbf{u}_{\vert\Omega^\ell} \Vert_{H^2(\Omega^\ell)}$,
$\Vert\sigma_{\vert\Omega^\ell}\Vert_{H^1(\Omega^\ell)^4}$, and
 $\Vert \mathbf{div}\,\sigma_{\vert\Omega^\ell}\Vert_{H^1(\Omega^\ell)^2}$,
 $\ell =1,2$.
\end{theorem}

\begin{remark} \rm
More technical arguments allow us to avoid the regularity
assumption $\text{div}\,\sigma\in (H^1(\Omega^\ell))^2$. In this case $\sigma$
is estimated in the $L^2$-norm.
\end{remark}

\section{Implementation details}

To perform the computations for the mixed hybrid problem,
the matrix formulation of discrete problem \eqref{eq:Fhybaldm} is
derived. Let $\mathbf{V}$, $\mathbf{U}$ denote the vectors with the entries
given by the nodal values of the functions $(v_h,\tau_h)$ and
$(u_h,\sigma_h)$, respectively. Let $M$ and $\Lambda$ be the
vectors with the entries given by the nodal values of $\mu_h$ and
$\lambda_h$, respectively, for the various choices of the space
$\mathbf{M}_h$ and the convex set $M_{h+}$. Therefore, the saddle-point
problem for the Lagrangian can be rewritten
in finite dimensional setting:\\
Find $\mathbf{U}=(u_h,\sigma_h)$ and $\Lambda$, defined by the following
max-min condition
\begin{equation}\label{eq:sadelm}
\max_{SM\geq 0}\bigl( \min_{\mathbf{V}} \frac{1}{2}{}^t\mathbf{V} \mathbf{K}\mathbf{V}
- {}^t\mathbf{V}\mathbf{F}
+({}^t\mathbf{V}\,\mathbf{L})SM \bigr),
\end{equation}
where $\mathbf{K}$ denotes the stiffness matrix,  $\mathbf{F}$ is the vector
corresponding to the external loading and the matrix $S$ expresses
the sign conditions for the multipliers.

The solution $(\mathbf{U},\Lambda)$ of \eqref{eq:sadelm} satisfies
the saddle-point conditions and we have
\begin{equation}\label{eq:bu}
\mathbf{U} = \mathbf{K}^{-1} (\mathbf{F}-L\,S\Lambda).
\end{equation}
Therefore, for
$\Phi=S\Lambda$, the saddle-point problem \eqref{eq:sadelm} can
be rewritten as a quadratic programming problem
\begin{equation}\label{eq:minQP}
\min_{\Phi\geq 0}\bigl( \frac{1}{2}{}^t\Phi{}^tL\mathbf{K}^{-1}L\Phi-
{}^t\Phi{}^tL\mathbf{K}^{-1}\mathbf{F}+\frac{1}{2}{}^t\mathbf{F}\mathbf{K}^{-1}\mathbf{F}\bigr).
\end{equation}
If $\overline\Phi$ is the solution of \eqref{eq:minQP} then
$\Lambda=S^{-1}\overline\Phi$. The solution $\mathbf{U}$ is
obtained by solving \eqref{eq:bu}.

Several numerical simulations and experiments are performed in
\cite{T} for this augmented Lagrangian method applied to the new
formulation of crack problems in the smooth domain. The
computations with various mixed finite element discretizations,
for both unsymmetric stress tensor formulation (PEERS, BDMS) and
the symmetric one considered in this article, give satisfactory
results confirming the theoretical estimates. The use of the
augmented Lagrangian formulation (i.e. $r\not= 0$) is unnecessary
for practical computations when compatible pairs of finite
elements for the displacements and stresses are chosen, however it
allows using efficient algorithms to solve the discrete problem.


\begin{thebibliography}{00}

\bibitem{AW} D. N. Arnold, R. Winther; Mixed finite elements for elasticity,
Numer. Math., \textbf{92} (2002), 401-419.

\bibitem {BSS1} Z. Belhachmi, J. M. Sac-Ep\'ee, J. Sokolowski; 
Approximation
par la m\'ethode des \'el\'ements finis de la formulation en
domaine r\'egulier de probl\`emes de fissures, C. R. Acad. Sci. Paris,
S\'erie I \textbf{338} (2004), 499-504.

\bibitem {BSS2} Z. Belhachmi, J.M. Sac-Ep\'ee, J. Sokolowski; 
Mixed finite element methods for a smooth domain formulation of a
crack problem. submitted.

\bibitem {BBZ} Z. Belhachmi, F. Ben Belgacem; Quadratic finite element
for Signorini problem, Math. Comp. \textbf{72}, (2003), 83-104.

\bibitem {BHL}F. Ben Belgacem, P. Hild, P. Laborde; Extention
of the mortar finite element method to a variational
inequality modeling unilateral contact, Math. Models Methods.
Appl. Sci., \textbf{9}, (1999), 287-303.

\bibitem{BBR} F. Ben Belgacem, Y. Renard; Hybrid finite element methods for the Signorini problem, Math. Comput. \textbf{72}, (2003), 1117-1145.

\bibitem {BHR} F. Brezzi, W. W. Hager and P. A. Raviart; 
Error estimates for the finite element solution  of variational
inequalities, Numer. Math., \textbf{28}, (1977), 431-443.

\bibitem{BF} F. Brezzi, M. Fortin; 
Mixed and hybrid finite element methods, Springer Verlag, New York, 
Springer Series in Computational Mathematics, \textbf{15}, (1991).

\bibitem{CHLS}P. Coorevits, P. Hild, K. Lhalouani, T. Sassi; 
Mixed finite element methods for unilateral
problems: convergence analysis and numerical studies.
Math. Comp., \text{71}, 237, (2001), 1-25.

\bibitem {DuL} G. Duvaut, J.-L. Lions; 
\textit{Les in{\'e}quations en m{\'e}canique et en physique}, Dunod,
(1972).

\bibitem{GLT} R. Glowinski, P. Le Tallec;  
Augmented Lagrangian and operator-splitting  methods in nonlinear mechanics. 
Siam, Philadelphia.
%

\bibitem{GL} R. Glowinski;  Lectures on numerical methods
for nonlinear variational problems, Springer, Berlin, (1980).

\bibitem {HHN1} J. Haslinger and I. Hlav\' a\v{c}ek; 
\textit{Contact between Elastic Bodies -2.Finite Element
Analysis}, Aplikace Matematiky, \textbf{26}, (1981), 263--290.

\bibitem {HHN} J. Haslinger, I. Hlav\' a\v{c}ek, J. Ne\v{c}as; 
 Numerical Methods  for Unilateral Problems in Solid Mechanics,
 in the Handbook of Numerical Analysis, Vol. \textbf{IV}, Part \textbf{2},
 P.G. Ciarlet \& J.-L. Lions eds, North-Holland, (1996).

\bibitem{KK} A. M. Khludnev, V. A. Kovtunenko;  Analysis of
cracks in solids, Southampton-Boston, WIT press, (2000).

\bibitem{SK} A. M. Khludnev, J. Sokolowski;  Smooth domain method
for crack problems. to appear in Quarterly of Applied Mathematics.

\bibitem{KO} N. Kikuchi, J. Oden;  Contact problems in
elasticity: A study  of variational inequalities and finite element methods,
SIAM, 1988.

\bibitem {KS} D. Kinderlehrer, G. Stamppachia; 
An introduction to variational inequalities and their
applications, Academic Press, (1980).

\bibitem {LM} J.-L. Lions, E. Magenes; 
Probl{\`e}mes aux limites non homog{\`e}nes, Dunod, (1968).

\bibitem {BLS} L. Slimane, A. Bendali, P. Laborde; 
Mixed formulations for a class of variational inequalities, 
M2AN, \textbf{38}, 1, (2004), 177-201.

\bibitem {S} L. Slimane; 
M\'ethodes mixtes et traitement du verouillage num\'erique pour la r\'esolution des in\'equations variationnelles,
 Ph.D. Thesis, INSA Toulouse, France (2001).

\bibitem{T} S. Tahir; M\'ethodes d'approximation et analyse a posteriori des  in\'equations variationnelles et applications \`a des probl\`emes de fissures en \'elasticit\'e, Ph.D. Thesis, universit\'e de Metz, to be published.

\bibitem {Z} Z.-H. Zhong;
\textit{Finite Element Procedures for Contact-Impact Problems},
Oxford University Press, 1993.

\end{thebibliography}


\end{document}
