A staggered finite element cell-centered scheme for nearly incompressible elasticity on general meshes
Name\corref{cor1
%% Text of abstract We apply the finite element cell-centered (FECC) scheme \cite{PO12 to the solution of the nearly incompressible elasticity problem. By applying a technique of dual mesh, such a low-order finite element scheme can be constructed from any given mesh and it is proved to satisfy the ``macroelement condition" \cite{RS84, meaning that the stability condition is fulfilled. Numerical results show that the method, which is simple to implement, is effective in terms of accuracy and computational cost compared with other methods. %to investigate the performance of the method and to compare it with the other methods in terms of accuracy and computational cost.
1406.2073.tex
%% 
%% Copyright 2007, 2008, 2009 Elsevier Ltd
%% 
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%% 
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%% 
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%% 

%% Template article for Elsevier's document class `elsarticle'
%% with numbered style bibliographic references
%% SP 2008/03/01

\documentclass[preprint,12pt]{elsarticle}

%% Use the option review to obtain double line spacing
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}

%% For including figures, graphicx.sty has been loaded in
%% elsarticle.cls. If you prefer to use the old commands
%% please give \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}
%% The amsthm package provides extended theorem environments
\usepackage{amsthm}
\usepackage{amsmath}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers.
%% \usepackage{lineno}

% Author's macros
%
\newtheorem{remark}{Remark}
\newtheorem{theorem}{Theorem}
\newcommand{\bE}{\pmb{E}}
\newcommand{\bB}{\pmb{B}}
\newcommand{\bV}{\pmb{V}}
\newcommand{\bv}{\pmb{v}}
\newcommand{\bx}{\pmb{x}}
\newcommand{\bA}{\pmb{A}}
\newcommand{\bC}{\pmb{C}}
\newcommand{\bF}{\pmb{F}}
\newcommand{\bU}{\pmb{U}}
\newcommand{\bP}{\pmb{P}}
\newcommand{\bsigma}{\pmb{\sigma}}
\newcommand{\bff}{\pmb{f}}
\newcommand{\bu}{\pmb{u}}
\newcommand{\iT}{\mathcal{T}}
\newcommand{\iB}{\mathcal{B}}
\newcommand{\iK}{\mathcal{K}}
\newcommand{\iN}{\mathcal{N}}
\newcommand{\iC}{\mathcal{C}}
\newcommand{\mR}{\mathbb{R}}
\newcommand{\Div}{\text{div}\;}
%



%\journal{} %{Applied Mathematics and Computation}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for theassociated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for theassociated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for theassociated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}

\title{A staggered finite element cell-centered scheme for nearly incompressible elasticity on general meshes}

%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{}
%% \address[label1]{}
%% \address[label2]{}

\author[Phuong]{T.T.P. Hoang}
\address[Phuong]{Department of Mathematics, Ho Chi Minh city University of Pedagogy, 280 An Duong Vuong street, ward 4, district 5, Ho Chi Minh city, Vietnam \vspace{5pt}}
\author[aHai]{Ong Thanh Hai}
\address[aHai]{Faculty of Mathematics and Computer Science, University of Science, VNU HCMC, 227 Nguyen Van Cu Street, District 5, Ho Chi Minh City, Vietnam.}  
\author[aHai]{ H. Nguyen-Xuan}
%\address[aHung]{Faculty of Mathematics and Computer Science, University of Science, VNU HCMC, 227 Nguyen Van Cu Street, District 5, Ho Chi Minh City, Vietnam.}
\begin{abstract}
%% Text of abstract
We apply the finite element cell-centered (FECC) scheme \cite{PO12} to the solution of the nearly incompressible elasticity problem. By applying a technique of dual mesh, such a low-order finite element scheme can be constructed from any given mesh and it is proved to satisfy the ``macroelement condition" \cite{RS84}, meaning that the stability condition is fulfilled. Numerical results show that the method, which is simple to implement, is effective in terms of accuracy and computational cost compared with other methods. 
%to investigate the performance of the method and to compare it with the other methods in terms of accuracy and computational cost. 
\end{abstract} 

\begin{keyword} nearly incompressible elasticity \sep finite elements \sep cell-centered scheme \sep macroelement condition \sep inf-sup condition
%% keywords here, in the form: keyword \sep keyword

%% PACS codes here, in the form: \PACS code \sep code

%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)

\end{keyword}

\end{frontmatter}

%% \linenumbers

%% main text
%
%
%
%
%
%
%
%
%
%
% ###################################################
%
%			INTRODUCTION
%
% ###################################################
\section{Introduction}
\label{Sec:intro}
%\begin{itemize}
%\item Motivation: new numerical method for nearly incompressible elasticity problems using low-order finite elements.  
%\item Important points: FECC scheme can deal with an arbitrary triangulation of the domain of calculation and give high accuracy with low computational cost; stability (the inf-sup condition) is obtained by using the ``macroelement condition".\\
%$\rightarrow$ compare our approach with Lamichane's paper: we do not use bubble functions.  
%\end{itemize}
In this paper, we present a new numerical method for nearly incompressible elasticity problems using low-order finite elements. 
The scheme, firstly introduced for solving stationary diffusion problems, has many advantages: i) it can deal with general meshes and it involves only cell unknowns; ii) it gives higher accuracy while the computational cost is the same as other cell-centered scheme; iii) the stability condition is obtained by using the ``macroelement condition" instead of using bubble functions (as in \cite{BL05}), which means that the scheme is simple to implement and very well-suited to the problem. 

The rest of this paper is organized as follows: in Section~\ref{Sec:model}, a two-dimensional model problem in mixed form %and its variational formulation 
is introduced. In Section~\ref{Sec:discrete}, we present the discretizations using FECC scheme (in which a technique of dual mesh is employed) and formulate the corresponding discrete problem. Some theoretical results concerning the stability of the new scheme are proved, and the advantages of using FECC scheme in terms of accuracy and computational cost are also discussed. Numerical experiments for two-dimensional problems comparing the performance of different schemes are shown in Section~\ref{Sec:NumRe}.  

% ###################################################
%
%			A MODEL PROBLEM
%
% ###################################################
\section{A model problem}
\label{Sec:model}
For a bounded domain $ \Omega $ in $ \mR^2 $ with Lipschitz boundary $ \partial \Omega $, we consider the following stationary linear elasticity equation:
\begin{equation} \label{primal}
-\Div \bsigma(\bu)=\bff \quad \text{in} \; \Omega,
\end{equation}
where $ \bu $ is the displacement of an elastic material, $ \bsigma $ the Cauchy stress and $ \bff $ the body forces. For simplicity, we impose a homogeneous Dirichlet boundary condition on $ \partial \Omega$:
\begin{equation} \label{BC}
\bu = 0 \quad \text{on} \; \partial \Omega.
\end{equation}
If the elastic material is isotropic and linear, the stress $ \bsigma (\bu) $ is defined by
\begin{equation*}
\bsigma (\bu)=2\mu \varepsilon (\bu) + \lambda \Div \bu \mathbf{Id},
\end{equation*}
where $ \mathbf{Id} $ is the identity matrix of size $ 2 $, $ \varepsilon (\bu) $ the infinitesimal strain tensor defined by
\begin{equation*}
\varepsilon (\bu) = \frac{1}{2}\left (\nabla \bu + (\nabla \bu)^{T}\right ),
\end{equation*}
(here $ A^{T} $ denotes the transpose of a matrix A),  $ \lambda $ and $ \mu $ are the Lam\'e constants defined by
\begin{equation*}
\lambda =\frac{\nu E}{(1+\nu)(1-2\nu)}, \quad \mu=\frac{E}{2(1+\nu)},
\end{equation*} 
with $ \nu $ the Poisson's ratio and $ E $ the Young's modulus. In this work, we study the case where $ \nu $ is close to $ 0.5 $ (or $\lambda$ is large), i.e. the material is nearly incompressible (for example rubber or rubber-like materials). It is well-known that for such a case, the standard finite elements might give inaccurate results due to volumetric locking and instability. In order to avoid this, one may use the mixed formulation by introducing an additional variable, the pressure $ p $, as follows:
\begin{equation*}
p:=\lambda \Div \bu \quad \text{in} \; \Omega,
\end{equation*} 
 and rewrite problem \eqref{primal}-\eqref{BC} equivalently in a mixed displacement-pressure form as
  \begin{subequations} \label{mixed}
 \begin{align}
 -\Div \left (2\mu \varepsilon (\bu) + p \mathbf{Id}\right ) &= \bff & \text{in} \; \Omega,  \label{1stmixed} \\
 \Div \bu - \frac{1}{\lambda} p & = 0 & \text{in} \; \Omega, \label{2ndmixed} \\
 \bu & =0 & \text{on} \; \partial \Omega.  \label{3rdmixed}
 \end{align}
\end{subequations}
% \begin{equation}\label{mixed}
% \begin{array}{rll}
% -\Div \left (2\mu \varepsilon (\bu) + p \mathbf{Id}\right ) &= \bff & \text{in} \; \Omega, \vspace{4pt}\\
% \Div \bu - \frac{1}{\lambda} p & = 0 & \text{in} \; \Omega, \vspace{4pt}\\
% \bu & =0 & \text{on} \; \partial \Omega. 
% \end{array}
%\end{equation}
To derive the variational formulation of \eqref{mixed}, we first introduce the following Sobolev spaces
\begin{equation*}
\bV_{0}=\left (H_{0}^{1}(\Omega)\right )^{2} \; \; \text{and} \; \; L_{0}^{2}(\Omega):=\left \{ q \in L^{2}(\Omega): \int_{\Omega} q d\Omega = 0\right \}.
\end{equation*}
Denote by $\| \cdot \|_{0}$ and $\| \cdot \|_{1}$ the norms defined on $L_{0}^{2}$ and $\bV_{0}$ respectively. We  shall seek for $ \bu \in \bV_{0} $ and $ p \in L_{0}^{2}(\Omega) $. The latter is obtained by integrating  equation \eqref{2ndmixed} over $ \Omega $ and using divergence theorem together with boundary condition \eqref{3rdmixed}. We also define the bilinear forms: 
\begin{equation*}
\begin{array}{lll}
a:& \bV_{0} \times \bV_{0} &\rightarrow \mR \\
& (\bu, \bv) & \mapsto a(\bu, \bv) = 2\mu \int_{\Omega} \varepsilon (\bu) : \varepsilon(\bv) \, d\Omega, \vspace{3pt}\\
%
b:& \bV_{0} \times L_{0}(\Omega) & \rightarrow \mR \\
& (\bu, q) & \mapsto b(\bu, q) = \int_{\Omega} q \Div \bv \, d\Omega, \vspace{3pt} \\
%
c:& L_{0}(\Omega) \times L_{0}(\Omega) & \rightarrow \mR \\
& (p, q) & \mapsto c(p, q) = \int_{\Omega}  pq \, d\Omega, \vspace{3pt}\\
%
L_{f}:&  \bV_{0} & \rightarrow \mR \\
& \bv & \mapsto L_{f}(\bv) = \int_{\Omega}  \bff \cdot \bv \, d\Omega,
\end{array}
\end{equation*}
where $ \bff \in \left (L^{2}(\Omega)\right )^{2} $. \vspace{4pt}

With these notations, the variational form of \eqref{mixed} is written as:  \vspace{-0.2cm}
\begin{eqnarray} \label{variational-mixed}
\mbox{Find $ \bu \in \bV_{0} $ and $ p \in L_{0}^{2}(\Omega) $ such that}
   \hspace{4cm}\nonumber\\    
\begin{array}{rll} a(\bu, \bv) + b(\bv, p)   & = L_{f}(\bv ), & \forall \bv \in \bV_{0}, \vspace{3pt}\\
 b(\bu, q) - \frac{1}{\lambda } c(p,q) & = 0, & \forall q \in L_{0}^{2}(\Omega).
\end{array}
\end{eqnarray}
%Next, we multiply equation \eqref{1stmixed} by a test function $ \bv \in \bV_{0} $, integrate over $ \Omega $ then apply Green's formula and boundary condition \eqref{3rdmixed} to obtain
%\begin{equation}
%2\mu \int_{\Omega} \varepsilon (\bu) : \varepsilon(\bv) + \int_{\Omega} p \Div \bv = \int_{\Omega} \bff \cdot \bv \quad \forall \bv \in \bV_{0}. 
%\end{equation}
In the next section, we present the finite element cell-centered (FECC) scheme \cite{PO12} to obtain a numerical solution to problem \eqref{variational-mixed}.
% ###################################################
%
%			THE DISCRETE PROBLEM
%
% ###################################################
\section{The discrete problem using FECC}
\label{Sec:discrete}
%Transitional paragraph
%\begin{itemize}
%\item The description for 2D problem: for simplicity, the extension to 3D can be easily followed (see \citep{Othesis12} for the case of diffusion problems). 
%\item The construction of 3 meshes $ \rightarrow $ explain ideas of FECC and its advantages (high accuracy but low cost). 
%\item The reason for choosing the approximate spaces for $\bu$ and $p$: the approximate space for $\bu$ is direct from FECC scheme while that for $p$ is chosen so that the ``macroelement condition" is satisfied $ \rightarrow $ macroelements in case of our scheme is a very natural approach. 
%\end{itemize}
We extend the idea of FECC scheme for stationary diffusion problems \citep{Othesis12} to the case of nearly incompressible elasticity problems. We recall that FECC scheme is based on a technique of dual mesh and it involves only cell unknowns. Three meshes are respectively constructed, the primal mesh, the dual mesh and the third mesh. The third mesh consisting of triangular elements can be seen as a refinement of the primal mesh. However, the number of unknowns does not increase compared with other cell-centered schemes. Thus, FECC scheme is effective in terms of accuracy and computational cost. Furthermore, by choosing appropriate approximate spaces for $\bu$ and $p$ we obtain a stable, low-order finite element scheme. This is done using macroelement techniques which turn out to be very natural in this case. 

In the following, we construct the meshes and introduce the approximate spaces to obtain the discrete problem associated with \eqref{variational-mixed}(cf. Subsection~\ref{Subsec:Discrete}). Then we prove that the discrete problem is well-posed and the scheme is stable (cf. Subsection~\ref{Subsec:proof}). Finally, we perform some calculations to obtain the algebra system corresponding to the discrete problem and discuss some issues concerning the efficiency of the scheme (cf. Subsection~\ref{Subsec:efficiency}).

%------------------------------------
%	The discretization
%------------------------------------
\subsection{The discretization} \label{Subsec:Discrete}
For a polygonal domain $ \Omega $, consider a triangulation $ \iT_{h} $ (of $\Omega$) that consists of non-empty connected close disjoint subsets of $\Omega$: $$ \overline{\Omega} = \bigcup_{K \in \iT_{h}} K.
$$
We assume that each element $K \in \iT_{h}$ is a star-shaped polygon in which we will choose a point $C_{K} \in \text{int}(K)$ and call it the mesh point of $K$. Throughout the paper, we refer to $\iT_{h}$ as the primal mesh. Next, we briefly recall the construction of the dual mesh $\iT_{h}^{*}$ and the third mesh $\iT_{h}^{**}$ which will be necessary to define the FECC scheme (see \citep{PO12, Othesis12} for a more detailed presentation). To define the dual mesh, we assume that the line joining two mesh points of any two neighboring elements is inside $\Omega$ but it doesn't need to intersect the common edge of the two element since the problem is homogeneous. 

The introduction of the dual mesh $\iT_{h}^{*}$ is based on the primal mesh so that each dual control volume  of $\iT_{h}^{*}$ corresponds to a vertex of $\iT_{h}$. Denote by $\iN$ the set of all nodes or vertices of $\iT_{h}$
$$ \iN : = \left \{ i: \; \text{i is a vertex of element $K \in \iT_{h}$ } \right \}.
$$
For each $i \in \iN$, denote by 
$$\iT_{i}:=\left \{ K \in \iT_{h}: \; \text{K shares the vertex i} \right \},$$ 
the set of primal elements that have $i$ as their vertex. We consider two cases (see Figure): 
\begin{enumerate}
\item[(a)] If $i$ is an interior vertex, then by connecting mesh points of neighboring elements in $\iK_{i}$ we obtain the dual control volume $M_{i}$ associated with the vertex $i$.
\item[(b)] If $i$ is on the boundary $\partial \Omega$, denote by $E_{i} \subset \partial K_{i}^{E}$ and $F_{i} \subset \partial K_{i}^{F}$ the two edges on the boundary that have $i$ as their vertex in which $K_{i}^{E}, K_{i}^{F} \in \iK_{i}$ can be two distinguished elements or just one element. 
The dual control volume $M_{i}$ is defined by joining mesh points of neighboring elements in $\iK_{i}$ and the mesh point of $K_{i}^{E}$ (and $K_{i}^{F}$) with the midpoint of $E_{i}$ (and $F_{i}$ respectively), note that in this case $M_{i}$ has $i$ as its vertex as well. 
\end{enumerate}
The collection of all $M_{i}$ defines a dual mesh $ \iT_{h}^{*} $ such that
$$ \overline{\Omega} = \bigcup_{i \in \iN} M_{i}.
$$
As for $ \iT_{h} $, we denote by $ C_{M} $ the mesh point of $ M \in \iT_{h}^{*}$. Note that if $ M_{i} $ has edges lying on $ \partial \Omega $, then $ C_{M_{i}} $ is chosen to be the corresponding vertex $i$ of the primal mesh (see Figure). Finally, we construct a third grid $ \iT_{h}^{**} $ as a triangular subgrid of the dual grid as follows: for an element $ M \in \iT_{h}^{*} $, we construct elements of $ \iT_{h}^{**} $ by connecting $ C_{M} $ to all vertices of $ \iT_{h}^{*} $ (see Figure):
$$ \overline{\Omega} = \bigcup_{T \in \iT_{h}^{**}} \overline{T}.
$$

\begin{remark} \label{rmk:macro}
By construction, each dual control volume $M \in \iT_{h}^{*}$ is indeed a macroelement - the union of a fixed number of adjacent elements of the third mesh. If one chooses $c_{M_{i}}$ to be the associated vertex $i$ (of the primal mesh), then $M_{i}$ consists of triangles $T \in \iT_{h}^{**}$ that have $i$ as their vertex. 
\end{remark}

By applying the FECC scheme, the finite element space for the displacement $ \bu $ is the standard finite elements of order 1 defined on the third mesh $ \iT_{h}^{**} $. For the pressure $p$, using macroelement techniques \cite{RS84}, \cite[pp.235-238] {BF91} and due to Remark~\ref{rmk:macro}, the finite element space for $p$ is chosen to be P0 functions on the dual mesh. Thus, $p_{h}$ is piecewise constant on each macroelement $M \in \iT_{h}^{*}$. The effect of this choice will be show in the next section when we prove the stability of the resulting scheme. In particular:
\begin{enumerate}
\item[(a)] For the pressure:
$$ p_{h} \in Q_{h}=\left \{ q_{h} \in L_{0}^{2}(\Omega): q_{h}\vert_{M} \in \mathbb{P}^{0}(M), \; \forall M \in \iT_{h}^{*} \right \} \subset L_{0}^{2}(\Omega).
$$
Thus
\begin{equation} \label{FEp}
 p_{h}(\bx) = \sum_{M \in \iT_{h}^{*}} p_{M} \chi_{M}(\bx),
\end{equation}
where $ \chi_{M} $ is the characteristic function of $ M \in \iT_{h}^{*} $. 
%%
\item[(b)] For the displacement: 
$$ \bu_{h} \in \bV_{h}=\left \{ \bv_{h} \in \bV_{0}: \bv_{h} \vert_{T} \in \left (\mathbb{P}^{1}(T)\right )^{2}, \; \forall T \in \iT_{h}^{**}  \right \} \subset \bV_{0}.
$$
The basis functions of $ \bV_{h} $ are defined at the nodes of elements of the third mesh. Since we impose homogeneous Dirichlet boundary conditions, we only need to deal with interior nodes. Denote by $\iN^{**}$ the set of interior nodes of elements of $\iT_{h}^{**}$.
%Denote by $\iC$ the set of mesh points of the primal mesh, $\iC^{*}$ the set of mesh points of interior dual control volumes of $\iT_{h}^{*}$ and denote by $\iN^{**}$ the set of interior nodes of elements of $\iT_{h}^{**}$. 
\begin{remark} \label{rmk:nodes}
 By construction, $\iN^{**}$ consists of mesh points of the primal mesh and mesh points of interior dual control volumes:
 $$\iN^{**} = \bigcup_{K \in \iT_{h}} C_{K} \cup \bigcup_{\substack{M \in \iT_{h}^{*} \\ \partial M \cap \partial \Omega = \emptyset}} C_{M}.$$
\end{remark}
Denote by $N_{P}$ the basis function of $\bV_{h}$ at node $P \in \iN^{**}$, we seek for $ u_{h} \in V_{h} $ of the form:
\begin{align} 
\bu_{h}(\bx)&= \sum_{P \in \iN^{**}} \left (u_{P}^{(1)}N_{P}(\bx), u_{P}^{(2)}N_{P}(\bx)\right ), \label{FEu} %\notag \\ 
%&=\sum_{P \in \iC} \left (u_{P}^{(1)}N_{P}(\bx), u_{P}^{(2)}N_{P}(\bx)\right )+\sum_{P^{*} \in \iC^{*}} \left (u_{P^{*}}^{(1)}N_{P^{*}}(\bx), u_{P^{*}}^{(2)} N_{P^{*}}(\bx)\right ), \label{FEu}
\end{align}
where $ \bu_{P} = \left (u_{P}^{(1)}, u_{P}^{(2)}\right ) $ is the nodal values of $ \bu_{h} $ at the vertex $ P \in \iN^{**}$. 
%The basis functions of $ \bV_{h} $, denoted by $ N_{T_{(l)}} \; \left (l=1,2,3\right )$, are defined at each node $ T_{(i)} $ of elements $ T \in \iT_{h}^{**} $. Hence, we seek for $ u_{h} \in V_{h} $ in the form:
%\begin{equation} \label{FEu}
%\bu_{h}(\bx)=\sum_{l=1}^{3} \, \text{diag}\, \left (N_{T_{(l)}}(\bx), \hdots,  N_{T_{(l)}}(\bx)\right ) \bu_{T_{(i)}},
%\end{equation}
%where $ \bu_{T_{(i)}}  $ is the nodal values of $ \bu_{h} $ at the vertex $ T_{(i)} $. 
\end{enumerate}
%
The discrete variational formulation of problem \eqref{variational-mixed} is then 
\begin{eqnarray} \label{discreteVM}
\mbox{Find $ \bu_{h} \in \bV_{h} $ and $ p_{h} \in Q_{h} $ such that}
   \hspace{4cm}\nonumber\\    
\begin{array}{rll} a(\bu_{h}, \bv_{h}) + b(\bv_{h}, p_{h})   & = L_{f}(\bv_{h} ), & \forall \bv_{h} \in \bV_{h}, \vspace{3pt}\\
 b(\bu_{h}, q_{h}) - \frac{1}{\lambda } c(p_{h},q_{h}) & = 0, & \forall q_{h} \in Q_{h}.
\end{array}
\end{eqnarray}
%
%
%
\subsection{Well-posedness of the discrete problem} \label{Subsec:proof}
According to the theory of mixed finite elements \cite[Chapter II, \S 2]{BF91}, the well-posedness of \eqref{discreteVM} is given by the following three conditions 
\begin{enumerate}
\item The bilinear form $a(\cdot, \cdot)$ is continuous, symmetric on $\bV_{h} \times \bV_{h}$ and is uniformly coercive on $\bV_{h}^{0}:=\left \{ \bv_{h} \in \bV_{h}: b(\bv_{h}, q_{h}) = 0, \forall q_{h} \in Q_{h} \right \}$, i.e. there exists $\alpha >0$ independent of the mesh size $h$ such that 
$$ a(\bv_{h}, \bv_{h}) \geq \alpha_{0} \| \bv_{h} \|_{1}, \; \forall \bv_{h} \in \bV_{h}^{0}. 
$$
%
\item The bilinear form $b(\cdot, \cdot)$ is continuous on $\bV_{h} \times Q_{h}$ and satisfies the uniform inf-sup condition (or Babu\v{s}ka-Brezzi stability condition), i.e. there exists $\beta >0$ independent of $h$ such that 
$$ \inf_{q_{h} \in Q_{h}} \, \sup_{\bv_{h} \in \bV_{h}} \frac{b(\bv_{h}, q_{h})}{\|\bv_{h}\|_{1} \| q_{h}\|_{0}} \geq \beta, \; \| \bv_{h}\|_{1} \neq 0, \, \| q_{h} \|_{0} \neq 0. 
$$
%
\item The bilinear form $c(\cdot, \cdot)$ is continuous, symmetric on $Q_{h} \times Q_{h}$ and is positive semi-definite:
$$  c(q_{h}, q_{h}) \geq 0, \; \forall q_{h} \in Q_{h}. 
$$ 
\end{enumerate}
We now check these three conditions: the continuity of the bilinear forms $a(\cdot, \cdot)$, $b(\cdot, \cdot)$ and $c(\cdot, \cdot)$ on their associated spaces is straightforward and so the symmetry of $a(\cdot, \cdot)$ and $c(\cdot, \cdot)$. It is simple to verify that $c(\cdot, \cdot)$ is positive semi-definite. The uniform coercivity of $a(\cdot, \cdot)$ is given by Korn's first inequality inequality \cite[Chapter III.3]{DL76}. There only remains to show that the $b(\cdot, \cdot)$ satisfies the uniform inf-sup condition. This is usually the main problem one has to deal with in order to prove the stability of the numerical scheme. However, in our case, this can be obtained directly by using the macroelement condition \cite{RS84}. In particular, for a macroelement $M \in \iT_{h}^{*}$, define the spaces
$$ \bV_{0,M}:= \left \{ \bv_{h} \in \bV_{h}: \bv_{h} = 0 \; \text{in} \; \Omega \setminus M \right \},
$$
and 
$$ N_{M}:=\left \{ q_{M}: q_{M} = q_{h} \vert_{M}, q_{h} \in Q_{h}, \int_{M} q_{h} \Div \bv_{h} d\bx = 0, \; \forall \bv_{h} \in \bV_{0,M}\right  \}. 
$$
By the definitions of $\bV_{h}$ and $Q_{h}$, we have that $N_{M}$ is one-dimensional (indeed, $q_{M}$ is piecewise constant on $M$ and $q_{M} \vert_{T_{1}} = q_{M} \vert_{T_{2}}, \; \forall T_{1}, T_{2} \subset M$).  
%In particular, assume that there is a fixed set of reference macroelements $\hat{M}_{j}, \, j=\overline{1,n}$ and denote by $\mathcal{E}_{\hat{M}_{j}}$ the family of macroelements equivalent with $\hat{M}_{j}$. The following three conditions can be verified straightforwardly:
%\begin{enumerate}
%\item[(i)] For each $M \in \mathcal{E}_{M_{j}}, \;j=\overline{1,n} $, the space 
%\end{enumerate}
  
Using classical results of the approximation of the saddle point problem, we have the following theorems:
\begin{theorem} \label{thrm:well-posed}
There exists a unique solution $(\bu_{h}, p_{h}) \in \bV_{h} \times Q_{h}$ of the discrete problem \eqref{discreteVM}. 
%The solution is uniformly stable with respect to $\bff$ in the sense that 
%$$ \| \bu_{h}\|_{1} + \| p_{h} \|_{0} \leq C \| \bff \|_{0},
%$$
%where $C$ is a constant independent of the mesh size. 
\end{theorem}
%
%
\begin{theorem}\label{thrm:estimates}
Let $(\bu, p) \in \bV_{0} \times L_{0}^{2}(\Omega)$ and $(\bu_{h}, p_{h}) \in \bV_{h} \times Q_{h}$ be the solutions to problems \eqref{variational-mixed} and \eqref{discreteVM} respectively. Then the following estimate, which is uniform with respect to $\lambda$, hold:
$$ \| \bu - \bu_{h} \|_{1} + \| p - p_{h} \|_{0} \leq C \left ( \inf_{\bv_{h} \in \bV_{h}} \| \bu - \bv_{h} \|_{1} + \inf_{q_{h} \in Q_{h}} \| p- q_{h} \|_{0} \right ),
$$
where $C$ is a constant independent of the mesh size. 
\end{theorem}  
In the rest of this section, we will present in detail the calculation of the scheme \eqref{discreteVM} for implementation purpose. It will be shown that the scheme is cell-centered in the sense that the unknowns for $\bu$ represent the average values of $\bu$ over primal elements $K \in \iT_{h}$ (after performing some linear transformation) while the unknowns for $p$ represent the average values of $p$ over dual control volumes $M \in \iT_{h}^{*}$. Thus, the number of unknowns for $\bu$ is the number of elements of the primal mesh and that for $p$ is the number of nodes of elements of the primal mesh. However, the scheme is more effective than other cell-centered scheme (using the same order of the approximation space) because we have approximated $\bu$ on a finer mesh (the third mesh) without requiring additional computational cost.
%
%
% 
\subsection{Accuracy and computational cost of FECC scheme}\label{Subsec:efficiency}
In order to simplify the calculation process, we firstly rewrite \eqref{discreteVM} equivalently as \cite[Chapter VI, p.201]{BF91}:
\begin{eqnarray} \label{discreteVM2}
\mbox{Find $ \bu_{h} \in \bV_{h} $ and $ p_{h} \in Q_{h} $ such that}
   \hspace{4cm}\nonumber\\    
\begin{array}{rll} \mu \int_{\Omega} \nabla \bu_{h} : \nabla \bv_{h} \, d \bx + (\mu +1) \int_{\Omega} p_{h} \Div \bv_{h}  \, d\bx  & = \int_{\Omega} \bff \cdot \bv_{h} \, d \bx, & \forall \bv_{h} \in \bV_{h}, \vspace{7pt}\\ 
 \int_{\Omega} q_{h} \Div \bu_{h} \,  d \bx - \frac{1}{\lambda } \int_{\Omega} p_{h} \, q_{h} \, d \bx & = 0, & \forall q_{h} \in Q_{h}.
\end{array}
\end{eqnarray}
%
In the following, we make use of the notation
$$ \partial_{x_{1}} g:= \frac{\partial g}{\partial x_{1}}, \; \partial_{x_{2}} g:= \frac{\partial g}{\partial x_{2}}, 
$$
for some function $g \in H^{1}(\Omega)$. 

Now to obtain the linear system from \eqref{discreteVM2}, we proceed as usual by choosing test functions of the first and the second equations of \eqref{discreteVM} as basis functions of $\bV_{h}$ and $Q_{h}$ respectively. 
%For a vector-valued function $\bv \in \left (L^{2}(\Omega)\right )^{2}$ $(d=2)$, we write $ \bv = (v^{(1)}, v^{(2)})$ where $v^{(1)}, v^{(2)} \in L^{2}(\Omega)$. If $\bv \in \left (H^{1}(\Omega)\right )^{2}$ then we denote by 
%$$ \nabla v^{j} = \left (\frac{\partial v^{j}}{\partial x_{1}}, \frac{\partial v^{j}}{\partial x_{2}}\right ), \; j=1,2.
%$$
%Finally we denote by $\be_{1} = (1,0)$ and $\be_{2} = (0,1)$ the unit vectors in $\mR^{(2)}$. 

\begin{itemize}
\item[(a)] From the representations of $\bu_{h}$ and $p_{h}$ in \eqref{FEu} and \eqref{FEp} respectively, letting $\bv_{h} = \left (N_{Q}(\bx), 0\right ), \; Q \in \iN^{**}$ in the first equation of \eqref{discreteVM2} to obtain
\begin{multline} \label{eq2:u1}
\mu \int_{\Omega} \sum_{P \in \iN^{**}}  \left ( u_{P}^{(1)} \partial_{x_{1}} N_{P}(\bx) \; \partial_{x_{1}} N_{Q}(\bx) + u_{P}^{(1)} \partial_{x_{2}} N_{P}(\bx) \; \partial_{x_{2}} N_{Q}(\bx) \right ) d\bx + \\
 (\mu +1) \int_{\Omega} \partial_{x_{1}} N_{Q}(\bx) \; \sum_{M \in \iT_{h}^{*}} p_{M} \chi_{M}(\bx)d\bx = \int_{\Omega} f^{(1)} \, N_{Q}(\bx) d\bx, \; \;  \forall Q \in \iN^{**}.
\end{multline}
Denote by $\iN_{M}$ the set of nodes of elements $T \in \iT_{h}^{**}$ that $T \subset M$. 
In the following, we will show how to condense out the unknowns $u_{C_{M}}^{(1)}, \; \forall M~\in~\iT_{h}^{*}, \; \partial M \cap \partial \Omega = \emptyset$ (i.e. interior dual control volumes) from the formulation. In fact, $u_{C_{M}}^{(1)}$ can be computed by a linear combination of unknowns at primal mesh points $u_{C_{K}}^{(1)}, \, K \in \iT_{h} $. \vspace{4pt}\\ To this purpose, let $\bv_{h} = \left (N_{C_{M}}, 0 \right ), \; M \in \iT_{h}^{*}, \; \partial M \cap \partial \Omega = \emptyset$ then \eqref{eq2:u1} becomes
\begin{multline} \label{eq2:u1bis}
\mu \int_{M} \sum_{P \in \iN_{M}}  \left ( u_{P}^{(1)} \partial_{x_{1}} N_{P}(\bx) \; \partial_{x_{1}} N_{C_{M}}(\bx) +  u_{P}^{(1)} \partial_{x_{2}} N_{P}(\bx) \; \partial_{x_{2}} N_{C_{M}}(\bx) \right ) + \\
(\mu +1) \int_{M} \partial_{x_{1}} N_{C_{M}}(\bx) \;  p_{M} = \int_{M} f^{(1)} \, N_{C_{M}}(\bx) d\bx.
\end{multline}
According to Remark~\ref{rmk:nodes}, $\iN_{M}$ consists of the mesh point $C_{M}$ of $M$ and the set of mesh points $C_{K}$ of elements $K \in \iT_{h}$ such that $K~\cap~M~\neq~\emptyset$.  %such that $C_{K}$ is a vertex of $M$. 
Thus, from \eqref{eq2:u1bis} we can compute $u_{C_{M}}^{(1)}$ using unknowns at primal mesh points of elements that intersect $M$:
%which is equal to $\Pi_{M}^{(1)} \left ( \left \{ u_{C_{K}}\right \}_{K \in \iT_{h}, \, K \cap M \neq \emptyset}, \bff \right )$ defined by 
\begin{multline} \label{eq2:uM1}
u_{C_{M}}^{(1)} = \Pi_{M}^{(1)} \left ( \left \{ u_{C_{K}}\right \}_{K \in \iT_{h}, \, K \cap M \neq \emptyset}, \bff \right )  =\\
 -\frac{1}{\Theta_{M}} \int_{M} \sum_{\substack{K \in \iK_{h} \\ K \cap M \neq \emptyset}}  \left ( u_{C_{K}}^{(1)} \partial_{x_{1}} N_{C_{K}}(\bx) \; \partial_{x_{1}} N_{C_{M}}(\bx) +  u_{C_{K}}^{(1)} \partial_{x_{2}} N_{C_{K}}(\bx) \; \partial_{x_{2}} N_{C_{M}}(\bx) \right ) + \\
 \frac{1}{\mu \, \Theta_{M}} \left ((\mu + 1) \int_{M} \partial_{x_{1}} N_{C_{M}}(\bx) \;  p_{M} - \int_{M} f^{(1)} \, N_{C_{M}}(\bx) d\bx \right ),
\end{multline}
where 
$$ \Theta_{M}:= \int_{M} \left (\partial_{x_{1}} N_{C_{M}}\right )^{2} + \left (\partial_{x_{2}} N_{C_{M}}\right )^{2}  d \bx. 
$$
For convenience, we write  $\Pi_{M}^{(1)}$ for 
 $\Pi_{M}^{(1)} \left ( \left \{ u_{C_{K}}\right \}_{K \in \iT_{h}, \, K \cap M \neq \emptyset}, \bff \right )$.
 %noting that this projection is only dependent on unknowns at primal mesh points. 
%
Now substituting \eqref{eq2:uM1} into \eqref{eq2:u1} we obtain
\begin{multline} \label{eq2:u1Final}
\mu \int_{\Omega} \sum_{K \in \iT_{h}}  \left ( u_{C_{K}}^{(1)} \partial_{x_{1}} N_{C_{K}}(\bx) \; \partial_{x_{1}} N_{C_{H}}(\bx) + u_{C_{K}}^{(1)} \partial_{x_{2}} N_{C_{K}}(\bx) \; \partial_{x_{2}} N_{C_{H}}(\bx) \right ) d\bx + \\
\mu \int_{\Omega} \sum_{\substack{M \in \iT_{h}^{*} \\ \partial M \cap \partial \Omega = \emptyset}}  \left ( \Pi_{M}^{(1)} \partial_{x_{1}} N_{C_{M}}(\bx) \; \partial_{x_{1}} N_{C_{H}}(\bx) + \Pi_{M}^{(1)} \partial_{x_{2}} N_{C_{M}}(\bx) \; \partial_{x_{2}} N_{C_{H}}(\bx) \right ) d\bx + \\
(\mu + 1)  \int_{\Omega} \partial_{x_{1}} N_{C_{H}}(\bx) \; \sum_{M \in \iT_{h}^{*}} p_{M} \chi_{M}(\bx)d\bx = \int_{\Omega} f^{(1)} \, N_{C_{H}}(\bx) d\bx, \; \;  \forall H \in \iT_{h}.
\end{multline}
This linear equation only involves unknowns of the primal elements. 
\item[(b)] Similar to (a), we now choose the test function $\bv_{h} = \left (0, N_{Q}(\bx) \right ), \; Q \in \iN^{**}$ in the first equation of \eqref{discreteVM} and obtain
\begin{multline} \label{eq2:u2Final}
\mu \int_{\Omega} \sum_{K \in \iT_{h}}  \left ( u_{C_{K}}^{(2)} \partial_{x_{2}} N_{C_{K}}(\bx) \; \partial_{x_{2}} N_{C_{H}}(\bx) + u_{C_{K}}^{(2)} \partial_{x_{1}} N_{C_{K}}(\bx) \; \partial_{x_{1}} N_{C_{H}}(\bx) \right ) d\bx + \\
\mu \int_{\Omega} \sum_{\substack{M \in \iT_{h}^{*} \\ \partial M \cap \partial \Omega = \emptyset}}  \left ( \Pi_{M}^{(2)} \partial_{x_{2}} N_{C_{M}}(\bx) \; \partial_{x_{2}} N_{C_{H}}(\bx) + \Pi_{M}^{(2)} \partial_{x_{1}} N_{C_{M}}(\bx) \; \partial_{x_{1}} N_{C_{H}}(\bx) \right ) d\bx + \\
 \int_{\Omega} \partial_{x_{2}} N_{C_{H}}(\bx) \; \sum_{M \in \iT_{h}^{*}} p_{M} \chi_{M}(\bx)d\bx = \int_{\Omega} f^{(2)} \, N_{C_{H}}(\bx) d\bx, \; \;  \forall H \in \iT_{h},
\end{multline}
where 
\begin{multline*}
 \Pi_{M}^{(2)} = \Pi_{M}^{(2)} \left ( \left \{ u_{C_{K}}\right \}_{K \in \iT_{h}, \, K \cap M \neq \emptyset}, \bff \right ) =\\
 -\frac{1}{\Theta_{M}} \int_{M} \sum_{\substack{K \in \iK_{h} \\ K \cap M \neq \emptyset}}  \left ( u_{C_{K}}^{(2)} \partial_{x_{1}} N_{C_{K}}(\bx) \; \partial_{x_{1}} N_{C_{M}}(\bx) +  u_{C_{K}}^{(2)} \partial_{x_{2}} N_{C_{K}}(\bx) \; \partial_{x_{2}} N_{C_{M}}(\bx) \right ) + \\
 \frac{1}{\mu \, \Theta_{M}} \left ((\mu + 1) \int_{M} \partial_{x_{2}} N_{C_{M}}(\bx) \;  p_{M} - \int_{M} f^{(2)} \, N_{C_{M}}(\bx) d\bx \right ).
\end{multline*}
\item[(c)] Letting $q_{h} = q_{M}, \, M \in \iT_{h}^{*}$, we rewrite the second equation of \eqref{discreteVM} as: 
\begin{equation*}
\sum_{P \in \iN_{M}} \int_{M} u_{P}^{(1)} \partial_{x_{1}} N_{P}(\bx) + u_{P}^{(2)} \partial_{x_{2}} N_{P}(\bx) d \bx -\frac{1}{\lambda} \int_{M} p_{M} d\bx = 0, \; \forall M \in \iT_{h}^{*}.
\end{equation*}
Applying the results in $(a)$ and $(b)$ we have
\begin{multline} \label{eq1:pFinal}
\sum_{\substack{K \in \iT_{h} \\ K \cap M \neq \emptyset}} \int_{M} u_{C_{K}}^{(1)} \partial_{x_{1}} N_{C_{K}}(\bx) + u_{C_{K}}^{(2)} \partial_{x_{2}} N_{C_{K}}(\bx) d \bx  + \\
\int_{M} \Pi_{M}^{(1)} \partial_{x_{1}} N_{C_{M}}(\bx) + \Pi_{M}^{(2)} \partial_{x_{2}} N_{C_{M}}(\bx) d \bx -\frac{1}{\lambda} \int_{M} p_{M} d\bx = 0, \; \forall M \in \iT_{h}^{*}.
\end{multline}
\end{itemize}
Finally, we end up with a linear system of problem \eqref{discreteVM} as follows:
\begin{equation}
\left [ \begin{array}{ll} \bA & \bB^{t} \\
\bB & -\frac{1}{\lambda} \bC
\end{array} \right ] \, \left [ \begin{array}{c} \bU \\
\bP
\end{array} \right ] =  \left [ \begin{array}{c} \bF \\
0
\end{array} \right ],
\end{equation}
where $\bU$ is a vector $\left ( u_{C_{T}}^{(1)}, u_{C_{T}}^{(2)} \right )_{T \in \iT_{h}}$ and $\bP$ is $\left (p_{M}\right )_{M \in \iT_{h}^{*}}$, which implies that the scheme is cell-centered. 


% ###################################################
%
%			NUMERICAL RESUTLS
%
% ###################################################
\section{Numerical experiments}
\label{Sec:NumRe}
We carry out several test cases in 2D to verify the performance of the constructed scheme and to compare its performance with other schemes for nearly incompressible elasticity problems. 

%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix
\section*{References}
%\label{}

%% If you have bibdatabase file and want bibtex to generate the
%% bibitems, please use
%%
%  \bibliographystyle{siam} %{elsarticle-num} 
%  \bibliography{bFECC}

%% else use the following coding to input the bibitems directly in the
%% TeX file.

\begin{thebibliography}{00}

\bibitem{BL05} {Lamichhane, Bishnu P.}, {\em Inf-sup stable finite-element pairs based on dual meshes and bases for nearly incompressible elasticity}, IMA Journal of Numerical Analysis, Volume 29, 2, pp. 404--420, (2009).
\bibitem{PO12} {Le Potier, C. and Ong, Thanh Hai.},	{International Journal on Finite Volumes}, {\em A cell-centered scheme for heterogeneous anisotropic diffusion problems on general meshes}, {8}, pp. 1-40, (2012).

\bibitem{RP87} {Pierre, R.}, {\em Regularization procedures of mixed finite element approximations of the {S}tokes problem}, {INRIA research report, N$^o $~673},(1987)

\bibitem{JM05} {Matsumoto, J.}, {\em A relationship between stabilization {FEM} and bubble function element stabilization method with orthogonal basis for incompressible flows}, {Journal of Applied Mechanics}, volume {8}, {(2005)}.

\bibitem{RS84} {Stenberg, Rolf}, {\em Analysis of mixed finite elements methods for the {S}tokes problem: a unified approach}. {Mathematics of Computation}, volume {42}, {165},
pp. {9--23}, {(1984)}.

\bibitem{Othesis12} {Ong, Thanh Hai.}, {\em Cell-centered scheme for heterogeneous anisotropic diffusion problems on general meshes}, {PhD thesis, Universit\'e Paris-Est},
 {(2012)}

\bibitem{BF91} {Brezzi, F and Fortin, M.}, {\em Mixed and hybrid finite elements methods}, {Springer-Verlag}, {(1991)}. 

\bibitem{DL76} {Duvaut, G and Lions, J.~L.}, {\em Inequalities in mechanics and physics},
 {Springer-Verlag}, {Berlin New York}, (1976).

\end{thebibliography}



\end{document}
\endinput
%%
%% End of file `elsarticle-template-num.tex'.

1406.2073.extracted_bib
{BL05 {Lamichhane, Bishnu P., {\em Inf-sup stable finite-element pairs based on dual meshes and bases for nearly incompressible elasticity, IMA Journal of Numerical Analysis, Volume 29, 2, pp. 404--420, (2009).
1406.2073.extracted_bib
{PO12 {Le Potier, C. and Ong, Thanh Hai., {International Journal on Finite Volumes, {\em A cell-centered scheme for heterogeneous anisotropic diffusion problems on general meshes, {8, pp. 1-40, (2012).
1406.2073.extracted_bib
{RP87 {Pierre, R., {\em Regularization procedures of mixed finite element approximations of the {Stokes problem, {INRIA research report, N$^o $~673,(1987)
1406.2073.extracted_bib
{JM05 {Matsumoto, J., {\em A relationship between stabilization {FEM and bubble function element stabilization method with orthogonal basis for incompressible flows, {Journal of Applied Mechanics, volume {8, {(2005).
1406.2073.extracted_bib
{RS84 {Stenberg, Rolf, {\em Analysis of mixed finite elements methods for the {Stokes problem: a unified approach. {Mathematics of Computation, volume {42, {165, pp. {9--23, {(1984).
1406.2073.extracted_bib
{Othesis12 {Ong, Thanh Hai., {\em Cell-centered scheme for heterogeneous anisotropic diffusion problems on general meshes, {PhD thesis, Universit\'e Paris-Est, {(2012)
1406.2073.extracted_bib
{BF91 {Brezzi, F and Fortin, M., {\em Mixed and hybrid finite elements methods, {Springer-Verlag, {(1991).
1406.2073.extracted_bib
{DL76 {Duvaut, G and Lions, J.~L., {\em Inequalities in mechanics and physics, {Springer-Verlag, {Berlin New York, (1976).