#  Expandable Local and Parallel Two-Grid Finite Element Scheme for the Stokes
Equations

Yanren Hou School of Mathematics and Statistics, Xi’an Jiaotong University,
Xi’an, Shaanxi 710049, China. (yrhou@mail.xjtu.edu.cn). Partially supported by
NSFC (Grant No. 11971378 & 11571274). Feng Shi College of Science, Harbin
Institute of Technology, Shenzhen 518055, China. (shi.feng@hit.edu.cn).
Partially supported by Foundation Research Project of Shenzhen (Grant No.
JCYJ2018030617181319 & ZDSYS201707280904031) Haibiao Zheng School of
Mathematical Sciences, East China Normal University, Shanghai Key Laboratory
of Pure Mathematics and Mathematical Practice, Shanghai, P.R. China.
(hbzheng13@gmail.com). Partially supported by NSFC (Grant No. 11971174) and
NSF of Shanghai (Grant No. 19ZR1414300) and Science and Technology Commission
of Shanghai Municipality (Grant No. 18dz2271000 & 19JC1420102).

###### Abstract

In this paper, we present a novel local and parallel two-grid finite element
scheme for solving the Stokes equations, and rigorously establish its a priori
error estimates. The scheme admits simultaneously small scales of subproblems
and distances between subdomains and its expansions, and hence can be
expandable. Based on the a priori error estimates, we provide a corresponding
iterative scheme with suitable iteration number. The resulting iterative
scheme can reach the optimal convergence orders within specific two-grid
iterations ($O(|\ln H|^{2})$ in 2-D and $O(|\ln H|)$ in 3-D) if the coarse
mesh size $H$ and the fine mesh size $h$ are properly chosen. Finally, some
numerical tests including 2-D and 3-D cases are carried out to verify our
theoretical results. ’

###### keywords:

two-grid method, domain decomposition method, superposition principle, local
and parallel, iterative scheme

###### AMS:

65N15, 65N30, 65N55

##  1 Introduction

Due to the limiting of computer resources, two-grid finite element
methods/nonlinear Galerkin schemes [1, 2] and domain decomposition methods [3]
are popular and powerful tools for numerical simulations of linear and
nonlinear PDEs nowadays. Such as two-grid/two-level post-processing schemes
for incompressible flow, we refer [4, 5, 6, 7, 8, 9] and the references
therein for details.

In the past decades, a local and parallel two-grid finite element method for
elliptic boundary value problems was initially proposed in [10]. The scheme
firstly solves the elliptic equation on a coarse mesh to get an initial lower
frequency guess of the solution. Then the whole computational domain is
divided into a series of disjoint subdomains $\\{D_{j}\\}$ and the driven term
of the error equation, namely the residual term, is split into several parts
defined only on such small subdomains. Finally the global error equation can
be transferred into a series of subproblems with local driven terms. Since the
solution of higher frequency to each subproblem decays very fast apart from
the support of the local driven term, by suitablely expanding each $D_{j}$ to
the domain $\Omega_{j}\supset D_{j}$ and imposing the homogeneous boundary
condition on $\partial\Omega_{j}$, each subproblem can be approximated in a
localized version defined in the corresponding expanded domain $\Omega_{j}$.
The most attractive feature of the scheme is that not any communication is
required between local fine grid subproblems, which makes the scheme a highly
effective parallel scheme. Such local and parallel two-grid scheme can be
found in [11] and has been extended to the Stokes equations in [12]. Error
estimates derived in [10, 12] show that the approximate solutions in such
schems can reach the optimal convergence orders in both $H^{1}$ and $L^{2}$
norms.

However, according to [13], the error estimates are limited by the usage of
the superapproximation property of finite element spaces, which makes the
error constant appeared in [10, 12] to be a form of $O(t^{-1})$, where $t$
denotes the distance between $\partial D_{j}$ and $\partial\Omega_{j}$. To
obtain the optimal error orders, usually $t=O(1)$ is required, which means the
distance of $\partial\Omega_{j}$ and $\partial D_{j}$ is almost a constant.
Therefore one can not expect $\Omega_{j}$ to be arbitrary small. This will
prevent the corresponding scheme from utilization in large parallel computer
systems. In the previously mentioned local and parallel schemes, computational
results for each $\Omega_{j}$ are usually removed outside $D_{j}$, and are
simply pasted together to form the final approximation, in which discontinuity
may appear along the boundaries of different $D_{j}$. Instead, based on the
method of partition of unity [14], a local and parallel two-grid scheme is
proposed for second order linear elliptic equations [15, 16] and is also
extended to dealing with the Stokes equations and Navier-Stokes equations[17,
18]. Although the partition of unity method can guarantee that the global
approximation is continuous, usage of superposition principle causes a crucial
requirement that the distance $t$ should be $O(1)$. To overcome such defects,
some research on linear elliptic equations has been done by the first author
and his collaborator in [19] by iterative method.

In this paper, based on the basic idea presented in [10, 12], as an important
extension of the idea in [19], we construct a local and parallel two-grid
iterative scheme for solving the Stokes equations, in which the scale of each
subproblem is $O(H)$, just requiring that $\mbox{diam}(D_{j})=O(H)$ and
$t=O(H)$ (much smaller than $O(1)$). Since only a small overlapping in each
two adjacent subproblems while $H$ tends to zero, the scale of each subproblem
can be arbitrary small as $O(H)$, that’s a main reason why we call the scheme
an expandable local and parallel two-grid scheme. Meanwhile, in each cycle of
two-grid iteration, to guarantee a better $L^{2}$ error estimate, we adopt a
coarse grid correction. Another main contribution in this paper is that, to
yield the globally continuous velocity in $\Omega$, we use the principle of
superposition based on a partition of unity to generate a series of local and
independent subproblems. In particular, for patches of given size, to obtain a
similar approximate accuracy as the one from the standard Galerkin method on
the fine mesh, we carry out rigorous analysis, and through the a priori error
estimate of the scheme, we show that a few number of iterations of order
$O(|\ln H|^{2})$ in 2-D and $O(|\ln H|)$ in 3-D is only needed, respectively.
Similar technique has been successfully applied for adaptive schemes with some
a posterior error estimates in [20, 21].

The remainder of this paper is organized as follows. In Section 2, we will
introduce the model problem and some preliminary materials. In Section 3, we
will present our expandable local and parallel scheme for the Stokes
equations. The a priori error estimates of the scheme are derived and then
suggested iterative scheme is presented in Section 4. Some numerical
experiments including 2-D and 3-D examples are carried out in a parallel
computer system to support our theoretical analysis in Section 5. Finally we
give some conclusions in Section 5 6.

##  2 Preliminaries

Consider the following Stokes equations in a smooth bounded convex domain
$\Omega\subset R^{d}$, $d=2,3$.

$\left\\{\begin{array}[]{rl}-\nu\Delta\mathbf{u}+\nabla p=\mathbf{f},&\mbox{
in }\Omega,\\\ \nabla\cdot\mathbf{u}=0,&\mbox{ in }\Omega,\\\
\mathbf{u}=0,&\mbox{ on }\partial\Omega.\end{array}\right.$ (1)

For a given bounded domain $D\subset R^{d}$, standard notations for Sobolev
spaces $W^{s,p}(D)$ and their associated norms will be used, see, e.g., [22,
23]. In particular, for $p=2$, we simply define $H^{s}(D)=W^{s,2}(D)$,
$H_{0}^{1}(D)=\\{v\in H^{1}(D):v|_{\partial D}=0\\}$, and their norms as
$\|\cdot\|_{s,D}=\|\cdot\|_{s,2,D}$ and $|\cdot|_{s,D}$. We denote by
$(\cdot,\cdot)_{D}$ the $L^{2}-$inner product on $D$. Therefore, we have
$\|\cdot\|_{0,D}=(\cdot,\cdot)_{D}^{\frac{1}{2}}$, and we use $\|\cdot\|_{D}$
to denote $\|\cdot\|_{0,D}$ in the rest of the paper.

Hereafter, we always use boldface characters to denote vector valued functions
or spaces, for instance, ${\boldsymbol{H}}^{s}(D)=(H^{s}(D))^{d}$. For
simplicity, we will also use the symbols $\lesssim$, $\gtrsim$ and $\approxeq$
in the following sense: $x_{1}\lesssim y_{1}$, $x_{2}\gtrsim y_{2}$ and
$x_{3}\approxeq y_{3}$ are equivalent to $x_{1}\leq C_{1}y_{1}$, $x_{2}\geq
c_{2}y_{2}$ and $c_{3}x_{3}\leq y_{3}\leq C_{3}x_{3}$ for some positive
constants $C_{1}$, $c_{2}$, $c_{3}$ and $C_{3}$ independent of mesh size. We
know that $\|\cdot\|_{1,D}\approxeq\|\nabla\cdot\|_{0,D}$ in $H^{1}_{0}(D)$.
For subdomains $D_{1}\subset D_{2}\subset\Omega$, we use $D_{1}\subset\subset
D_{2}$ to express that $\mbox{dist}(\partial D_{1},\partial
D_{2})\stackrel{{\scriptstyle\triangle}}{{=}}\mbox{dist}(\partial
D_{2}\backslash\partial\Omega,\partial D_{1}\backslash\partial\Omega)>0$.

Furthermore, let us denote

$\displaystyle{\boldsymbol{X}}(D)={\boldsymbol{H}}^{1}(D),\quad{\boldsymbol{X}}_{0}(D)={\boldsymbol{H}}_{0}^{1}(D),\quad
M(D)=L^{2}(D)$ $\displaystyle Q(D)=L_{0}^{2}(D)=\\{q\in
L^{2}(D):\;\int_{D}q=0\\},$

and we introduce the following bilinear forms

$\displaystyle
a_{D}({\boldsymbol{u}},{\boldsymbol{v}})=\nu(\nabla{\boldsymbol{u}},\nabla{\boldsymbol{v}})_{D},\;\;d_{D}(p,{\boldsymbol{v}})=-(p,\nabla\cdot{\boldsymbol{v}})_{D}\;\;\forall{\boldsymbol{u}},{\boldsymbol{v}}\in{\boldsymbol{X}}(D),\;p\in
M(D).$

By these notations, when $D=\Omega$, we simply denote by
${\boldsymbol{X}}={\boldsymbol{X}}(\Omega)$,
${\boldsymbol{X}}_{0}={\boldsymbol{X}}_{0}(\Omega)$, $M=M(\Omega)$,
$Q=Q(\Omega)$, and $a(\cdot,\cdot)=a_{\Omega}(\cdot,\cdot)$,
$d(\cdot,\cdot)=d_{\Omega}(\cdot,\cdot)$ and
$(\cdot,\cdot)=(\cdot,\cdot)_{\Omega}$. Then the weak formulation of (1) is:
find $[{\boldsymbol{u}},p]\in{\boldsymbol{X}}_{0}\times Q$ such that

$\left\\{\begin{array}[]{ll}a({\boldsymbol{u}},{\boldsymbol{v}})+d(p,{\boldsymbol{v}})=({\boldsymbol{f}},{\boldsymbol{v}})&\forall{\boldsymbol{v}}\in{\boldsymbol{X}}_{0},\\\
d(q,{\boldsymbol{u}})=0&\forall q\in Q.\end{array}\right.$

If we further introduce the bilinear form

$B_{D}([{\boldsymbol{u}},p],[{\boldsymbol{v}},q])=a_{D}({\boldsymbol{u}},{\boldsymbol{v}})+d_{D}(p,{\boldsymbol{v}})-d_{D}(q,{\boldsymbol{u}}),$

and $B(\cdot,\cdot)=B_{\Omega}(\cdot,\cdot)$, then the above weak form for the
Stokes equations can be rewritten as

$B([{\boldsymbol{u}},p],[{\boldsymbol{v}},q])=({\boldsymbol{f}},{\boldsymbol{v}})\quad\forall[{\boldsymbol{v}},q]\in{\boldsymbol{X}}_{0}\times
Q.$ (2)

Assume that $T^{H}(\Omega)=\\{\tau_{\Omega}^{H}\\}$ is a regular triangulation
of $\Omega$, with the mesh size defined as $H=\max_{\,\tau_{\Omega}^{H}\in
T^{H}(\Omega)}\\{\mbox{diam}(\tau_{\Omega}^{H})\\}$. For certain given
positive integer $r\geq 1$, let ${\boldsymbol{X}}^{H}\subset{\boldsymbol{X}}$
and $M^{H}\subset M$ be $C^{0}-$finite element spaces defined on $\Omega$ with
approximation order $r+1$ in ${\boldsymbol{H}}^{r+1}(\Omega)$ and
approximation order $r$ in $H^{r}(\Omega)$, respectively. We denote
${\boldsymbol{X}}_{0}^{H}={\boldsymbol{X}}^{H}\cap{\boldsymbol{H}}^{1}_{0}(\Omega)$,
$Q^{H}=M^{H}\cap L_{0}^{2}(\Omega)$ and assume ${\boldsymbol{X}}_{0}^{H}\times
Q^{H}$ is a stable finite element pair. Given $D\subset\Omega$, which aligns
with $T^{H}(\Omega)$, we denote $T^{H}(D)$ as the restriction of
$T^{H}(\Omega)$ on $D$, ${\boldsymbol{X}}^{H}(D)$ and $M^{H}(D)$ as the
restriction of ${\boldsymbol{X}}^{H}$ and $M^{H}$ on $D$ and
$Q^{H}(D)=M^{H}(D)\cap L_{0}^{2}(D)$.

Based on the finite element spaces defined above and weak formulation (2), we
make the following assumptions.

  * A1.

_Interpolant_. Let $I_{H}$ be a Lagrange finite element interpolation of
${\boldsymbol{X}}(D)$ onto ${\boldsymbol{X}}^{H}(D)$ and
$\hat{I}_{H}=I-I_{H}$. There holds for any
${\boldsymbol{w}}\in{\boldsymbol{H}}^{s}(\tau_{\Omega}^{H})$, $0\leq m\leq s$
and $s>d/2$,

$\|\hat{I}_{H}{\boldsymbol{w}}\|_{m,\tau_{\Omega}^{H}}\lesssim
H^{s-m}|{\boldsymbol{w}}|_{s,\tau_{\Omega}^{H}}.$

  * A2.

_Inverse Inequality_. For any ${\boldsymbol{w}}\in{\boldsymbol{X}}^{H}(D)$,
$0\leq m\leq s$, $1\leq p,q\leq\infty$,

$\|{\boldsymbol{w}}\|_{s,p,\tau_{\Omega}^{H}}\lesssim
H^{m-s+\frac{d}{p}-\frac{d}{q}}\|{\boldsymbol{w}}\|_{m,q,\tau_{\Omega}^{H}}.$

  * A3.

_Regularity_. Assume $D$ is smooth such that for any
${\boldsymbol{f}}\in{\boldsymbol{L}}^{2}(D)$, any solution to

$B_{D}([{\boldsymbol{u}},p],[{\boldsymbol{v}},q])=({\boldsymbol{f}},{\boldsymbol{v}})_{D}\quad\forall[{\boldsymbol{v}},q]\in{\boldsymbol{X}}_{0}(D)\times
Q(D),$

satisfies

$\|{\boldsymbol{u}}\|_{2,D}+\|p\|_{1,D}\lesssim\|{\boldsymbol{f}}\|_{D}.$

Now we can state the standard Galerkin approximation of (2) as follows: find
$[{\boldsymbol{u}}_{H},p_{H}]\in{\boldsymbol{X}}_{0}^{H}\times Q^{H}$ such
that

$B([{\boldsymbol{u}}_{H},p_{H}],[{\boldsymbol{v}}_{H},q_{H}])=({\boldsymbol{f}},{\boldsymbol{v}}_{H})\quad\forall\,[{\boldsymbol{v}}_{H},q_{H}]\in{\boldsymbol{X}}_{0}^{H}\times
Q^{H}.$ (3)

For this problem, a classical result holds, namely assuming that
$[{\boldsymbol{u}},p]\in({\boldsymbol{H}}^{r+1}(\Omega)\times
H^{r}(\Omega))\cap({\boldsymbol{X}}_{0}\times Q)$, one have

$\|{\boldsymbol{u}}-{\boldsymbol{u}}_{H}\|_{\Omega}+H(\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H})\|_{\Omega}+\|p-p_{H}\|_{\Omega})=O(H^{r+1}).$
(4)

##  3 Local and Parallel Two-Grid Scheme

In this section, we will introduce our local and parallel two-grid scheme for
solving Stokes equations. Firstly, we define the error functions by

$\hat{\boldsymbol{u}}={\boldsymbol{u}}-{\boldsymbol{u}}_{H}\in{\boldsymbol{X}}_{0},\quad\hat{p}=p-p_{H}\in
Q,$

where $[{\boldsymbol{u}},p]$ and $[{\boldsymbol{u}}_{H},p_{H}]$ are solutions
to (2) and (3), respectively. Then $[\hat{\boldsymbol{u}},\hat{p}]$ satisfies
the error equation

$B([\hat{\boldsymbol{u}},\hat{p}],[{\boldsymbol{v}},q])=({\boldsymbol{f}},{\boldsymbol{v}})-B([{\boldsymbol{u}}_{H},p_{H}],[{\boldsymbol{v}},q])\quad\forall[{\boldsymbol{v}},q]\in{\boldsymbol{X}}_{0}\times
Q.$ (5)

For any given partition of unity of $\Omega$, namely
$\\{\phi_{j}\\}_{j=1}^{N}$ with $N\geq 1$ an integer,
$\Omega\subset\bigcup\limits_{j=1}^{N}\mbox{supp }\phi_{j}$ and
$\sum\limits_{j=1}^{N}\phi_{j}\equiv 1$ on $\Omega$, we denote
$D_{j}=\mbox{supp}\,\phi_{j}$ and always assume that $D_{j}$ aligns with
$T^{H}(\Omega)$. With this patition of unity, the error equation (5) can be
rewritten as:

$B([\hat{\boldsymbol{u}},\hat{p}],[{\boldsymbol{v}},q])=({\boldsymbol{f}},\sum\limits_{j=1}^{N}\phi_{j}{\boldsymbol{v}})-B([{\boldsymbol{u}}_{H},p_{H}],\sum\limits_{j=1}^{N}\phi_{j}[{\boldsymbol{v}},q])\quad\forall[{\boldsymbol{v}},q]\in{\boldsymbol{X}}_{0}\times
Q.$

Clearly, by the principle of superposition, this equation is equivalent to the
summation of the following $N$-subproblems:

$B([\hat{\boldsymbol{u}}^{j},\hat{p}^{j}],[{\boldsymbol{v}},q])=({\boldsymbol{f}},\phi_{j}{\boldsymbol{v}})-B([{\boldsymbol{u}}_{H},p_{H}],\phi_{j}[{\boldsymbol{v}},q])\quad\forall[{\boldsymbol{v}},q]\in{\boldsymbol{X}}_{0}\times
Q.$ (6)

Namely,
$[\hat{\boldsymbol{u}},\hat{p}]=\sum\limits_{j=1}^{N}[\hat{\boldsymbol{u}}^{j},\hat{p}^{j}]$.
The most important feature for such subproblems is that all the subproblems
are mutually independent when $[{\boldsymbol{u}}_{H},p_{H}]$ is known.
Meanwhile each subproblem is globally defined with homogeneous Dirichlet
boundary condition for the velocity, however is driven by right-hand-side term
with very small compact support, namely $D_{j}$. As is pointed out in [21]
(used also in [19]) for elliptic problems, solution $\hat{\boldsymbol{u}}^{j}$
to each subproblem may decay very fast away from $D_{j}$. So for the present
Stokes problems, we will localize each subproblem to a small extension domain
of $D_{j}$, namely $\Omega_{j}$, and impose homogeneous Dirichlet boundary
condition for the velocity and discretize the subproblem on some triangulation
$T^{h}(\Omega_{j})$ for finite element approximation.

######  Remark 3.1.

For each $\Omega_{j}$, here we give a much natural choice. The partition of
unity functions $\\{\phi_{j}\\}_{j=1}^{N}$ of $\Omega$ is selected as the
piecewise linear Lagrange basis functions associated with the coarse
triangulation $T^{H}(\Omega)$, where $N$ is the number of all vertices in
$T^{H}(\Omega)$ including the boundary ones. For each vertex $x_{j}$ of
$T^{H}(\Omega)$, we can denote $D_{j}=\mbox{supp }\phi_{j}$, and construct
$\Omega_{j}$ by one layer extension of $D_{j}$ in $T^{H}(\Omega)$, namely
$\Omega_{j}=\bigcup_{x_{i}\in D_{j}}D_{i}.$ We refer [17] for more details.
Clearly, we know the scale esimates

$\mbox{diam}(D_{j}),\;\mbox{dist}(\partial D_{j},\partial\Omega_{j})\approxeq
H.$ (7)

Generally, we need assume that each $D_{j}$ and its extension domain
$\Omega_{j}$ are aligned with $T^{H}(\Omega)$. For each local domain
$\Omega_{j}$, its regular triangulation is denoted by
$T^{h}(\Omega_{j})=\\{\tau^{h}_{\Omega_{j}}\\}$, with mesh scale over all $j$
as $h=\max_{1\leq j\leq N}\left\\{\max_{\,\tau^{h}_{\Omega_{j}}\in
T^{h}(\Omega_{j})}\\{\mbox{diam}(\tau^{h}_{\Omega_{j}})\\}\right\\}$. As to
the natural choices of partition of unity functions and subdomains as
discussed in Remark 3.1, we can simply select local fine mesh
$T^{h}(\Omega_{j})$ as the limitation of the global regular triangulation of
$\Omega$, $T^{h}(\Omega)=\\{\tau^{h}_{\Omega}\\}$, namely,
$T^{h}(\Omega_{j})=T^{h}(\Omega)|_{\Omega_{j}}$, which will be used in the
rest of this paper.

Base on the fine mesh triangulations $T^{h}(\Omega_{j})$ and $T^{h}(\Omega)$
for each $\Omega_{j}$ and $\Omega$, we can introduce corresponding finite
element spaces ${\boldsymbol{X}}^{h}(\Omega_{j})$,
${\boldsymbol{X}}_{0}^{h}(\Omega_{j})$ and ${\boldsymbol{X}}^{h}$,
${\boldsymbol{X}}_{0}^{h}$, similarly defined as ${\boldsymbol{X}}^{H}$ and
${\boldsymbol{X}}_{0}^{H}$ previously. We also introduce $M^{h}(\Omega_{j})$,
$Q^{h}(\Omega_{j})$, $M^{h}$ and $Q^{h}$ for the pressure accordingly. Noting
that any function in ${\boldsymbol{X}}_{0}^{h}(\Omega_{j})$ can be extended to
a function in ${\boldsymbol{X}}_{0}^{h}$ through zero value assignment over
$\Omega\backslash\Omega_{j}$, we regard ${\boldsymbol{X}}_{0}^{h}(\Omega_{j})$
as a subspace of ${\boldsymbol{X}}_{0}^{h}$ in the sense of such zero
extension. In addition, we always assume that

$\left\\{\begin{array}[]{l}{\boldsymbol{X}}^{H}\subset{\boldsymbol{X}}^{h},\quad{\boldsymbol{X}}_{0}^{H}\subset{\boldsymbol{X}}_{0}^{h}=\bigcup\limits_{1\leq
j\leq N}{\boldsymbol{X}}_{0}^{h}(\Omega_{j}),\\\ M^{H}\subset M^{h},\quad
M^{H}|_{\Omega_{j}}\subset M^{h}(\Omega_{j}).\end{array}\right.$ (8)

Besides, we define discontinuous finite element spaces associated with the
pressure

$M_{j}^{h}=\\{q_{h}:\;q_{h}|_{\Omega_{j}}\in
M^{h}(\Omega_{j}),\;q_{h}|_{\Omega\backslash\Omega_{j}}\in
M^{h}(\Omega\backslash\Omega_{j})\\},\quad Q^{h}_{j}=M_{j}^{h}\cap
L_{0}^{2}(\Omega).$

Clearly, $M^{h}\subset\bigcup\limits_{j=1}^{N}M_{j}^{h}\subset
M_{H}^{h}=\\{q_{h}\in L^{2}(\Omega):q_{h}|_{\tau_{H}}\in M^{h}(\tau_{H})\\}$.
Then for $Q_{H}^{h}=M_{H}^{h}\cap L_{0}^{2}(\Omega)$, we have that
$Q^{h}\subset Q_{H}^{h}$ and the finite element pair
${\boldsymbol{X}}_{0}^{h}\times Q_{H}^{h}$ is stable (see [24]).

Right now, we can present the localized finite element approximation to error
equation (6) as follows: find
$[\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}]\in{\boldsymbol{X}}_{0}^{h}(\Omega_{j})\times
Q^{h}(\Omega_{j})$ such that
$\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}(\Omega_{j})\times
Q^{h}(\Omega_{j})$

$B_{\Omega_{j}}([\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},\phi_{j}{\boldsymbol{v}}_{h})_{\Omega_{j}}-B_{\Omega_{j}}([{\boldsymbol{u}}_{H},p_{H}],\phi_{j}[{\boldsymbol{v}}_{h},q_{h}]),$
(9)

in which $\hat{\boldsymbol{u}}^{j}_{H,h}$ and $\hat{p}^{j}_{H,h}$ are defined
on $\Omega_{j}$. In the rest, we always use the same symbols, that is
$[\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}]$, to denote their zero
extension over $\Omega\backslash\Omega_{j}$. In such sense,
$[\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}]\in{\boldsymbol{X}}_{0}^{h}\times
Q_{j}^{h}(\Omega)$. Then we construct totally

$\hat{\boldsymbol{u}}_{H,h}=\sum\limits_{j=1}^{N}\hat{\boldsymbol{u}}^{j}_{H,h}\in{\boldsymbol{X}}_{0}^{h},\quad\hat{p}_{H,h}=\sum\limits_{j=1}^{N}\hat{p}_{H,h}^{j}\in
Q_{H}^{h},$

and the intermediate approximate solution

$[{\boldsymbol{u}}_{H,h},p_{H,h}]=[{\boldsymbol{u}}_{H}+\hat{\boldsymbol{u}}_{H,h},p_{H}+\hat{p}_{H,h}].$
(10)

Since this approximation $[{\boldsymbol{u}}_{H,h},p_{H,h}]$ is derived by
solving a series of subproblems on local fine mesh with homogeneous boundary
conditions, its high frequency error may be suppressed apparently. To balance
the lower and higher frequency errors, a smooth step associated with the
intermediate approximation $[{\boldsymbol{u}}_{H,h},p_{H,h}]$ is resorted
through coarse grid correction as follows: find
$[{\boldsymbol{E}}_{\boldsymbol{u}}^{H},E_{p}^{H}]\in{\boldsymbol{X}}_{0}^{H}\times
Q^{H}$ such that
$\forall[{\boldsymbol{v}}_{H},q_{H}]\in{\boldsymbol{X}}_{0}^{H}\times Q^{H}$

$B([{\boldsymbol{E}}_{u}^{H},E_{p}^{H}],[{\boldsymbol{v}}_{H},q_{H}])=({\boldsymbol{f}},{\boldsymbol{v}}_{H})-B([{\boldsymbol{u}}_{H,h},p_{H,h}],[{\boldsymbol{v}}_{H},q_{H}]).$
(11)

To this end, an expected more accurate approximation than the coarse grid
approximation is obtained,

$[{\boldsymbol{u}}_{H}^{h},p_{H}^{h}]=[{\boldsymbol{u}}_{H,h}+{\boldsymbol{E}}_{\boldsymbol{u}}^{H},p_{H,h}+E_{p}^{H}]=[{\boldsymbol{u}}_{H}+\hat{\boldsymbol{u}}_{H,h}+{\boldsymbol{E}}_{\boldsymbol{u}}^{H},p_{H}+\hat{p}_{H,h}+E_{p}^{H}].$
(12)

In summary, we propose our local and parallel two-grid scheme for solving (1)
in the following.

Local and parallel two-grid scheme:

Step 0.

Deriving $[{\boldsymbol{u}}_{H},p_{H}]$ by solving (3);

Step 1.

Solving the equation (9) to get
$\\{[\hat{\boldsymbol{u}}_{H,h}^{j},\hat{p}_{H,h}^{j}]\\}_{j=1}^{N}$ for each
$j$, then constructing $[{\boldsymbol{u}}_{H,h},p_{H,h}]$ by formula (10);

Step 2.

Deriving $[{\boldsymbol{E}}_{\boldsymbol{u}}^{h},E_{p}^{h}]$ by solving (11)
and finally constructing $[{\boldsymbol{u}}_{H}^{h},p_{H}^{h}]$ by (12).

##  4 Theoretical Analysis and Suggested Local and Parallel Two-Grid Iterative
Scheme

As the basic step for analyzing the scheme above, we follow the idea discussed
in [19], and extend the local sub-problem (9) to the original domain $\Omega$.
For this purpose, we denote $\Gamma=\partial\Omega$,
$\Gamma_{j}=\partial\Omega_{j}\backslash\Gamma$, and introduce a trace space
${\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})={\boldsymbol{H}}_{0}^{1}(\Omega)|_{\Gamma_{j}}$
on $\Gamma_{j}$, which can be defined by interpolation (for example, see [26])

${\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})=\left\\{\begin{array}[]{ll}[{\boldsymbol{L}}^{2}(\Gamma_{j}),{\boldsymbol{H}}^{1}(\Gamma_{j})]_{\frac{1}{2}},&\mbox{when}\,\Gamma_{j}\,\mbox{is
a closed curve or surface},\\\\[0.0pt]
[{\boldsymbol{L}}^{2}(\Gamma_{j}),{\boldsymbol{H}}_{0}^{1}(\Gamma_{j})]_{\frac{1}{2}},&\mbox{when}\,\Gamma_{j}\,\mbox{is
a non-closed curve or surface}.\end{array}\right.$

We also denote by
${\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})={\boldsymbol{X}}_{0}^{h}(\Omega)|_{\Gamma_{j}}\subset{\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})$,
a finite dimensional trace space, equipped with the same norm as in
${\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})$, and its dual space by
${\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})=({\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j}))^{\prime}$,
which is equipped with the norm

$\|{\boldsymbol{\mu}}\|_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}=\sup\limits_{{\boldsymbol{v}}_{h}\in{\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})}\frac{\int_{\Gamma_{j}}{\boldsymbol{v}}_{h}{\boldsymbol{\mu}}}{\|{\boldsymbol{v}}_{h}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})}}.$

Then based on the fictitious domain method (see e.g., [25]), we can construct
a problem with multiplier as finding
$([\hat{\boldsymbol{u}}_{H,h}^{j},\hat{p}_{H,h}^{j}],{\boldsymbol{\xi}}^{j})\in{\boldsymbol{X}}_{0}^{h}\times
Q_{j}^{h}\times{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})$, such that

$\displaystyle
B([\hat{\boldsymbol{u}}_{H,h}^{j},\hat{p}_{H,h}^{j}],[{\boldsymbol{v}}_{h},q_{h}])+<{\boldsymbol{\xi}}^{j},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}+<{\boldsymbol{\mu}},\hat{\boldsymbol{u}}_{H,h}^{j}>_{\Gamma_{j}}$
(13)
$\displaystyle\qquad=({\boldsymbol{f}},\phi_{j}{\boldsymbol{v}}_{h})-B([{\boldsymbol{u}}_{H},p_{H}],[\phi_{j}{\boldsymbol{v}}_{h},\phi_{j}q_{h}]),$

holds for every
$([{\boldsymbol{v}}_{h},q_{h}],{\boldsymbol{\mu}})\in{\boldsymbol{X}}_{0}^{h}\times
Q_{j}^{h}\times{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})$, where

$<{\boldsymbol{\mu}},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}=\int_{\Gamma_{j}}{\boldsymbol{\mu}}{\boldsymbol{v}}_{h}ds,\quad\forall{\boldsymbol{\mu}}\in{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j}),\;{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}.$

In the following, we will show the equivalence between problems (9) and (13)
in two steps.

Firstly, by introducing a subspace of ${\boldsymbol{X}}_{0}^{h}$

${\boldsymbol{X}}_{j0}^{h}=\\{{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}:\;<{\boldsymbol{\mu}},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}=0\;\;\forall{\boldsymbol{\mu}}\in{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})\\}=\mbox{ker}(\gamma).$

the above system (13) reduces to finding
$[\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}]\in{\boldsymbol{X}}_{j0}^{h}\times
Q_{j}^{h}$ such that

$B([\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},\phi_{j}{\boldsymbol{v}}_{h})-B([{\boldsymbol{u}}_{H},p_{H}],[\phi_{j}{\boldsymbol{v}}_{h},\phi_{j}q_{h}]),\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{j0}^{h}\times
Q_{j}^{h},$

which is obviously well-posed, since it actually consists of two independent
Stokes problems defined in $\Omega_{j}$ and $\Omega\backslash\Omega_{j}$
respectively, with homogeneous boundary conditions for the velocity.

Secondly, for any given
${\boldsymbol{g}}\in{\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})$, we
introduce two auxiliary elliptic problems

$a_{\Omega_{j}}({\boldsymbol{u}}_{h}^{1},{\boldsymbol{v}}_{h})=0,\quad{\boldsymbol{u}}_{h}^{1}|_{\Gamma_{j}}={\boldsymbol{g}},\quad{\boldsymbol{u}}_{h}^{1}|_{\partial\Omega_{j}/\Gamma_{j}}=0,\quad\forall{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}(\Omega_{j}),$

and

$a_{\Omega\backslash\Omega_{j}}({\boldsymbol{u}}_{h}^{2},{\boldsymbol{v}}_{h})=0,\quad{\boldsymbol{u}}_{h}^{2}|_{\Gamma_{j}}={\boldsymbol{g}},\quad{\boldsymbol{u}}_{h}^{2}|_{\partial(\Omega/\Omega_{j})/\Gamma_{j}}=0,\quad\forall{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}(\Omega\backslash\Omega_{j}).$

Clearly, these two problems establish two mappings $\gamma_{1}^{-1}$ and
$\gamma_{2}^{-1}$ from ${\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})$ into
${\boldsymbol{X}}_{E}^{h}(\Omega_{j})$ and
${\boldsymbol{X}}_{E}^{h}(\Omega\backslash\Omega_{j})$, respectively, with

$\displaystyle{\boldsymbol{X}}_{E}^{h}(\Omega_{j})=\\{{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}^{h}(\Omega_{j}):{\boldsymbol{v}}_{h}|_{\partial\Omega_{j}\cap\partial\Omega}=0\\},$
$\displaystyle{\boldsymbol{X}}_{E}^{h}(\Omega\backslash\Omega_{j})=\\{{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}^{h}(\Omega\backslash\Omega_{j}):{\boldsymbol{v}}_{h}|_{\partial\Omega\backslash\partial\Omega_{j}}=0\\}.$

We can simply write as

${\boldsymbol{u}}_{h}^{1}=\gamma_{1}^{-1}{\boldsymbol{g}},\quad{\boldsymbol{u}}_{h}^{2}=\gamma_{2}^{-1}{\boldsymbol{g}},$

which have the following estimates

$\|\gamma_{1}^{-1}{\boldsymbol{g}}\|_{{\boldsymbol{H}}^{1}(\Omega_{j})},\;\|\gamma_{2}^{-1}{\boldsymbol{g}}\|_{H^{1}(\Omega\backslash\Omega_{j})}\lesssim\|{\boldsymbol{g}}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})}.$

We can also define a lifting operator $\gamma^{-1}$ from
${\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})$ into
${\boldsymbol{X}}_{0}^{h}$: for any given
${\boldsymbol{g}}\in{\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})$

$\gamma^{-1}{\boldsymbol{g}}=\left\\{\begin{array}[]{ll}\gamma_{1}^{-1}{\boldsymbol{g}},&\mbox{in}\;\Omega_{j},\\\
\gamma_{2}^{-1}{\boldsymbol{g}},&\mbox{in}\;\Omega\backslash\Omega_{j},\end{array}\right.$

which is the right inverse operator of the trace operator $\gamma$ from
${\boldsymbol{X}}_{0}^{h}$ onto
${\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})$, with the property of

$\|\gamma^{-1}{\boldsymbol{g}}\|_{{\boldsymbol{H}}^{1}(\Omega)}\lesssim\|{\boldsymbol{g}}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})}\quad\forall{\boldsymbol{g}}\in{\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j}).$

Then we have for any
${\boldsymbol{\mu}}\in{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})$

$\displaystyle\|{\boldsymbol{\mu}}\|_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}=\sup\limits_{{\boldsymbol{g}}\in{\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})}\frac{<{\boldsymbol{\mu}},{\boldsymbol{g}}>_{\Gamma_{j}}}{\|{\boldsymbol{g}}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})}}$
$\displaystyle\qquad\lesssim\sup\limits_{{\boldsymbol{g}}\in{\boldsymbol{H}}_{h}^{\frac{1}{2}}(\Gamma_{j})}\frac{<{\boldsymbol{\mu}},\gamma^{-1}{\boldsymbol{g}}>_{\Gamma_{j}}}{\|\gamma^{-1}{\boldsymbol{g}}\|_{{\boldsymbol{H}}^{1}(\Omega)}}\leq\sup\limits_{{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}(\Omega)}\frac{<{\boldsymbol{\mu}},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}}{\|{\boldsymbol{v}}_{h}\|_{{\boldsymbol{H}}^{1}(\Omega)}},$
(14)

which verifies that the bilinear form $<\cdot,\cdot>_{\Gamma_{j}}$ on
${\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})\times{\boldsymbol{X}}_{0}^{h}$
satisfies the inf-sup condition. Therefore, for any solution
$[\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}]\in{\boldsymbol{X}}_{0}^{h}\times
Q_{j}^{h}$ to (9), there exists a unique
${\boldsymbol{\xi}}^{j}\in{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})$
such that
$([\hat{\boldsymbol{u}}^{j}_{H,h},\hat{p}^{j}_{H,h}],{\boldsymbol{\xi}}^{j})$
satisfies the problem (13). In such sense, the problems (9) and (13) are
equivalent.

Now let us turn back to global and local problems (5) and (6) related with the
residuals of approximate solution. For the fine triangulation $T^{h}(\Omega)$
and the associated finite element space ${\boldsymbol{X}}_{0}^{h}$, their
corresponding Galerkin approximations are finding
$[\hat{\boldsymbol{u}}_{H},\hat{p}_{H}]\in{\boldsymbol{X}}_{0}^{h}\times
Q^{h}$ and
$[\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}]\in{\boldsymbol{X}}_{0}^{h}\times
Q^{h}$, $j=1,2,\cdots,N$, such that
$\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times Q^{h}$

$B([\hat{\boldsymbol{u}}_{H},\hat{p}_{H}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},{\boldsymbol{v}}_{h})-B([{\boldsymbol{u}}_{H},p_{H}],[{\boldsymbol{v}}_{h},q_{h}]),$
(15)

and

$B([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},\phi_{j}{\boldsymbol{v}}_{h})-B([{\boldsymbol{u}}_{H},p_{H}],[\phi_{j}{\boldsymbol{v}}_{h},\phi_{j}q_{h}]).$
(16)

We know that
$[\hat{\boldsymbol{u}}_{H},\hat{p}_{H}]=[\sum\limits_{j=1}^{N}\hat{\boldsymbol{u}}_{H}^{j},\sum\limits_{j=1}^{N}\hat{p}_{H}^{j}]$
is the Galerkin approximation of $[\hat{\boldsymbol{u}},\hat{p}]$ in
${\boldsymbol{X}}_{0}^{h}\times Q^{h}$, and meanwhile
$[{\boldsymbol{u}}_{h},p_{h}]=[{\boldsymbol{u}}_{H},p_{H}]+[\hat{\boldsymbol{u}}_{H},\hat{p}_{H}]\in{\boldsymbol{X}}_{0}^{h}\times
Q^{h}$ is the finite element approximation of $[{\boldsymbol{u}},p]$ on fine
mesh.

Then if denoting by

$\displaystyle{\boldsymbol{e}}^{j}_{{\boldsymbol{u}},h}=\hat{\boldsymbol{u}}^{j}_{H}-\hat{\boldsymbol{u}}^{j}_{H,h},\quad{\boldsymbol{e}}_{{\boldsymbol{u}},h}=\sum\limits_{j=1}^{N}{\boldsymbol{e}}^{j}_{{\boldsymbol{u}},h}={\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H,h},$
$\displaystyle e^{j}_{p,h}=\hat{p}^{j}_{H}-\hat{p}^{j}_{H,h},\quad
e_{p,h}=\sum\limits_{j=1}^{N}e^{j}_{p,h}=p_{h}-p_{H,h},$

the local and the global error of $[{\boldsymbol{u}}_{H,h},p_{H,h}]$,
respectively, subtracting (6) from (13) gives

$a({\boldsymbol{e}}^{j}_{{\boldsymbol{u}},h},{\boldsymbol{v}}_{h})+d(e_{p,h}^{j},{\boldsymbol{v}}_{h})+<{\boldsymbol{\xi}}^{j},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}=0,\quad\forall{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h},$
(17)

and further summing over all $j$’s yields

$a({\boldsymbol{e}}_{{\boldsymbol{u}},h},{\boldsymbol{v}}_{h})+d(e_{p,h},{\boldsymbol{v}}_{h})+\sum\limits_{j=1}^{N}<{\boldsymbol{\xi}}^{j},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}=0,\quad\forall{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}.$
(18)

In the following, we will carry out estimates on the local quantities
${\boldsymbol{e}}^{j}_{{\boldsymbol{u}},h},e_{p,h}^{j}$ and
${\boldsymbol{\xi}}^{j}$, which are very important in the error analysis of
our scheme. Firstly, for ${\boldsymbol{\xi}}^{j}$, by the previously defined
operator $\gamma^{-1}$ and the fact that for every point $x\in\Omega$, there
exists a positive integer $\kappa$, which is independent of $N$ and $x$, such
that each $x$ belongs to at most $\kappa$ different $\Omega_{j}$, we can
easily derive the estimates as follows.

######  Lemma 1.

The multiplier ${\boldsymbol{\xi}}^{j}$ in (13), also in (17), satisfies

$\|{\boldsymbol{\xi}}^{j}\|_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}\lesssim\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega}+\|e_{p,h}^{j}\|_{\Omega},$

and

$\sum\limits_{j=1}^{N}<{\boldsymbol{\xi}}^{j},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}\lesssim\kappa^{\frac{1}{2}}\left(\sum\limits_{j=1}^{N}\|{\boldsymbol{\xi}}^{j}\|^{2}_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}\right)^{\frac{1}{2}}\|\nabla{\boldsymbol{v}}_{h}\|_{\Omega}\quad\forall{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}.$

###### Proof.

The first estimate is a direct result of property (14). For the second
estimate, by the definition of
$\|\cdot\|_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}$, we have
$\forall{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}$

$\displaystyle\sum\limits_{j=1}^{N}<{\boldsymbol{\xi}}^{j},{\boldsymbol{v}}_{h}>_{\Gamma_{j}}\leq\sum\limits_{j=1}^{N}\|{\boldsymbol{\xi}}^{j}\|_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}\|{\boldsymbol{v}}_{h}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\Gamma_{j})}\lesssim\sum\limits_{j=1}^{N}\|{\boldsymbol{\xi}}^{j}\|_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}\|{\boldsymbol{v}}_{h}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\partial\Omega_{j})}$
$\displaystyle\qquad\lesssim\sum\limits_{j=1}^{N}\|{\boldsymbol{\xi}}^{j}\|_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})}\|{\boldsymbol{v}}_{h}\|_{{\boldsymbol{H}}^{1}(\Omega_{j})}\leq(\sum\limits_{j=1}^{N}\|{\boldsymbol{\xi}}^{j}\|^{2}_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})})^{\frac{1}{2}}(\sum\limits_{j=1}^{N}\|{\boldsymbol{v}}_{h}\|^{2}_{{\boldsymbol{H}}^{1}(\Omega_{j})})^{\frac{1}{2}}$
$\displaystyle\qquad\leq\kappa^{\frac{1}{2}}(\sum\limits_{j=1}^{N}\|{\boldsymbol{\xi}}^{j}\|^{2}_{{\boldsymbol{H}}_{h}^{-\frac{1}{2}}(\Gamma_{j})})^{\frac{1}{2}}\|\nabla{\boldsymbol{v}}_{h}\|_{\Omega}.$

This is the second estimate of the lemma. ∎

As to the quantity of
$\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega}+\|e_{p,h}^{j}\|_{\Omega}$,
combining the identities

$\displaystyle\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega}^{2}=\|\nabla(\hat{\boldsymbol{u}}_{H}^{j}-\hat{\boldsymbol{u}}_{H,h}^{j})\|_{\Omega}^{2}=\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega\backslash\Omega_{j}}^{2}+\|\nabla(\hat{\boldsymbol{u}}_{H}^{j}-\hat{\boldsymbol{u}}_{H,h}^{j})\|_{\Omega_{j}}^{2},$
$\displaystyle\|e_{p,h}^{j}\|_{\Omega}^{2}=\|\hat{p}_{H}^{j}-\hat{p}_{H,h}^{j}\|_{\Omega}^{2}=\|\hat{p}_{H}^{j}\|_{\Omega\backslash\Omega_{j}}^{2}+\|\hat{p}_{H}^{j}-\hat{p}_{H,h}^{j}\|_{\Omega_{j}}^{2},$

and the equations satisfied by
$[{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}|_{\Omega_{j}},e_{p,h}^{j}|_{\Omega_{j}}]$

$\left\\{\begin{array}[]{ll}a_{\Omega_{j}}({\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}|_{\Omega_{j}},{\boldsymbol{v}}_{h})+d_{\Omega_{j}}(e_{p,h}^{j}|_{\Omega_{j}},{\boldsymbol{v}}_{h})=0,&\forall{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}(\Omega_{j}),\\\
d_{\Omega_{j}}(q_{h},{\boldsymbol{e}}^{j}_{{\boldsymbol{u}},h}|_{\Omega_{j}})=0,&\forall
q_{h}\in Q^{h}(\Omega_{j}),\\\
{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}|_{\partial\Omega_{j}}=\hat{\boldsymbol{u}}_{H}^{j}|_{\partial\Omega_{j}},\end{array}\right.$

we have by virtue of $\hat{\boldsymbol{u}}_{H}^{j}\in{\boldsymbol{X}}_{0}^{h}$

$\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega_{j}}+\|e_{p,h}^{j}\|_{\Omega_{j}}\lesssim\|\hat{\boldsymbol{u}}_{H}^{j}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\partial\Omega_{j})}\leq\|\hat{\boldsymbol{u}}_{H}^{j}\|_{{\boldsymbol{H}}^{\frac{1}{2}}(\partial(\Omega\backslash\Omega_{j}))}{\leq}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega\backslash\Omega_{j}},$

which is arriving at

$\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega_{j}}+\|e_{p,h}^{j}\|_{\Omega_{j}}\lesssim\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega\backslash\Omega_{j}}.$
(19)

Then we present the following two results related with the Sobolev space
${\boldsymbol{H}}_{0}^{1}(\Omega)$ to make further analysis concerning on
(19).

######  Lemma 2.

Let ${\cal D}\subset\Omega$ with $\mbox{diam}({\cal D})\approxeq H$ be any
convex subdomain of $\Omega$. Then we have

$\|{\boldsymbol{w}}\|^{2}_{{\boldsymbol{L}}^{2}(\partial{\cal D})}\lesssim
H^{d-1}\beta_{d}(H)\|\nabla{\boldsymbol{w}}\|^{2}_{\Omega\backslash{\cal
D}}\quad\forall{\boldsymbol{w}}\in{\boldsymbol{H}}_{0}^{1}(\Omega),$

where $\beta_{d}(H)=H^{-1}$ when $d=3$ and $\beta_{d}(H)=|\ln H|$ when $d=2$.

###### Proof.

Since ${\cal D}$ and $\Omega$ are convex domains, there exists a point
$P\in{\cal D}$ such that $\mbox{dist}(P,\partial{\cal D})\approxeq H$ and
$\overline{PP^{\prime}}\subset{\cal D}$ for any $P^{\prime}\in\partial{\cal
D}$ and $\overline{PP^{\prime}}\subset\Omega$ for any
$P^{\prime}\in\partial\Omega$. If we denote $\Omega^{c}=\Omega\backslash{\cal
D}$, we can establish a local polar or spherical coordinate $(\rho,\omega)$
with origin at this point $P$ corresponding to 2-D to 3-D case, see Fig.1 for
the sketch of such setting in the 2-D case.

Fig. 1: Local polar coordinate in 2-dimensional domain.

Here $\rho$ is the radius distant between any point $P^{\prime}$ in $\Omega$
and the origin $P$, and $\omega=\omega_{1}$ in 2-D case and
$\omega=(\omega_{1},\omega_{2})$ in 3-D case. The Jacobi determinant of the
transformation between the Cartesian coordinate $(x_{1},\cdots,x_{d})$ and
$(\rho,\omega)$ is

$J=\frac{D(x_{1},\cdots,x_{d})}{D(\rho,\omega_{1}\cdots,\omega_{d-1})}=\rho^{d-1}\delta(\omega),\quad\mbox{with}\quad\delta(\omega)=\sin^{d-2}\omega_{1}.$

Moreover, we use $\rho_{\Omega}(\omega)$ and $\rho_{\cal D}(\omega)$ to
characterize the boundary points of $\Omega$ and ${\cal D}$, respectively.

For any ${\boldsymbol{w}}\in{\boldsymbol{H}}_{0}^{1}(\Omega)$, since
$\rho_{\cal D}(\omega)\approxeq H$ and $H\lesssim\rho_{\Omega}(\omega)\lesssim
1$, and ${\boldsymbol{w}}(\rho_{\Omega}(\omega),\omega)=0$, we have

$\displaystyle|{\boldsymbol{w}}(\rho_{\cal
D}(\omega),\omega)|=|-\int_{\rho_{\Omega}(\omega)}^{\rho_{\cal
D}(\omega)}\frac{\partial{\boldsymbol{w}}}{\partial\rho}d\rho|$
$\displaystyle\qquad\leq(\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D}(\omega)}\frac{1}{\rho^{d-1}}d\rho)^{\frac{1}{2}}(\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D}(\omega)}\rho^{d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho}|^{2}d\rho)^{\frac{1}{2}}$
$\displaystyle\qquad\lesssim\beta_{d}^{\frac{1}{2}}(H)(\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D}(\omega)}\rho^{d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho}|^{2}d\rho)^{\frac{1}{2}}.$

Thus, if we denote by $\omega_{d}$ the unit ball in $R^{d}$, we have

$\displaystyle\|{\boldsymbol{w}}\|^{2}_{{\boldsymbol{L}}^{2}(\partial{\cal
D})}=\int_{\partial\omega_{d}}(\rho_{\cal
D}(\omega))^{d-1}|{\boldsymbol{w}}(\rho,\omega)|^{2}\delta(\omega)d\omega$
$\displaystyle\qquad\lesssim\int_{\partial\omega_{d}}(\rho_{\cal
D}(\omega))^{d-1}\beta_{d}(H)(\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D}(\omega)}\rho^{d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho}|^{2}d\rho)\delta(\omega)d\omega$
$\displaystyle\qquad\lesssim
H^{d-1}\beta_{d}(H)\|\nabla{\boldsymbol{w}}\|^{2}_{\Omega\backslash{\cal D}}.$

The proof is complete. ∎

######  Lemma 3.

If ${\cal D}^{\prime}\subset{\cal D}\subset\Omega$ are any two convex
subdomains of $\Omega$ with $\mbox{diam}({\cal D})\approxeq H$,
$\mbox{dist}(\partial{\cal D},\partial{\cal D}^{\prime})\approxeq h$ and
$h<H$. Then we have

$\|{\boldsymbol{w}}\|^{2}_{{\cal D}\backslash{\cal D}^{\prime}}\lesssim
hH^{d-1}\beta_{d}(H)\|\nabla{\boldsymbol{w}}\|^{2}_{\Omega\backslash{\cal
D}^{\prime}}\quad\forall{\boldsymbol{w}}\in{\boldsymbol{H}}_{0}^{1}(\Omega).$

###### Proof.

Like what has been done in the proof of Lemma 2, for ${\cal
D}^{\prime}\subset{\cal D}\subset\Omega$, we can choose some $P\in{\cal
D}^{\prime}$ such that $\mbox{dist}(P,\partial{\cal
D}^{\prime}),\mbox{dist}(P,\partial{\cal D})\approx H$, then we can establish
a local polar or spherical coordinate $(\rho,\omega)$ with origin at this
point $P$ corresponding to 2-D to 3-D case. And from the proof of Lemma 2, we
know for any point $(\rho,\omega)\in{\cal D}\backslash{\cal D}^{\prime}$:

$\displaystyle|{\boldsymbol{w}}(\rho,\omega)|=|-\int_{\rho_{\Omega}(\omega)}^{\rho}\frac{\partial{\boldsymbol{w}}}{\partial\rho^{\prime}}d\rho^{\prime}|\leq(\int^{\rho_{\Omega}(\omega)}_{\rho}\frac{1}{\rho^{\prime
d-1}}d\rho^{\prime})^{\frac{1}{2}}(\int^{\rho_{\Omega}(\omega)}_{\rho}\rho^{\prime
d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho^{\prime}}|^{2}d\rho^{\prime})^{\frac{1}{2}}$
$\displaystyle\qquad\lesssim\beta_{d}^{\frac{1}{2}}(H)(\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D^{\prime}}(\omega)}\rho^{\prime
d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho^{\prime}}|^{2}d\rho^{\prime})^{\frac{1}{2}}.$

Then we have

$\displaystyle\|{\boldsymbol{w}}\|_{D\backslash
D^{\prime}}^{2}=\int_{\partial\omega_{d}}\int_{\rho_{D^{\prime}(\omega)}}^{\rho_{D}(\omega)}{\boldsymbol{w}}^{2}\rho^{d-1}\delta(\omega)d\rho
d\omega$
$\displaystyle\quad\lesssim\int_{\partial\omega_{d}}\int_{\rho_{D^{\prime}(\omega)}}^{\rho_{D}(\omega)}\beta_{d}(H)[\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D^{\prime}}(\omega)}\rho^{\prime
d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho^{\prime}}|^{2}d\rho^{\prime}]\rho^{d-1}\delta(\omega)d\rho
d\omega$ $\displaystyle\quad\lesssim
H^{d-1}\beta_{d}(H)\int_{\partial\omega_{d}}\int_{\rho_{D^{\prime}(\omega)}}^{\rho_{D}(\omega)}[\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D^{\prime}}(\omega)}\rho^{\prime
d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho^{\prime}}|^{2}d\rho^{\prime}]\delta(\omega)d\rho
d\omega$ $\displaystyle\quad\leq
H^{d-1}\beta_{d}(H)\int_{\partial\omega_{d}}(\int_{\rho_{D^{\prime}(\omega)}}^{\rho_{D}(\omega)}[\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D^{\prime}}(\omega)}\rho^{\prime
d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho^{\prime}}|^{2}d\rho^{\prime}]^{2}\delta^{2}(\omega)d\rho)^{\frac{1}{2}}(\int_{\rho_{D^{\prime}(\omega)}}^{\rho_{D}(\omega)}d\rho)^{\frac{1}{2}}d\omega$
$\displaystyle\quad\lesssim
hH^{d-1}\beta_{d}(H)\int_{\partial\omega_{d}}\int^{\rho_{\Omega}(\omega)}_{\rho_{\cal
D^{\prime}}(\omega)}\rho^{\prime
d-1}|\frac{\partial{\boldsymbol{w}}}{\partial\rho^{\prime}}|^{2}d\rho^{\prime}\delta(\omega)d\omega=hH^{d-1}\beta_{d}(H)\|\nabla{\boldsymbol{w}}\|^{2}_{\Omega\backslash
D^{\prime}},$

which completes the proof. ∎

Right now, we can give a more rigorous estimate for
$\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega}^{2}+\|e_{p,h}^{j}\|_{\Omega}^{2}$
based on (19), which will play a crucial role in the analysis of this section.

######  Lemma 4.

Defining

$\alpha_{d}=\left\\{\begin{array}[]{ll}\frac{c}{|\ln H|^{2}},&d=2,\\\
\frac{c}{|\ln H|},&d=3,\end{array}\right.$

where $c>0$ is a positive constant independent of $H$, $h$ and any subdomain
$\Omega_{j}$. Then we have

$\displaystyle\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega}+\|e_{p,h}^{j}\|_{\Omega}\lesssim
H^{\alpha_{d}}(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}+\|\hat{p}_{H}^{j}\|_{\Omega}).$
(20)

###### Proof.

The proof of this lemma is divided into four steps.

_Step 1._ Firstly, we use the same argument to split the region
$\Omega_{j}\backslash D_{j}$ for each $j$ as in [19]. We recall the sketch of
a domain partition in Fig. 2.

Fig. 2: Division of the region $\Omega_{j}\backslash D_{j}$

Right now, we can define a series of disjoint annular zones as
$\Omega_{j}^{0},\Omega_{j}^{1},\cdots,\Omega_{j}^{\cal M}$, and construct a
sequence of subdomains as

$\tilde{\Omega}_{j}^{k}=\bigcup\limits_{i=0}^{k}\Omega_{j}^{i},\quad
k=0,1,2,\cdots,{\cal M},$

with the property

$\tilde{\Omega}_{j}^{0}\subset\subset\tilde{\Omega}_{j}^{1}\subset\subset\cdots\subset\subset\tilde{\Omega}_{j}^{\cal
M}.$

Then we define

$\gamma_{j}^{k}=\partial\tilde{\Omega}_{j}^{k}\backslash\partial\Omega,\quad\Gamma_{j}^{k}=\partial\tilde{\Omega}_{j}^{k}\backslash\gamma_{j}^{k},\quad\mbox{and}\quad\partial\tilde{\Omega}_{j}^{k}=\gamma_{j}^{k}\cup\Gamma_{j}^{k},\quad\quad
k=0,1,2,\cdots,{\cal M}.$

_Step 2._ Note that, since $\texttt{supp}\phi_{j}=D_{j}$, the right hand sides
in (16) will be equal to 0 for the testing functions in
${\boldsymbol{X}}_{0}^{h}(\tilde{\Omega}_{j}^{k})\times
M^{h}(\tilde{\Omega}_{j}^{k})$ defined on the domain
$\tilde{\Omega}_{j}^{k}\subset\subset\Omega\backslash D_{j}$. By denoting
$M_{0}^{h}(\tilde{\Omega}_{j}^{k})=M^{h}(\tilde{\Omega}_{j}^{k})\cap
H^{1}_{0}(\tilde{\Omega}_{j}^{k})$, we clearly have

$\displaystyle
B_{\tilde{\Omega}_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}],[{\boldsymbol{v}}_{h},q_{h}])=0\quad\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}(\tilde{\Omega}_{j}^{k})\times
M_{0}^{h}(\tilde{\Omega}_{j}^{k}).$ (21)

For any $\theta\in R$, since
${\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}(\tilde{\Omega}_{j}^{k})$, we
have

$d_{\tilde{\Omega}_{j}^{k}}(\theta,{\boldsymbol{v}}_{h})=\int_{\tilde{\Omega}_{j}^{k}}\theta\nabla\cdot{\boldsymbol{v}}_{h}dx=\theta\int_{\tilde{\Omega}_{j}^{k}}\nabla\cdot{\boldsymbol{v}}_{h}dx=\theta\int_{\partial\tilde{\Omega}_{j}^{k}}{\boldsymbol{v}}_{h}\cdot
nds=0.$

By the definition of $B$, we know

$0=B_{\tilde{\Omega}_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}],[{\boldsymbol{v}}_{h},q_{h}])=B_{\tilde{\Omega}_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}],[{\boldsymbol{v}}_{h},q_{h}])-d_{\tilde{\Omega}_{j}^{k}}(\theta,{\boldsymbol{v}}_{h})=B_{\tilde{\Omega}_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta],[{\boldsymbol{v}}_{h},q_{h}]).$

Then, we can split it into two parts since
${\tilde{\Omega}_{j}^{k}}=\tilde{\Omega}_{j}^{k-1}\cup\Omega_{j}^{k}$ as

$0=B_{\tilde{\Omega}_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta],[{\boldsymbol{v}}_{h},q_{h}])=B_{\tilde{\Omega}_{j}^{k}-1}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta],[{\boldsymbol{v}}_{h},q_{h}])+B_{\Omega_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta],[{\boldsymbol{v}}_{h},q_{h}]),$

therefore, by simply using

$\displaystyle\tilde{\boldsymbol{v}}_{h}^{k}={\boldsymbol{v}}_{h}|_{\tilde{\Omega}_{j}^{k}},\quad\tilde{\boldsymbol{v}}_{h}^{k-1}={\boldsymbol{v}}_{h}|_{\tilde{\Omega}_{j}^{k-1}},\quad{\boldsymbol{v}}_{h}^{k}={\boldsymbol{v}}_{h}|_{\Omega_{j}^{k}},$
$\displaystyle\tilde{q}_{h}^{k}=q_{h}|_{\tilde{\Omega}_{j}^{k}},\quad\tilde{q}_{h}^{k-1}=q_{h}|_{\tilde{\Omega}_{j}^{k-1}},\quad
q_{h}^{k}=q_{h}|_{\Omega_{j}^{k}},$

we have for any $\theta\in R$

$B_{\tilde{\Omega}_{j}^{k-1}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta],[\tilde{\boldsymbol{v}}_{h}^{k-1},\tilde{q}_{h}^{k-1}])=-B_{\Omega_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta],[{\boldsymbol{v}}_{h}^{k},q_{h}^{k}]).$
(22)

_Step 3._ We firstly recall the following classical result (see for instance
(2.11) in [24]):

$\displaystyle||q||_{0,\Omega}=\inf\limits_{c\in R}||q+c||_{0,\Omega}\ \ \
\forall q\in L^{2}_{0}(\Omega).$

Then, for convenience of expression, we denote by
${\boldsymbol{X}}_{E}^{h}(\tilde{\Omega}_{j}^{k})={\boldsymbol{X}}_{0}^{h}(\Omega\backslash
D_{j})|_{\tilde{\Omega}_{j}^{k}}$,
$M_{E}^{h}(\tilde{\Omega}_{j}^{k})=M_{0}^{h}(\Omega\backslash
D_{j})|_{\tilde{\Omega}_{j}^{k}}$, and introduce

$\tilde{\theta}^{k}_{j,q_{h}}=\frac{1}{|\tilde{\Omega}_{j}^{k}|}\int_{\tilde{\Omega}_{j}^{k}}q_{h},\quad\theta^{k}_{j,q_{h}}=\frac{1}{|\Omega_{j}^{k}|}\int_{\Omega_{j}^{k}}q_{h}.$

Since
${\boldsymbol{X}}_{0}^{h}(\tilde{\Omega}_{j}^{k})\subset{\boldsymbol{X}}_{E}^{h}(\tilde{\Omega}_{j}^{k})$,
$M_{E}^{h}(\tilde{\Omega}_{j}^{k})\subset M^{h}(\tilde{\Omega}_{j}^{k})$ and
${\boldsymbol{X}}_{0}^{h}(\tilde{\Omega}_{j}^{k})\times
M^{h}(\tilde{\Omega}_{j}^{k})/R$ is a stable finite element pair, we know
$\Sigma_{j}^{k}={\boldsymbol{X}}_{E}^{h}(\tilde{\Omega}_{j}^{k})\times
M_{E}^{h}(\tilde{\Omega}_{j}^{k})/R$ is also a stable finite element pair,
that is

$\sup\limits_{\tilde{\boldsymbol{v}}_{h}^{k}\in{\boldsymbol{X}}_{E}^{h}(\tilde{\Omega}_{j}^{k})}\frac{d_{\tilde{\Omega}_{j}^{k}}(\tilde{q}_{h}^{k},\tilde{\boldsymbol{v}}_{h}^{k})}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}\gtrsim\inf\limits_{\theta\in
R}\|\tilde{q}_{h}^{k}-\theta\|_{\tilde{\Omega}_{j}^{k}}=\|\tilde{q}_{h}^{k}-\tilde{\theta}_{j,q_{h}}^{k}\|_{\tilde{\Omega}_{j}^{k}}.$

We define a smooth function $\phi\in{C^{\infty}(\tilde{\Omega}_{j}^{k})}$ such
that $\phi|_{\gamma_{j}^{k}}=0$ and

$\mbox{supp}\,\phi=\tilde{\Omega}_{j}^{k},\quad\phi(x)\equiv 1\;\forall
x\in\tilde{\Omega}_{j}^{k-1},\quad 0\leq\phi\leq
1,\cancel{\quad\mbox{and}\quad|\phi(x)|\lesssim C}.$

In the following, we still use $I_{h}$ to denote the scalar valued Lagrange
finite element interpolation, which share the same property in A1. Firstly, we
can arrive at

$\displaystyle\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{k}}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{k}\|_{\tilde{\Omega}_{j}^{k}}=\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{k}}+\inf\limits_{\theta\in
R}\|\hat{p}_{H}^{j}-\theta\|_{\tilde{\Omega}_{j}^{k}}$
$\displaystyle\lesssim\sup\limits_{[\tilde{\boldsymbol{v}}_{h}^{k},\tilde{q}_{h}^{k}]\in\Sigma_{j}^{k}}\frac{B_{\tilde{\Omega}_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}],[\tilde{\boldsymbol{v}}_{h}^{k},\tilde{q}_{h}^{k}])}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}+\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\lesssim\sup\limits_{[\tilde{\boldsymbol{v}}_{h}^{k},\tilde{q}_{h}^{k}]\in\Sigma_{j}^{k}}\frac{B_{\tilde{\Omega}_{j}^{k-1}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}],[\tilde{\boldsymbol{v}}_{h}^{k-1},\tilde{q}_{h}^{k-1}])+B_{\Omega_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}],[{\boldsymbol{v}}_{h}^{k},q_{h}^{k}])}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}+\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\qquad=\sup\limits_{[\tilde{\boldsymbol{v}}_{h}^{k},\tilde{q}_{h}^{k}]\in\Sigma_{j}^{k}}\frac{-B_{\Omega_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}],[I_{h}(\phi
v_{h}^{k}),I_{h}(\phi
q_{h}^{k})])+B_{\Omega_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}],[{\boldsymbol{v}}_{h}^{k},q_{h}^{k}])}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}+\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\lesssim\sup\limits_{[\tilde{\boldsymbol{v}}_{h}^{k},\tilde{q}_{h}^{k}]\in\Sigma_{j}^{k}}\frac{-B_{\Omega_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}],[I_{h}(\phi
v_{h}^{k}),I_{h}(\phi
q_{h}^{k})])}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}+\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}+\sup\limits_{[\tilde{\boldsymbol{v}}_{h}^{k},\tilde{q}_{h}^{k}]\in\Sigma_{j}^{k}}\frac{B_{\Omega_{j}^{k}}([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}],[{\boldsymbol{v}}_{h}^{k},q_{h}^{k}])}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}+\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle:=I_{1}+I_{2}.$

In the following, we will carry out further estimates for the two items.

$\displaystyle I_{1}$
$\displaystyle\qquad\lesssim\sup\limits_{{\tilde{\boldsymbol{v}}_{h}^{k}\in{\boldsymbol{X}}_{E}^{h}(\tilde{\Omega}_{j}^{k})}}\frac{|a_{\Omega_{j}^{k}}(\hat{\boldsymbol{u}}_{H}^{j},I_{h}(\phi{\boldsymbol{v}}_{h}^{k}))+d_{\Omega_{j}^{k}}(\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k},I_{h}(\phi{\boldsymbol{v}}_{h}^{k}))|}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\qquad\quad+\sup\limits_{\tilde{q}_{h}^{k}\in{M_{E}^{h}(\tilde{\Omega}_{j}^{k})/R}}\frac{|d_{\Omega_{j}^{k}}(I_{h}(\phi
q_{h}^{k}),\hat{u}_{H}^{j})|}{\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\qquad\lesssim{\sup\limits_{{\tilde{\boldsymbol{v}}_{h}^{k}\in{\boldsymbol{X}}_{E}^{h}(\tilde{\Omega}_{j}^{k})}}}\frac{\|\nabla
I_{h}(\phi{\boldsymbol{v}}_{h}^{k})\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}+\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}\|\nabla\cdot(I_{h}(\phi{\boldsymbol{v}}_{h}^{k}))\|_{\Omega_{j}^{k}}}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\qquad\quad+{\sup\limits_{\tilde{q}_{h}^{k}\in
M_{E}^{h}(\tilde{\Omega}_{j}^{k})}}\frac{\|I_{h}\phi
q_{h}^{k}\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}}{\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\qquad\lesssim{\sup\limits_{{\tilde{\boldsymbol{v}}_{h}^{k}\in{\boldsymbol{X}}_{E}^{h}(\tilde{\Omega}_{j}^{k})}}}\frac{h^{-1}\|{\boldsymbol{v}}_{h}^{k}\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}+h^{-1}\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}\|{\boldsymbol{v}}_{h}^{k}\|_{\Omega_{j}^{k}}}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\qquad\quad+{\sup\limits_{\tilde{q}_{h}^{k}\in
M_{E}^{h}(\tilde{\Omega}_{j}^{k})}}\frac{\|q_{h}^{k}\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}}{\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}},$

and

$\displaystyle I_{2}$
$\displaystyle\qquad\lesssim\sup\limits_{\tilde{\boldsymbol{v}}_{h}^{k}\in{\boldsymbol{X}}_{0}(\Omega\backslash
D_{j})|_{\tilde{\Omega}_{j}^{k}}}\frac{|a_{\Omega_{j}^{k}}(\hat{\boldsymbol{u}}_{H}^{j},{\boldsymbol{v}}_{h}^{k})+d_{\Omega_{j}^{k}}(\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k},{\boldsymbol{v}}_{h}^{k})|}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}+\sup\limits_{\tilde{q}_{h}^{k}\in
M^{h}(\tilde{\Omega}_{j}^{k})}\frac{|d_{\Omega_{j}^{k}}(q_{h}^{k},\hat{u}_{H}^{j})|}{\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}$
$\displaystyle\qquad\lesssim\frac{h^{-1}\|{\boldsymbol{v}}_{h}^{k}\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}+h^{-1}\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}\|{\boldsymbol{v}}_{h}^{k}\|_{\Omega_{j}^{k}}}{\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}+\frac{\|q_{h}^{k}\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}}{\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}}.$

We can choose ${\cal D}=\Omega\backslash\tilde{\Omega}_{j}^{k-1},\ {\cal
D}\backslash{\cal D^{\prime}}=\Omega_{j}^{k},\ \Omega\backslash{\cal
D}^{\prime}=\tilde{\Omega}_{j}^{k}$ in Lemma 3. Noting that $\mbox{diam}({\cal
D})\approxeq H$, $\mbox{dist}(\partial{\cal D},\partial{\cal
D}^{\prime})\approxeq h$ by the definitions of $\tilde{\Omega}_{j}^{k-1}$ and
$\Omega_{j}^{k}$, then we can derive by virtue of (A2) and Lemma 3:

$\displaystyle
h^{-1}\|{\boldsymbol{v}}_{h}^{k}\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}\lesssim
h^{-\frac{1}{2}}H^{\frac{d-1}{2}}\beta_{d}^{\frac{1}{2}}(H)\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}},$
$\displaystyle
h^{-1}\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}\|{\boldsymbol{v}}_{h}^{k}\|_{\Omega_{j}^{k}}\lesssim
h^{-\frac{1}{2}}H^{\frac{d-1}{2}}\beta_{d}^{\frac{1}{2}}(H)\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}\|\nabla\tilde{\boldsymbol{v}}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}.$

From $q_{h}\in M_{0}^{h}(\tilde{\Omega}_{j}^{k})$, Lemma 3 and the inverse
inequality, we also have

$\displaystyle\|q_{h}^{k}\|_{\Omega_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}\lesssim
h^{\frac{1}{2}}H^{\frac{d-1}{2}}\beta_{d}^{\frac{1}{2}}(H)\|\nabla\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}\lesssim
h^{-\frac{1}{2}}H^{\frac{d-1}{2}}\beta_{d}^{\frac{1}{2}}(H)\|\tilde{q}_{h}^{k}\|_{\tilde{\Omega}_{j}^{k}}\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}.$

Combining the estimates above leads to

$\displaystyle\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{k}}+\inf\limits_{\theta\in
R}\|\hat{p}_{H}^{j}-\theta\|_{\tilde{\Omega}_{j}^{k}}=\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{k}}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{k}\|_{\tilde{\Omega}_{j}^{k}}$
$\displaystyle\qquad\lesssim
h^{-\frac{1}{2}}H^{\frac{d-1}{2}}\beta_{d}^{\frac{1}{2}}(H)(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}+\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}).$

Therefore we have

$\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{k}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{k}\|_{\tilde{\Omega}_{j}^{k}}^{2}\lesssim
h^{-1}H^{d-1}\beta_{d}(H)(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}^{2}+\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}^{2}).$

It is obvious that

$\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{k-1}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{k-1}\|_{\tilde{\Omega}_{j}^{k-1}}^{2}\lesssim
h^{-1}H^{d-1}\beta_{d}(H)(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}^{2}+\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}^{2}),$

or equivalently, there exists a mesh independent constant $c>0$ so that

$\displaystyle\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{k}}^{2}+\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{k}\|_{\Omega_{j}^{k}}^{2}\geq{chH^{1-d}\beta_{d}^{-1}(H)(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{k-1}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{k-1}\|_{\tilde{\Omega}_{j}^{k-1}}^{2})}.$
(23)

_Step 4._ By using the above inequality successively, we can deduce that

$\displaystyle\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{\cal
M}}^{2}+\|\hat{p}_{H}^{j}\|_{\tilde{\Omega}_{j}^{\cal M}}^{2}$
$\displaystyle\qquad\geq(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{{\cal
M}-1}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{{\cal
M}-1}\|_{\tilde{\Omega}_{j}^{{\cal
M}-1}}^{2})+(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega_{j}^{\cal
M}}^{2}+\|\hat{p}_{H}^{j}-\theta_{j,\hat{p}_{H}^{j}}^{\cal
M}\|_{\Omega_{j}^{\cal M}}^{2})$
$\displaystyle\qquad\geq(1+chH^{1-d}\beta_{d}^{-1}(H))(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{{\cal
M}-1}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{{\cal
M}-1}\|_{\tilde{\Omega}_{j}^{{\cal M}-1}}^{2})\geq\cdots$
$\displaystyle\qquad\geq(1+chH^{1-d}\beta_{d}^{-1}(H))^{\cal
M}(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{0}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{0}\|_{\tilde{\Omega}_{j}^{0}}^{2})$
$\displaystyle\qquad=(1+chH^{1-d}\beta_{d}^{-1}(H))^{\cal
M}(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega\backslash\Omega_{j}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{0}\|_{\Omega\backslash\Omega_{j}}^{2}).$

Noting that

$\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\tilde{\Omega}_{j}^{\cal
M}}^{2}+\|\hat{p}_{H}^{j}\|_{\tilde{\Omega}_{j}^{\cal
M}}^{2}=\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega\backslash
D_{j}}^{2}+\|\hat{p}_{H}^{j}\|^{2}_{\Omega\backslash
D_{j}}\leq\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}^{2}+\|\hat{p}_{H}^{j}\|^{2}_{\Omega},$

we can get

$\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega\backslash\Omega_{j}}^{2}+\|\hat{p}_{H}^{j}-\tilde{\theta}_{j,\hat{p}_{H}^{j}}^{0}\|_{\Omega\backslash\Omega_{j}}^{2}\leq(1+chH^{1-d}\beta_{d}^{-1}(H))^{-{\cal
M}}(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}^{2}+\|\hat{p}_{H}^{j}\|^{2}_{\Omega}).$

Finally, due to (19), we can arrive at

$\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega_{j}}^{2}+\|e_{p,h}^{j}\|_{\Omega_{j}}^{2}\lesssim(1+chH^{1-d}\beta_{d}^{-1}(H))^{-{\cal
M}}(\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}^{2}+\|\hat{p}_{H}^{j}\|^{2}_{\Omega}).$

Since ${\cal M}\approxeq\frac{H}{h}$, for small $H$, simple calculation leads
to the result of this lemma. ∎

######  Remark 4.1.

For small $H$, due to a basic inequality
$\lambda/2\leq\ln(1+\lambda)\leq\lambda$ for any $0\leq\lambda\leq 1$, we have

$\displaystyle(1+chH^{1-d}\beta_{d}^{-1}(H))^{-{\cal M}}$ $\displaystyle=$
$\displaystyle\exp\left(-{\cal M}\ln(1+chH^{1-d}\beta_{d}^{-1}(H))\right)$
$\displaystyle\approxeq$ $\displaystyle\exp\left(-{\cal
M}(chH^{1-d}\beta_{d}^{-1}(H))\right)$ $\displaystyle\approxeq$
$\displaystyle\exp\left(\alpha_{d}|\ln H|\right)$ $\displaystyle=$
$\displaystyle H^{\alpha_{d}},$

with

$\displaystyle\alpha_{d}=\frac{c{\cal M}hH^{1-d}\beta_{d}^{-1}(H)}{|\ln
H|}=\left\\{\begin{array}[]{ll}({\cal M}\frac{h}{H})\frac{c}{|\ln
H|^{2}}\approxeq\frac{c}{|\ln H|^{2}},&d=2,\\\ ({\cal
M}\frac{h}{H})\frac{c}{|\ln H|}\approxeq\frac{c}{|\ln
H|},&d=3.\end{array}\right.$

Since $\beta_{d}(H)=H^{-1}$ when $d=3$ and $\beta_{d}(H)=|\ln H|$ when $d=2$
(defined in Lemma 2), we can easily check that $chH^{1-d}\beta_{d}^{-1}(H)\leq
1$ holds when $h\leq H/c$ for $d=3$ or $h\leq H|\ln H|/c$ for $d=2$, which is
a very trivial condition.

Similar with the estimate of
$\|\nabla{\boldsymbol{e}}_{{\boldsymbol{u}},h}^{j}\|_{\Omega}+\|e_{p,h}^{j}\|_{\Omega}$
in Lemma 4, we can derive an estimate for the approximate ”local” error of the
coarse mesh Galerkin approximation, say
$\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}+\|\hat{p}_{H}^{j}\|_{\Omega}$
in the following.

######  Lemma 5.

Under the assumptions of A1, A2 and A3, we have for $j=1,2,\cdots,N$

$\displaystyle\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}+\|\hat{p}_{H}^{j}\|_{\Omega}\lesssim\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H})\|_{D_{j}}+\|p-p_{H}\|_{D_{j}}.$
(25)

###### Proof.

First we rewrite the problem (16) as

$\displaystyle
B([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},\phi_{j}{\boldsymbol{v}}_{h})-a({\boldsymbol{u}}_{H},\phi_{j}{\boldsymbol{v}}_{h})$
$\displaystyle\qquad-d(p_{H},\phi_{j}{\boldsymbol{v}}_{h})+d(\phi_{j}q_{h},{\boldsymbol{u}}_{H}),\quad\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times
Q_{h},$

or equivalently,

$B([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}],[{\boldsymbol{v}}_{h},q_{h}])=a({\boldsymbol{u}}-{\boldsymbol{u}}_{H},\phi_{j}{\boldsymbol{v}}_{h})+d(p-p_{H},\phi_{j}{\boldsymbol{v}}_{h})-d(\phi_{j}q_{h},{\boldsymbol{u}}-{\boldsymbol{u}}_{H}).$

Since ${\boldsymbol{X}}_{0}^{h}\times Q^{h}$ is a stable pair of the finite
element spaces, we have

$\displaystyle\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}+\|\hat{p}_{H}^{j}\|_{\Omega}\lesssim\sup\limits_{[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times
M_{h}}\frac{B([\hat{\boldsymbol{u}}_{H}^{j},\hat{p}_{H}^{j}],[{\boldsymbol{v}}_{h},q_{h}])}{\|\nabla{\boldsymbol{v}}_{h}\|_{\Omega}+\|q_{h}\|_{\Omega}}$
$\displaystyle\qquad=\sup\limits_{[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times
M_{h}}\frac{a({\boldsymbol{u}}-{\boldsymbol{u}}_{H},\phi_{j}{\boldsymbol{v}}_{h})+d(p-p_{H},\phi_{j}{\boldsymbol{v}}_{h})-d(\phi_{j}q_{h},{\boldsymbol{u}}-{\boldsymbol{u}}_{H})}{\|\nabla{\boldsymbol{v}}_{h}\|_{\Omega}+\|q_{h}\|_{\Omega}}$
$\displaystyle\qquad\lesssim\sup\limits_{{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}}\frac{a({\boldsymbol{u}}-{\boldsymbol{u}}_{H},\phi_{j}{\boldsymbol{v}}_{h})}{\|\nabla{\boldsymbol{v}}_{h}\|_{\Omega}}+\sup\limits_{{\boldsymbol{v}}_{h}\in{\boldsymbol{X}}_{0}^{h}}\frac{d(p-p_{H},\phi_{j}{\boldsymbol{v}}_{h})}{\|\nabla{\boldsymbol{v}}_{h}\|_{\Omega}}+\sup\limits_{q_{h}\in
M^{h}}\frac{d(\phi_{j}q_{h},{\boldsymbol{u}}-{\boldsymbol{u}}_{H})}{\|q_{h}\|_{\Omega}}.$

Note that $I_{H}$ is the Lagrange finite element interpolation operator of
${\boldsymbol{X}}$ onto ${\boldsymbol{X}}^{H}$ and $\hat{I}_{H}=I-I_{H}$,
especially, $\mbox{supp }\phi_{j}=D_{j}$, and
$a({\boldsymbol{u}}-{\boldsymbol{u}}_{H},\phi_{j}{\boldsymbol{v}}_{h})+d(p-p_{H},\phi_{j}{\boldsymbol{v}}_{h})=a({\boldsymbol{u}}-{\boldsymbol{u}}_{H},\hat{I}_{H}[\phi_{j}{\boldsymbol{v}}_{h}])+d(p-p_{H},\hat{I}_{H}\phi_{j}{\boldsymbol{v}}_{h})$,
it’s easy to derive

$\displaystyle
a({\boldsymbol{u}}-{\boldsymbol{u}}_{H},\hat{I}_{H}[\phi_{j}{\boldsymbol{v}}_{h}])=a({\boldsymbol{u}}-{\boldsymbol{u}}_{H},\hat{I}_{H}[\phi_{j}I_{H}{\boldsymbol{v}}_{h}]+\hat{I}_{H}[\phi_{j}\hat{I}_{H}{\boldsymbol{v}}_{h}])$
$\displaystyle\qquad\lesssim\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H})\|_{D_{j}}(\|\nabla[\hat{I}_{H}(\phi_{j}I_{H}{\boldsymbol{v}}_{h})]\|_{D_{j}}+\|\nabla\hat{I}_{H}[\phi_{j}\hat{I}_{H}{\boldsymbol{v}}_{h}]\|_{D_{j}}),$
$\displaystyle
d(p-p_{H},\hat{I}_{H}\phi_{j}{\boldsymbol{v}}_{h})=d(p-p_{H},\hat{I}_{H}[\phi_{j}I_{H}{\boldsymbol{v}}_{h}]+\hat{I}_{H}[\phi_{j}\hat{I}_{H}{\boldsymbol{v}}_{h}])$
$\displaystyle\qquad\lesssim\|p-p_{H}\|_{D_{j}}(\|\nabla[\hat{I}_{H}(\phi_{j}I_{H}{\boldsymbol{v}}_{h})]\|_{D_{j}}+\|\nabla\hat{I}_{H}[\phi_{j}\hat{I}_{H}{\boldsymbol{v}}_{h}]\|_{D_{j}}),$
$\displaystyle
d(\phi_{j}q_{h},{\boldsymbol{u}}-{\boldsymbol{u}}_{H})\lesssim\|q_{h}\|_{D_{j}}\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H})\|_{D_{j}}.$

In view of A1, A2 and noting that $\phi_{j}$ is a linear function on each mesh
of $\tau_{\Omega}^{H}$ with $|D\phi_{j}|\lesssim H^{-1}$,

$\displaystyle\|\nabla[\hat{I}_{H}(\phi_{j}I_{H}{\boldsymbol{v}}_{h})]\|_{D_{j}}=(\sum\limits_{\tau_{\Omega}^{H}\subset
D_{j}}\|\nabla\hat{I}_{H}(\phi_{j}I_{H}{\boldsymbol{v}}_{h})\|_{\tau_{\Omega}^{H}}^{2})^{\frac{1}{2}}\leq
H(\sum\limits_{\tau_{\Omega}^{H}\subset
D_{j}}\|D^{2}(\phi_{j}I_{H}{\boldsymbol{v}}_{h})\|_{\tau_{\Omega}^{H}}^{2})^{\frac{1}{2}}$
$\displaystyle\qquad\lesssim H(\sum\limits_{\tau_{\Omega}^{H}\subset
D_{j}}[\|\phi_{j}D^{2}(I_{H}{\boldsymbol{v}}_{h})\|_{\tau_{\Omega}^{H}}^{2}+\|D\phi_{j}D(I_{H}{\boldsymbol{v}}_{h})\|_{\tau_{\Omega}^{H}}^{2}])^{\frac{1}{2}}$
$\displaystyle\qquad\lesssim H(\sum\limits_{\tau_{\Omega}^{H}\subset
D_{j}}[H^{-2}\|D(I_{H}{\boldsymbol{v}}_{h})\|_{\tau_{\Omega}^{H}}^{2}+H^{-2}\|D(I_{H}{\boldsymbol{v}}_{h})\|_{\tau_{\Omega}^{H}}^{2}])^{\frac{1}{2}}\lesssim\|\nabla{\boldsymbol{v}}_{h}\|_{D_{j}},$
$\displaystyle\|\nabla\hat{I}_{H}[\phi_{j}\hat{I}_{H}{\boldsymbol{v}}_{h}]\|_{D_{j}}\lesssim\|D(\phi_{j}\hat{I}_{H}{\boldsymbol{v}}_{h})\|_{D_{j}}\lesssim\|D\phi_{j}\hat{I}_{H}{\boldsymbol{v}}_{h}\|_{D_{j}}+\|\phi_{j}D(\hat{I}_{H}{\boldsymbol{v}}_{h})\|_{D_{j}}\lesssim\|\nabla{\boldsymbol{v}}_{h}\|_{D_{j}}.$

Combining these estimates immediately yield

$\|\nabla\hat{\boldsymbol{u}}_{H}^{j}\|_{\Omega}+\|\hat{p}_{H}^{j}\|_{\Omega}\lesssim\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H})\|_{D_{j}}+\|p-p_{H}\|_{D_{j}}.$

The proof is complete. ∎

Now we can state the main theorem concerning the approximate accuracy of
solutions derived by the scheme (9)-(12) as follows.

######  Theorem 6.

Assume that A1, A2, A3 and (4) hold and
$[{\boldsymbol{u}},p]\in{\boldsymbol{H}}^{r+1}(\Omega)\times{\boldsymbol{H}}^{r}(\Omega)$.
Then we have

$\displaystyle\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p_{h}-p_{H}^{h}\|_{\Omega}$
$\displaystyle\lesssim
H^{\alpha_{d}}(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H})\|_{\Omega}+\|p_{h}-p_{H}\|_{\Omega}),$
(26) $\displaystyle\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega}$
$\displaystyle\lesssim
H(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p_{h}-p_{H}^{h}\|_{\Omega}),$
(27)
$\displaystyle\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p-p_{H}^{h}\|_{\Omega}$
$\displaystyle\lesssim
h^{r}+H^{\alpha_{d}}(\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H})\|_{\Omega}+\|p-p_{H}\|_{\Omega}),$
(28) $\displaystyle\|{\boldsymbol{u}}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega}$
$\displaystyle\lesssim
h^{r+1}+H(\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p-p_{H}^{h}\|_{\Omega}),$
(29)

where $\alpha_{d}>0$ is defined in Lemma 4.

###### Proof.

We know from (11) that

$\|\nabla{\boldsymbol{E}}_{\boldsymbol{u}}^{H}\|_{\Omega}+\|E_{p}^{H}\|_{\Omega}\lesssim\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H,h})\|_{\Omega}+\|p_{h}-p_{H,h}\|_{\Omega}.$

Since

$[{\boldsymbol{u}}-{\boldsymbol{u}}_{H}^{h},p-p_{H}^{h}]=[{\boldsymbol{u}}-{\boldsymbol{u}}_{h},p-p_{h}]+[{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H,h},p_{h}-p_{H,h}]+[{\boldsymbol{E}}_{\boldsymbol{u}}^{H},E_{p}^{H}],$

and

$\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{h})\|_{\Omega}+\|p-p_{h}\|_{\Omega}\lesssim
h^{r},$

the conclusions (26) and (28) follows immediately due to Lemma 4 and Lemma 5
and the triangle inequality.

To estimate the $L^{2}$-norms of the velocity, we firstly introduce a Stokes
projection $P_{H}=[P_{H}^{\boldsymbol{u}},P_{H}^{p}]$ defined from
${\boldsymbol{X}}_{0}\times Q$ onto ${\boldsymbol{X}}_{0}^{H}\times Q^{H}$:
for given $[{\boldsymbol{w}},r]\in{\boldsymbol{X}}_{0}\times Q$, find
$[P_{H}^{\boldsymbol{u}}{\boldsymbol{w}},P_{H}^{p}r]\in{\boldsymbol{X}}_{0}^{H}\times
Q^{H}$ such that

$B([{\boldsymbol{v}}_{H},q_{H}],[{\boldsymbol{w}}-P_{H}^{\boldsymbol{u}}{\boldsymbol{w}},r-P_{H}^{p}r])=0,\quad\forall[{\boldsymbol{v}}_{H},q_{H}]\in{\boldsymbol{X}}_{0}^{H}\times
Q^{H}.$

In view of the definitions of $\hat{\boldsymbol{u}}_{H,h}$ and
$\hat{p}_{H,h}$, and $\sum\limits_{j=1}^{N}\phi_{j}=1$ in $\Omega$, adding up
all the equations of (13) over $j$ with $\mu=0$ gives

$\displaystyle
B([\hat{\boldsymbol{u}}_{H,h},\hat{p}_{H,h}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},{\boldsymbol{v}}_{h})-B([{\boldsymbol{u}}_{H},p_{H}],[{\boldsymbol{v}}_{h},q_{h}])$
(30)
$\displaystyle\qquad-\sum\limits_{j=1}^{N}\int_{\Gamma_{j}}{\boldsymbol{\xi}}^{j}{\boldsymbol{v}}_{h}ds,\quad\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times
Q^{h}.$

Since
$[{\boldsymbol{u}}_{H,h},p_{H,h}]=[{\boldsymbol{u}}_{H}+\hat{\boldsymbol{u}}_{H,h},p_{H}+\hat{p}_{H,h}]$,
then we have

$B([{\boldsymbol{u}}_{H,h},p_{H,h}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},{\boldsymbol{v}}_{h})-\sum\limits_{j=1}^{N}\int_{\Gamma_{j}}{\boldsymbol{\xi}}^{j}{\boldsymbol{v}}_{h}ds,\quad\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times
Q^{h}.$

On the another aspect, by the definition of the Stokes projection $P_{H}$
above, we can rewrite (11) associated with the coarse mesh correction as

$\displaystyle
B([{\boldsymbol{E}}_{\boldsymbol{u}}^{H},E_{p}^{H}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},P_{H}^{\boldsymbol{u}}{\boldsymbol{v}}_{h})-B([{\boldsymbol{u}}_{H,h},p_{H,h}],[P_{H}^{\boldsymbol{u}}{\boldsymbol{v}}_{h},P_{H}^{p}q_{h}]).$

Combining these two equations leads to

$\displaystyle
B([{\boldsymbol{u}}_{H}^{h},p_{H}^{h}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},{\boldsymbol{v}}_{h})+({\boldsymbol{f}},P_{H}^{\boldsymbol{u}}{\boldsymbol{v}}_{h})-B([{\boldsymbol{u}}_{H,h},p_{H,h}],[P_{H}^{\boldsymbol{u}}{\boldsymbol{v}}_{h},P_{H}^{p}q_{h}])$
$\displaystyle\qquad\quad-\sum\limits_{j=1}^{N}\int_{\Gamma_{j}}{\boldsymbol{\xi}}^{j}{\boldsymbol{v}}_{h}ds,\quad\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times
Q^{h}.$

Finally, we obtain
$\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times Q^{h}$

$B([{\boldsymbol{u}}_{H}^{h},p_{H}^{h}],[{\boldsymbol{v}}_{h},q_{h}])=({\boldsymbol{f}},{\boldsymbol{v}}_{h})-\sum\limits_{j=1}^{N}\int_{\Gamma_{j}}{\boldsymbol{\xi}}^{j}(I-P_{H}^{\boldsymbol{u}}){\boldsymbol{v}}_{h}ds,$

and $\forall[{\boldsymbol{v}}_{h},q_{h}]\in{\boldsymbol{X}}_{0}^{h}\times
Q^{h}$

$B([{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h},p_{h}-p_{H}^{h}],[{\boldsymbol{v}}_{h},q_{h}])=\sum\limits_{j=1}^{N}\int_{\Gamma_{j}}{\boldsymbol{\xi}}^{j}(I-P_{H}^{\boldsymbol{u}}){\boldsymbol{v}}_{h}ds.$
(31)

From this error equation about
$[{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h},p_{h}-p_{H}^{h}]$, and by
virtue of Lemma 1, 4 and 5, we can easily get

$\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p_{h}-p_{H}^{h}\|_{\Omega}\lesssim
H^{\alpha_{d}}(\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H})\|_{\Omega}+\|p-p_{H}\|_{\Omega}).$

Then we can derive the first result by using the triangle inequality.

Next, we also use the Aubin-Nitsche duality argument to show $L^{2}-$error
estimate of the velocity. By A3, for
${\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h}\in{\boldsymbol{L}}^{2}(\Omega)$,
there exists $[{\boldsymbol{w}},r]\in{\boldsymbol{H}}^{2}(\Omega)\times
H^{1}(\Omega)$ safistying

$B([{\boldsymbol{v}},q],[{\boldsymbol{w}},r])=({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h},{\boldsymbol{v}}),\quad\forall[{\boldsymbol{v}},q]\in{\boldsymbol{X}}_{0}(\Omega)\times
Q,$

and

$\|{\boldsymbol{w}}\|_{2,\Omega}+\|r\|_{1,\Omega}\lesssim\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega}.$

If taking
$[{\boldsymbol{v}},q]=[{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h},p_{h}-p_{H}^{h}]$
here and in view of (31), we have

$\displaystyle\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega}^{2}=B([{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h},p_{h}-p_{H}^{h}],[{\boldsymbol{w}},r])$
$\displaystyle\qquad=B([{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h},p_{h}-p_{H}^{h}],[{\boldsymbol{w}}-P_{H}^{\boldsymbol{u}}{\boldsymbol{w}},r-P_{H}^{p}r]).$

Then we can proceeds the estimate as

$\displaystyle\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega}^{2}\lesssim(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p_{h}-p_{H}^{h}\|_{\Omega})\|\nabla(I-P_{H}^{\boldsymbol{u}}){\boldsymbol{w}}\|_{\Omega}$
$\displaystyle\qquad\quad+\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}\|(I-P_{H}^{p})r\|_{\Omega}$
$\displaystyle\quad\lesssim
H(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p_{h}-p_{H}^{h}\|_{\Omega})(\|{\boldsymbol{w}}\|_{2,\Omega}+\|r\|_{1,\Omega})$
$\displaystyle\quad\lesssim
H(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p_{h}-p_{H}^{h}\|_{\Omega})\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega},$

which immediately yields

$\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega}\lesssim
H(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p_{h}-p_{H}^{h}\|_{\Omega}).$

By using the triangle inequality again, we prove $L^{2}-$error estimate for
the velocity. ∎

######  Remark 4.2.

By the results in Theorem 6 above, especially (26) and (27), we can clearly
improve the convergence orders by introducing a two-grid iteration. Actually,
we can verify that all the lemmas and theorems above still hold when replacing
$[{\boldsymbol{u}}_{H},p_{H}]$ by $[{\boldsymbol{u}}_{H}^{h},p_{H}^{h}]$. This
suggests us to introducing the following two-grid iterative scheme with a
iteration number of K.

$K=[\alpha_{d}^{-1}+0.5]=\left\\{\begin{array}[]{ll}O(|\ln H|^{2}),&d=2,\\\
O(|\ln H|),&d=3.\end{array}\right.$ (32)

Local and parallel two-grid iterative scheme:

Step 0.

Setting $k=0$; Deriving $[{\boldsymbol{u}}_{H},p_{H}]$ by solving (3) and
denoting by
$[{\boldsymbol{u}}_{H}^{0,h},p_{H}^{0,h}]=[{\boldsymbol{u}}_{H},p_{H}]$;

Step 1.

For $k\geq 0$, solving the equations (9) with
$[{\boldsymbol{u}}_{H},p_{H}]=[{\boldsymbol{u}}_{H}^{k,h},p_{H}^{k,h}]$ to get
$\\{[\hat{\boldsymbol{u}}_{H,h}^{j},\hat{p}_{H,h}^{j}]\\}_{j=1}^{N}$ for each
$j$, which are denoted as
$\\{[\hat{\boldsymbol{u}}_{H,h}^{k+1,j},\hat{p}_{H,h}^{k+1,j}]\\}_{j=1}^{N}$.
Then we construct $[{\boldsymbol{u}}_{H,h}^{k+1},p_{H,h}^{k+1}]$ by formula
(10);

Step 2.

Deriving $[{\boldsymbol{E}}_{\boldsymbol{u}}^{k+1},E_{p}^{k+1}]$ by solving
(11) with
$[{\boldsymbol{u}}_{H,h},p_{H,h}]=[{\boldsymbol{u}}_{H,h}^{k+1},p_{H,h}^{k+1}]$,
and constructing

$[{\boldsymbol{u}}_{H}^{k+1,h},p_{H}^{k+1,h}]=[{\boldsymbol{u}}_{H,h}^{k+1}+{\boldsymbol{E}}_{\boldsymbol{u}}^{k+1},p_{H,h}^{k+1}+E_{p}^{k+1}].$

Step 3.

Checking whether $k+1>K$, if yes, terminating the iteration and deriving
$[{\boldsymbol{u}}_{H}^{h},p_{H}^{h}]=[{\boldsymbol{u}}_{H}^{k+1,h},p_{H}^{k+1,h}]$;
otherwise, letting $k:=k+1$ and turning to Step 1.

######  Theorem 7.

Suppose $[{\boldsymbol{u}},p]\in{\boldsymbol{H}}^{r+1}(\Omega)\times
H^{r}(\Omega)$, the final approximation $[{\boldsymbol{u}}_{H}^{h},p_{H}^{h}]$
for the iterative scheme with $K$ defined by (32) has the following error
bounds

$\displaystyle\|\nabla({\boldsymbol{u}}-{\boldsymbol{u}}_{H}^{h})\|_{\Omega}+\|p-p_{H}^{h}\|_{\Omega}\lesssim
h^{r}+H^{r+1},$ (33)
$\displaystyle\|{\boldsymbol{u}}-{\boldsymbol{u}}_{H}^{h}\|_{\Omega}\lesssim
h^{r+1}+H^{r+2}.$ (34)

###### Proof.

Note that, the Ritz-Projection of $[{\boldsymbol{u}}_{H}^{h},p_{H}^{h}]$ onto
the corresponding course grid finite element spaces is nothing but the
standard Galerkin approximation, see (31), when selecting
$[{\boldsymbol{v}}_{h},q_{h}]=[{\boldsymbol{v}}_{H},q_{H}]$, we have

$\displaystyle
B([{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{h},p_{h}-p_{H}^{h}],[{\boldsymbol{v}}_{H},q_{H}])=0,$

then,

$\displaystyle
B([{\boldsymbol{u}}_{H}^{h},p_{H}^{h}],[{\boldsymbol{v}}_{H},q_{H}])=B([{\boldsymbol{u}}_{h},p_{h}],[{\boldsymbol{v}}_{H},q_{H}])=({\boldsymbol{f}},{\boldsymbol{v}}_{H})=B([{\boldsymbol{u}}_{H},p_{H}],[{\boldsymbol{v}}_{H},q_{H}]).$

Then, we can update $[{\boldsymbol{u}}_{H}^{0,h},p_{H}^{0,h}]$ by
$[{\boldsymbol{u}}_{H}^{1,h},p_{H}^{1,h}]$, and execute the previous local and
parallel two-grid iterative scheme (without Step 0) again to get
$[{\boldsymbol{u}}_{H}^{2,h},p_{H}^{2,h}]$, then follow (26) and (27) to
derive the error estimates of the form, equivalently, there exists a mesh
independent constant $c>0$,

$\displaystyle\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{2,h})\|_{\Omega}+\|p_{h}-p_{H}^{2,h}\|_{\Omega}$
$\displaystyle\lesssim
H^{\alpha_{d}}(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{1,h})\|_{\Omega}+\|p_{h}-p_{H}^{1,h}\|_{\Omega})$
$\displaystyle\lesssim
H^{2\alpha_{d}}(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{0,h})\|_{\Omega}+\|p_{h}-p_{H}^{0,h}\|_{\Omega}),$
$\displaystyle\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{2,h}\|_{\Omega}$
$\displaystyle\lesssim
H(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{2,h})\|_{\Omega}+\|p_{h}-p_{H}^{2,h}\|_{\Omega}).$

Therefore, cyclically executing our two-grid iterative scheme with updating
data by a iteration of $K$ times, we can finally derive

$\displaystyle\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{K,h})\|_{\Omega}+\|p_{h}-p_{H}^{K,h}\|_{\Omega}$
$\displaystyle\lesssim
H^{K\alpha_{d}}(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{0,h})\|_{\Omega}+\|p_{h}-p_{H}^{0,h}\|_{\Omega})$
$\displaystyle\lesssim
H^{K\alpha_{d}}(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H})\|_{\Omega}+\|p_{h}-p_{H}\|_{\Omega}),$
$\displaystyle\|{\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{K,h}\|_{\Omega}$
$\displaystyle\lesssim
H(\|\nabla({\boldsymbol{u}}_{h}-{\boldsymbol{u}}_{H}^{K,h})\|_{\Omega}+\|p_{h}-p_{H}^{K,h}\|_{\Omega}).$

By the choice of $K$ within (32) and the triangle inequality, since
$[{\boldsymbol{u}}_{H}^{h},p_{H}^{h}]=[{\boldsymbol{u}}_{H}^{K,h},p_{H}^{K,h}]$
now, we finally get (33) and (34).

∎

######  Remark 4.3.

By the estimates in Theorem 7, we can clearly know that in order to get the
optimal convergence orders of the velocity - pressure pair in the sense of
${\boldsymbol{H}}^{1}-L^{2}$-norms, or that of the velocity in
${\boldsymbol{L}}^{2}$-norm, we should choose $H$ and $h$ such that

$h\sim H^{\frac{r+1}{r}}\quad\mbox{or}\quad h\sim H^{\frac{r+2}{r+1}},$ (35)

respectively.

##  5 Numerical Tests

In all the numerical experiments below, the algorithms are implemented using
public domain finite element software Freefem++[27]. All simulations were
performed on a Dawning parallel cluster composed of 64 nodes (each node
consists of eight-core 2.0 GHz CPU, 8$\times$2 GB DRAM, and all the nodes are
connected via 20Gbps InfiniBand). The message-passing interface is supported
by MPICH.

We introduce the following notations for convenience:

SFEM means the standard finite element method.

EPLP denotes the expandable local and parallel two-grid finite element
iterative scheme.

Wall time covers the maximal CPU time among all processors used for EPLP,
including the CPU time for solving the two global coarse grid problems and the
parallel time for solving all the subproblems.

###  5.1 Problem 1

Firstly, to verify the theoretical results, we consider the following 2D
numerical example (referred as Problem 1) with exact solution

$\displaystyle\mathbf{u}=(10x^{2}(x-1)^{2}y(y-1)(2y-1),-10x(x-1)(2x-1)y^{2}(y-1)^{2}),$
$\displaystyle p=10(2x-1)(2y-1).$

The domain is selected as the unit square $\Omega=[0,1]\times[0,1]$ with a
uniform triangulation $T^{H}$. The Taylor-Hood finite element pairs are used
in solving the Stokes equation.

Fig. 3:  Structured mesh for Problem 1, $D_{j}$: local domain (red) and
$\Omega_{j}$: expanded domain of $D_{j}$ with one mesh layer (both red and
green).

According to Theorem 7 for EPLP, we know that the estimate has the form of

$\displaystyle\|\nabla(\mathbf{u}-{\mathbf{u}}_{H}^{h})\|_{\Omega}+\|p-{p}^{h}_{H}\|_{\Omega}=O(h^{2}+H^{3}),$

which suggests us to select $H$ and $h$ such that $h\approxeq H^{\frac{3}{2}}$
when verifying the optimal convergence orders for $H^{1}$-norm of velocity and
$L^{2}$-norm of pressure.

To this end, for deriving the approximate $H^{1}$-error of velocity and
$L^{2}$-error of pressure, we implement EPLP in parallel with coarse meshes of
the sizes H=1/16, 1/25, 1/36, 1/49, 1/64 and corresponding fine meshes with
the sizes $h=1/64,\ 1/125,\ 1/216,\ 1/343,\ 1/512$, under the parallel
environment of fixed 64 processors. Meanwhile, SFEM is executed for the same
fine mesh on a single processor. The results for both methods are presented in
Table 2, 1, respectively.

Table 1: The errors of SFEM for Problem 1.

$h$ $||\mathbf{u}-\mathbf{u}_{h}||_{1,\Omega}$ Order $||p-p_{h}||_{0,\Omega}$
Order CPU $1/64$ 0.000205741 0.000630369 2.718 $1/125$ 5.396e-05 1.99927
0.000165247 2.00000 68.49 $1/216$ - - - - - $1/343$ - - - - - $1/512$ - - - -
-

Table 2: The errors of EPLP for Problem 1, $h\approxeq H^{\frac{3}{2}}$, 64
processors.

$H$ $h$ $K$ $||\mathbf{u}-{\mathbf{u}}_{H}^{h}||_{1,\Omega}$ Order
$||p-{p}_{H}^{h}||_{0,\Omega}$ Order Wall time $1/16$ $1/64$ 2 0.000299829
0.00077599 10.54 $1/25$ $1/125$ 3 5.44032e-05 2.54961 0.000165716 2.30623
91.51 $1/36$ $1/216$ 3 1.81385e-05 2.00815 5.53793e-05 2.00391 465.08 $1/49$
$1/343$ 3 7.19776e-06 1.99862 2.19737e-05 1.99882 2084.38 $1/64$ $1/512$ 4
3.25815e-06 1.97859 1.01485e-05 1.92844 7892.11

Table 3: The $L^{2}$-error of Velocity by EPLP for Problem 1, $h\approxeq
H^{\frac{4}{3}}$, 64 processors.

$H$ $h$ $K$ $||\mathbf{u}-{\mathbf{u}}_{H}^{h}||_{0,\Omega}$ Order Wall time
$1/32$ $1/96$ 3 1.0881e-07 80.44 $1/64$ $1/256$ 4 5.71659e-09 3.00382 2032.04
$1/96$ $1/384$ 4 1.74889e-09 2.92107 10435.2

From Table 2, the optimal orders for the $H^{1}$-error of velocity and
$L^{2}$-error of pressure by EPLP are observed, which can verify the
theoretical result. Also the approximate accuracy by EPLP are comparable with
those of SFEM under same fine meshes in Table 1. Noting that, in the present
computational environment, SFEM fails for mesh with$h=1/216$ or finer,
meanwhile our EPLP still works well for much finer meshes including $h=1/512$.

Similarly, by Theorem 7, we choose mesh pairs with a relation of $h\approxeq
H^{\frac{4}{3}}$ to check the order for $L^{2}$-error of the velocity computed
by our EPLP. The parameters and computational results are shown in Table 3,
from which the optimal order of $O(h^{3})$ is verified. This also supports the
result in the theoretical analysis.

Table 4: Wall time $T(J)$ in seconds, speedup $S_{p}$ and parallel efficiency
$E_{p}$ for Example 1, with $H=1/36$.

$J$ 2 4 8 16 32 64 $T(J)$ 9204.77 4684.47 2520.25 1392.31 747.51 465.08
$S_{p}=\frac{T(2)}{T(J)}$ 1.00 1.96 3.65 6.61 12.31 19.79 $E_{p}=\frac{2\times
T(2)}{J\times T(J)}$ 1.00 0.98 0.91 0.83 0.77 0.62

The performance of a parallel algorithm in a homogeneous parallel environment
is usually measured by speedup and parallel efficiency, commonly defined as

$\displaystyle S_{p}=\frac{T(J_{1})}{T(J_{2})},\ E_{p}=\frac{J_{1}\times
T(J_{1})}{J_{2}\times T(J_{2})},$

where $T(J_{1})$ and $T(J_{2})$ ($J_{1}\leq J_{2}$) are wall times of the
parallel program when using $J_{1}$ and $J_{2}$ processors, respectively.

Table 4 reports the wall time of EPLP in a parallel environment using
processors of number $J=2,\ 4,\ 8,\ 16,\ 32,\ 64$, and presents the
corresponding speedup and parallel efficiency, which are computed by
comparison with $J_{1}=2$. These results show good parallel performance of our
EPLP.

###  5.2 Problem 2

Fig. 4:  A regular unstructured mesh for Problem 2, $D_{j}$: a local domain
(red) and $\Omega_{j}$: an expanded domain of $D_{j}$ (both red and green).
Table 5: The errors of EPLP for Problem 2 with Taylor-Hood elements,
$h\approxeq H^{\frac{3}{2}}$, with 64 processors.

$H$ $h$ $K$ $||\mathbf{u}-{\mathbf{u}}^{h}||_{1,\Omega}$ Order
$||p-{p}^{h}||_{0,\Omega}$ Order Wall time $1/16$ $1/64$ 2 0.00237151
9.76556e-05 10.75 $1/32$ $1/160$ 3 0.000415357 1.90130 1.51911e-05 2.03073
196.61 $1/48$ $1/288$ 3 0.000137697 1.87837 4.88635e-06 1.92972 1137.06 $1/64$
$1/448$ 4 4.91655e-05 2.33089 2.05054e-06 1.96532 6254.84

Table 6: The errors of EPLP for Problem 2 with Mini-elements, $h\approxeq
H^{2}$, with 64 processors.

$H$ $h$ $K$ $||\mathbf{u}-{\mathbf{u}}^{h}||_{1,\Omega}$ Order
$||p-{p}^{h}||_{0,\Omega}$ Order Wall time $1/8$ $1/64$ 2 0.151405 0.0550405
14.67 $1/16$ $1/256$ 2 0.0361013 1.03415 0.0131228 1.03421 538.43 $1/32$
$1/1024$ 3 0.00924041 0.983011 0.00212624 1.31285 52038.4

To further test our EPLP, we also consider another smooth problem (referred as
Problem 2) with exact solution of the trigonometric form

$\displaystyle\mathbf{u}=(\sin(\pi x)^{2}\sin(2\pi y),-\sin(2\pi x)\sin(\pi
y)^{2}),$ $\displaystyle p=\cos(\pi x)\cos(\pi y),$

on the domain $\Omega=[0,1]\times[0,1]$. The regular unstructured
triangulations, and Taylor-Hood elements and Mini-elements are used for our
EPLP. The corresponding approximate results are listed in Table 5 and 6,
respectively. From these two tables, one can observe that, for EPLP with $K$
iterations, the $H^{1}$-error of velocity and $L^{2}$-error of pressure
achieve the optimal orders, namely, $O(h^{2})$ by Taylor-Hood elements, and
$O(h)$ by Mini-elements, which support theoretical results in Theorem 7.

Table 7: The errors of EPLP for Problem 3 with Taylor-Hood element pair,
$h\approxeq H^{\frac{3}{2}}$, with 64 processors.

$H$ $h$ $K$ $||\mathbf{u}-{\mathbf{u}}^{h}||_{1,\Omega}$ Order
$||p-{p}^{h}||_{0,\Omega}$ Order $1/4$ $1/8$ 1 0.763 0.0408459 $1/8$ $1/24$ 1
0.104595 1.80879 0.00564176 1.80192 $1/16$ $1/64$ 2 0.0154352 1.95084
0.000769727 2.03085

###  5.3 Problem 3

As the final experiment, we will test our EPLP using a three dimensional
problem (referred as Problem 3) with $\nu=0.05$ and exact solution

$\displaystyle\mathbf{u}=\left(\sin(\pi x)^{2}\sin(2\pi y)\sin(2\pi
z),-\sin(2\pi x)\sin(\pi y)^{2}\sin(2\pi z),\sin(2\pi x)\sin(2\pi y)\right),$
$\displaystyle p=\cos(2\pi x)\cos(2\pi y)\cos(2\pi z).$

The uniform triangulation and Taylor-Hood elements are used for $T^{H}$ with
mesh size $H$ and the computational domain is chosen as the unit cube
$\Omega=[0,1]\times[0,1]\times[0,1]$. The computational results are shown in
Table 7, from which one can observe that, while $h$ deceases, by suitable
iteration of $K$, both the $H^{1}$-error of velocity and $L^{2}$-error of
pressure can reach the optimal orders of $O(h^{2})$, which also support our
theoretical analysis.

##  6 Conclusions

In this paper, we have designed an expandable local and parallel two-grid
finite element iterative scheme based on superposition principle for the
Stokes problem. The optimal convergence orders of the scheme are analyzed and
obtained within suitable two-grid iterations while numerical tests in 2D and
3D are carried out to show the flexible and high efficiency of the scheme. The
extension of the scheme to the time-dependent problems or nonlinear problems,
e.g., Navier-Stokes equations, will be our further work.

## References

  * [1] J Xu, Two-Grid Discretization Techniques for Linear and Nonlinear PDEs, SIAM J. Numer. Anal., 33 (5) (1996), 1759-1777
  * [2] A. Ait Ou Ammi and M. Marion, _Nonlinear Galerkin Methods and Mixed Finite Elements: Two-Grid Algorithms for the Navier-Stokes Equations_ , Numer. Math., 68(1994), 189–213.
  * [3] T. Chan and T. Mathew, _Domain Decomposition Algorithms_ , Acta Numer., 3(1994), 61–143.
  * [4] W. Layton and L. Tobiska,_A Two-Level Method with Backtracking for the Navier-Stokes Equations_ , SIAM J. Numer. Anal., 35 (5)(1998), 2035–2054.
  * [5] Y. Hou and K. Li, A Small Eddy Correction Method for Nonlinear Dissipative Evolutionary Equations, SIAM J. Numer. Anal., 41 (3), 1101–1130
  * [6] Y. Hou and K. Li, Post-processing Fourier Galerkin Method for the Navier-Stokes Equations, SIAM J. Numer. Anal., 47:3(2009), 1909–1922.
  * [7] K. Li and Y. Hou, An AIM and One-step Newton Method for the Navier-Stokes Equations, Comput. Meth. Appl. Mech. Eng., 190(2001), 6141–6155.
  * [8] Q. Liu and Y. Hou, A Post-processing Mixed Finite Element Method for Navier-Stokes Equations, Int. J. Comput. Fluid Dyn., 23:6(2009), 461–475.
  * [9] J. Xu, _Iterative Methods by Space Decomposition and Subspace Correction_ , SIAM Rev., 34(1992), 581–613.
  * [10] J. Xu and A. Zhou, Local and Parallel Finite Element Algorithms based on Two-Grid Discretizations, Math. Comput., 69(2000), 881–909.
  * [11] R. E. Bank and M. Holst, A New Paradigm for Parallel Adaptive Meshing Algorithms, SIAM Rev., 45 (2003), 291–323.
  * [12] Y. He, J. Xu and A. Zhou, Local and Parallel Finite Element Algorithms for the Stokes Problem, Numer. Math. 109:3(2008), 415–434.
  * [13] J. Nitsche and A.H. Schatz, _Interior Estimates for Ritz-Galerkin Methods_ , Math. Comput., 28(1974), 937–955.
  * [14] I. Babuska and J.M. Melenk, _The Partition of Unity Method_ , Int. J. Numer. Meth. Eng., 40(1997), 727–758.
  * [15] C. Wang, Z. Huang and L. Li, _Two-Grid Partition of Unity Method for Second Order Elliptic Problems_ , Appl. Math. Mech., 29(2008), 527–533.
  * [16] H. Zheng, L. Song, Y. Hou and Y. Zhang, _The Partition of Unity Parallel Finite Element Algorithm_ , Adv. Comput. Math., 41(2015), 937–951.
  * [17] J. Yu, F. Shi and H. Zheng, Local and parallel finite element algorithms Based on the partition of unity for the Stokes problem, SIAM J. Sci. Comput., 36:5(2014), C547–C567.
  * [18] H. Zheng, J. Yu and F. Shi,Local and Parallel Finite Element Algorithm Based on the Partition of Unity for Incompressible Flows, J. Sci. Comput., 65 (2)(2015), 512–532.
  * [19] Y. Hou and G. Du, An Expandable Local and Parallel Two-Grid Finite Element Scheme, Comput. Math. App., 71:12(2016), 2541–2556.
  * [20] M.G. Larson and A. Målqvist, Adaptive Variational Multi-scale Methods Based on a Posteriori Error Estimation: Energy Norm Estimates for Elliptic Problems, Comput. Meth. Appl. Mech. Eng., 196(2007), 2313–2324.
  * [21] M.G. Larson and A. Målqvist, _An Adaptive Variational Multi-scale Method for Convection-Diffusion Problems_ , Comm. Numer. Meth. Eng., 25(2009), 65–79.
  * [22] R.A. Adams, _Sobolev Spaces_ , Academic Press, New York, 1975.
  * [23] P.G. Ciarlet and J.L. Lions, _Handbook of Numerical Analysis, Vol. II, Finite Element Methods (Part I)_ , North-Holland, Amsterdam, 1991.
  * [24] V. Girault and P.A. Raviart, _Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms_ , Springer-Verlag, Berlin Heidelberg New York Tokyo, 1986.
  * [25] V. Girault, R. Glowinski, H. López and J.-P. Vila, _A Boundary Multiplier/Fictitious Domain Method for the Steady Incompressible Navier-Stokes Equations_ , Numer. Math., 88(2001), 75–103.
  * [26] J.L. Lions and E. Magenes, Non-Homogeneous Boundary Value Problems and Applications, Springer-Verlag, Berlin Heidelberg New York, 1972.
  * [27] F. Hecht, New development in freefem++. J. Numer. Math., 20(2012), 251–265.