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

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 337, pp. 1--16.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2016/337\hfil Stokes problem in a domain with corners]
{Solving the Stokes problem in a domain with corners by the mortar
spectral element method}

\author[A. Al Salem, N. Chorfi \hfil EJDE-2016/337\hfilneg]
{Azhar Al Salem, Nejmeddine Chorfi}

\address{Azhar Al Salem \newline
Department of Mathematics,
College of Science,
King Saud University, P.O. Box 2455,
Riyadh 11451, Saudi Arabia}
\email{azhar2002@msn.com}

\address{Nejmeddine Chorfi \newline
Department of Mathematics,
College of Science,
King Saud University, P.O. Box 2455,
Riyadh 11451, Saudi Arabia}
\email{nchorfi@ksu.edu.sa}

\thanks{Submitted September 3, 2016. Published December 28, 2016.}
\subjclass[2010]{35J15, 78M22}
\keywords{Mortar spectral element method; Stokes problem; singularity;
\hfill\break\indent  Uzawa algorithm}

\begin{abstract}
 In this article, we implement the mortar spectral element method for the
 Stokes problem on a domain within corners.
 We consider the Strang and Fix algorithm, which permits to enlarge the
 discrete space of the velocity by the first singular function.
 The usefulness of this method is confirmed by the numerical
 results presented here.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks

\section{Introduction}

The solution of the Stokes equation in a domain of $\mathbb{R}^2$ with corners
is divided into a regular part and a linear combination of singular functions
 \cite{G,K}. To take into account these singularities, we propose to decompose
the domain.
The domain decomposition method consists in dividing the domain of resolution
into sub-domains of smaller sizes and simpler geometries.
The emergence of parallel computing and the development of effective codes
(Message Passing Interface) have motivated the use of the decomposition domain
methods. For these methods the matching conditions at the interfaces are
decisive for a good approximation of the solution. One of the most used method
is the \emph{Mortar Element Method} \cite{BMP}. The advantage of this method
is that it admits a weak matching condition type, which gives flexibility
 on the geometry and on the choice of the functional spaces.

The approximation of singular functions which appear on the decomposition
of the solution was done for the first time by Babu\v{s}ka et al \cite{BR,BS}
for the $p$-version of the finite element method. This work was adapted
by Bernardi et al \cite{BM1} for the spectral method.
They proved that the approximation order of these singular functions by
polynomials is better than that given by the general theory of the spectral
approximation \cite{BM2}. We use these approximation results to implement
the Strang and Fix \cite{SF} algorithm in the case of the Stokes problem.
This algorithm consists in enlarging the discrete space of the velocity by
the first singular function. We are interested in this work by the implementation
of the mortar spectral element method for the Stokes problem in a domain
with corners. We first deduce the matrix system which is solved using Uzawa
algorithm \cite{FG}. Next, we present some numerical results showing the interest
 of using the Strand and Fix algorithm for this type of problems.

An outline of this paper is as follows. In section 2, we present the geometry
of the domain and the continuous problem, we give the singular functions and
some regularity results. In section 3, we present the discrete problem and
the error result obtained from the discretization of the Stokes problem by
the mortar spectral method. Section 4 is devoted to the implementation of
the mortar spectral element method. We describe the matrix system and its
resolution algorithm. Finally, we present some numerical results which
confirm the interest of the method.

\section{Continuous problem and singular functions}

Let $\Omega$ be a polygonal domain of $\mathbb{R}^2$ simply connected and $\Gamma$
be its connected boundary. The generic point in $\Omega$ is denoted by
$\mathbf{x}=(x,y)$. We suppose that there exists a finite number of vertex
$\Gamma_j$ for $j\in\{1,\dots ,J\}$, $J$ is a positive integer such that
$$
\Gamma=\cup_{j=1}^{J} \Gamma_j.
$$
We consider $\mathbf{c}_j$ the corner of $\Omega$ made by $\Gamma_j$ and
$\Gamma_{j+1}$, and $\alpha_j$ its measure.
Let $H^s(\Omega)$ be the Sobolev space of order $s\in {\mathbb R}$,
$H_0^1(\Omega)$ the subspace of $H^1(\Omega)$ of functions which vanish
on the boundary $\Gamma$ and the space
$$
L^2_0(\Omega)=\{q\in L^2(\Omega): \int_\Omega q(\mathbf{x})\,d\mathbf{x}=0\}.
$$

We introduce the Stokes problem on the domain $\Omega$ in the velocity and
pressure formulation:


For a forcing term $\mathbf{f}$ in $[H^{-1}(\Omega)]^2$,
we try to find the velocity $\mathbf{u}\in [H^1_0(\Omega)]^2$ and the pressure
$p\in L^2_0(\Omega)$ such that
\begin{equation}
\begin{gathered}
 -\nu\Delta \mathbf{u}+\nabla p=\mathbf{f}\quad\text{in } \Omega\\
 \operatorname{div}\mathbf{u}=0 \quad\text{in } \Omega \\
 \mathbf{u} = 0 \quad\text{on } \Gamma,
\end{gathered} \label{e2.1}
\end{equation}
where $\nu$ is the viscosity of the fluid that is supposed to be positive and
constant. Problem \eqref{e2.1} has the following variational formulation:

For $\mathbf{f}$ in $[H^{-1}(\Omega)]^2$, find $\mathbf{u}$ in $[H^1_0(\Omega)]^2$
and $p$ in $L_0^2(\Omega)$ such that for all $\mathbf{v}$ in $[H^1_0(\Omega)]^2$
and for all $q$ in $L^2(\Omega)$:
\begin{equation}
\begin{gathered}
 a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)=\langle\mathbf{f},\mathbf{v}\rangle\\
b(\mathbf{u},q)=0,
\end{gathered}\label{e2.2}
\end{equation}
where
\begin{gather*}
 a(\mathbf{u},\mathbf{v})
=\nu\int_\Omega\nabla \mathbf{u}\nabla \mathbf{v} d\mathbf{x}, \\
b(\mathbf{u},q)=- \int_\Omega (\operatorname{div}\mathbf{u}) q d\mathbf{x}.
\end{gather*}
 We denote by $\langle\cdot,\cdot\rangle$ the duality pairing between
$H^{-1}(\Omega)$ and $H^{1}_0(\Omega)$.

 We use the Cauchy-Schwarz and Poincar\'e inequalities to prove that the bilinear
form $ a(\cdot,\cdot)$ is continuous on the space
$[H^1_0(\Omega)]^2\times[H^1_0(\Omega)]^2$, elliptic on $[H^1_0(\Omega)]^2$,
and that the bilinear form $b(\cdot,\cdot)$ is continuous on the space
$[H^1_0(\Omega)]^2\times L^2(\Omega)$. The bilinear form $b(\cdot,\cdot)$
satisfies the following Inf-Sup condition \cite{BRR}:
There exists a positive constant $\gamma$ such that
$$
\forall t\in L^2_0(\Omega),\quad
\sup_{\mathbf{v}\in[H^1_0(\Omega)]^2}\frac{b(\mathbf{v},t)}{\| \mathbf{v}
\|_{[H^1(\Omega)]^2}}\geq \gamma\| t \|_{L^2(\Omega)}.
$$
Then we deduce that for all $\mathbf{f}$ in the space $[H^{-1}(\Omega)]^2$,
the problem \eqref{e2.2} has a unique solution $(\mathbf{u},p)$ in
$ [H^1_0(\Omega)]^2\times L^2_0(\Omega)$, such that
$$
\| \mathbf{u} \|_{[H^1(\Omega)]^2} + \gamma\| p \|_{ L^2(\Omega)} \leq
C \| \mathbf{f} \|_{[H^{-1}(\Omega)]^2},
$$
where $C$ is a positive constant \cite{GR,T}.

The incompressibility condition $\operatorname{div} \mathbf{u}=0$ on the
connected domain $\Omega$ allows us to deduce that there exists a stream
function $\psi$ in the space $H^2_0(\Omega)$ such that \cite{GR,rad1}
$$
\mathbf{u}=\operatorname{curl}(\psi).
$$
This allows us to deduce that the problem \eqref{e2.2} is equivalent to the
following problem:
For $\mathbf{f} \in [L^2(\Omega)]^2$, find $\psi \in H^2_0(\Omega)$ such that
\begin{gather*}
-\nu\Delta^2 \psi= \operatorname{curl}(\mathbf{f}) \quad\text{in }\Omega\\
 \psi =0 \quad\text{on } \Gamma \\
  \frac{\partial\psi}{\partial n} = 0 \quad\text{on } \Gamma\,,
\end{gather*}
where $n$ is the outward normal unit vector to $\Omega$ on $\Gamma$.

Because of its fundamental importance in the study of singularities and the 
regularity of the solution $\psi$, we consider the characteristic equation
 of the bi-Laplacian operator \cite{K,rad2}:
\begin{equation}
\sin ^2\alpha_j z=z^2\sin ^2\alpha_j. \label{e2.3}
\end{equation}

Next, we consider a partition of the domain $\Omega$ in rectangles 
$\Omega_i$, $1 \leq i \leq I$, satisfying
$$
\bar\Omega= \cup_{i=1}^{I}\bar \Omega_i,\quad 
\Omega_i\cap \Omega_j=\emptyset \text{ for } i\neq j.
$$
We suppose that the intersection of each $\bar \Omega_i$, $1 \leq i \leq I$,
with the boundary $\Gamma$ is either empty or a corner or one of several entire 
edges of $\Omega_i$. The edge of the $\Omega_i$ are parallel to the coordinate axes.
We choose a non-convex domain, and we assume that the non-convex angle 
$\alpha$ is equal to $3\pi\over 2$ or to $2\pi$ (case of the crack). 
The treatment of the singular function is processed locally, then we suppose 
that the non-convex corner is unique.

\subsection*{Assumption 1}
Let $\mathbf{c}$ be the corresponding corner of $\alpha$. We define $\sum$ 
the open domain in $\Omega$ such that $\bar\sum$ is the union of the 
$\bar\Omega_k$ which contains $\mathbf{c}$. We consider that the origin of the 
coordinate axes is at the point $\mathbf{c}$ and we introduce a system of polar 
coordinates $(r,\theta)$. For technical reasons, we lead to assume the 
following conforming property: if the intersection of $\bar \Omega_i$ and 
$\bar \Omega_j$, $i \neq j$, it contains either $\mathbf{c}$, or both an 
edge of $\Omega_i$ and $\Omega_j$.
\smallskip

If the angle $\alpha$ is equal to $3\pi/ 2$,  equation \eqref{e2.3} 
has two roots in the band $0<\operatorname{Re}(z)<1$. When we approximate 
those roots by the Newton method we find $z_1\simeq 0,544484$ and 
$z_2\simeq 0,908529$.

Let $V$ be a neighborhood of the corner $\mathbf{c}$, we introduce the
functions
$$
\eta_1(r, \theta)=r^{1+z_1}\beta_1(\theta), \quad
 \eta_2(r, \theta)=r^{1+z_2}\beta_2(\theta),
$$
where
\begin{gather*}
\begin{aligned}
\beta_1(\theta)&=2.093\big(\cos(0.459\theta)
-\cos(1.544\theta)\big) + 1.093\big(2.193\sin(0.459\theta)\\
&\quad -0.647\sin(1.544\theta)\big),
\end{aligned} \\
\begin{aligned}
\beta_2(\theta)&=4.302\big(\cos(0.092\theta(-\cos(1.908\theta)\big)
 - 1.815\big(10.869\sin(0.092\theta) \\
&\quad -0.524\sin(1.908\theta)\big).
\end{aligned}
\end{gather*}
If $\mathbf{f}$ belongs to $[L^{2}(\Omega)]^2$, the velocity $\mathbf{u}$ 
is decomposed as:
$$
\mathbf{u}=\mathbf{u}_r+\mathbf{u}_s
$$
such that $\mathbf{u}_r$ is in the space $H^2(\Omega)\cap H^1_0(\Omega)$ and
$$
\mathbf{u}_s= \lambda_1 \kappa_1 + \lambda_2 \kappa_2
$$
where $ \kappa_i(r,\theta)=\operatorname{curl}(\eta_i(r,\theta)) 
\in [H^{(1+z_i)-\epsilon}(V)]^2$, $\lambda_i\in \mathbb{R}$ are the two
singularity coefficients, $i \in \{1,2\}$, for all $\epsilon >0$.

In the case of the crack, for $\mathbf{f}$ in $[L^2(\Omega)]^2$, the 
velocity is written in the form
$$
\mathbf{u}=\mathbf{u}_r+\mathbf{u}_s,
$$
$\mathbf{u}_r$ is in $[H^2(\Omega)\cap H_0^1(\Omega)]^2$ and there exists two real
constants $\lambda_1$ and $\lambda_2$ such that
$$
\mathbf{u}_s=\lambda_1 \kappa_1+\lambda_2 \kappa_2
$$
with
\begin{gather*}
\kappa_1(r,\theta)=r^{1/2}(3\sin\theta\sin {\theta\over 2}, 3(1-\cos\theta)\sin
{\theta\over 2}), \\
\kappa_2(r,\theta)=r^{1/2}(2 \sin {\theta\over 2}+\sin\theta\cos{\theta\over 2}, 
(1-\cos\theta)\cos {\theta\over 2}),
\end{gather*}
where $\kappa_1$ and $\kappa_2$ belong to the space $[H^{{3\over 2}-\epsilon}(V)]^2$,
for all $\epsilon$ positive \cite{C2,K}.

\section{Discrete problem and the error estimate}

We start by recalling the basic results of the mortar spectral element method. 
Since the discretization is based on the Galerkin method, we have to define 
the discrete problem and give the quadrature formula which is used for the 
numerical integration.

We denote by $\delta=(N_1,N_2,\dots,N_I)$ the discretization parameter where 
$N_i\geq 2$, $1\leq i\leq I$. Let $\mathbb{P}_N(\Omega_i)$ be the space of 
polynomials on $\Omega_i$. The restriction of the discrete functions to 
$\Omega_i$ will belong to $\mathbb{P}_{N_i}(\Omega_i)$.

Let us recall the Gauss-Lobatto quadrature formula: there exists a unique set 
of nodes such that $\xi_0=-1$, $\xi_N=1$, $\xi_j^N \in ]-1,1[$,
$1\leq j\leq (N-1)$, and a set of $(N+1)$ positive weights $\rho_j^N$,
$0\leq j\leq N$, such that for any polynomial $\varphi_{2N-1}$ with degree 
less than or equal to $2N-1$,
\begin{equation}
\int_{-1}^1\varphi_{2N-1}(\varsigma)d\varsigma
=\sum_{j=0}^N\varphi_{2N-1}(\xi_j^N)\rho_j^N.\label{e3.1}
\end{equation}
Let $F^i$ be the affine bijection from $]-1,1[^2$ into $\Omega_i$.
We consider the local discrete scalar product:
For $u, v$ continuous functions on $\bar\Omega_i$,
\begin{equation}
(u,v)_{N_i}=\frac{|\Omega_i|}{4}\sum_{j=0}^{N_i}\sum_{l=0}^{N_i}
(uoF^i)(\xi_j^{N_i},\xi_l^{N_i})(voF^i)(\xi_j^{N_i},\xi_l^{N_i})
\rho_j^{N_i}\rho_l^{N_i},\label{e3.2}
\end{equation}
where $|\Omega_i|$ is the measure of $\Omega_i$.

We consider $\Gamma^{i,j}$, $1\leq i\leq I$, $1\leq j\leq 4$, the edges of 
$\Omega_i$ and we denote by $\gamma^{il}=\bar\Omega_i\cap\bar\Omega_l$,
$i\neq l$, $\gamma^{il}$ is not necessarily an entire edge $\Gamma^{i,j}$ 
since the decomposition is in general not conforming. We suppose that the 
boundary $\partial \Omega$ is composed of entire edges of the $\Omega_i$.

We define the skeleton of the decomposition:
$$
S=\cup_{k=1}^K\gamma_k\quad \text{and}\quad
\gamma_k\cap\gamma_{k'}=\emptyset\text{ for } k\neq k',
$$
where $\gamma_k$ is called mortar. 
$\gamma_k$ is assumed to be an entire edge of one rectangle $\Omega_i$, 
denoted by $\Omega_{i(k)}$. For any nonnegative integer n and for any segment 
$\gamma$, $\mathbb{P}_n(\gamma)$ is the space of polynomials with degree less 
or equal to $n$ on $\gamma$. We define $W_\delta$ the mortar functions space by:
$$
W_\delta=\{\psi\in L^2(S): \forall k, 1 \leq k \leq K,\psi/_{\gamma_k}\in
 \mathbb{P}_{N_{i(k)}}(\gamma_k)\}.
$$
The discrete space $X_\delta$ represents the space of functions 
$w_\delta$ in $L^2(\Omega)$ such that we have the following properties: 
\cite[Chap 3, \S 1]{BMP}


$\bullet$ the restriction of $w_\delta$ to $\Omega_i$, $1\leq i\leq I$, 
belongs to $\mathbb{P}_{N_i}(\Omega_i)$;

$\bullet$ $w_\delta$ vanishes on $\partial\Omega$;

$\bullet$ the mortar function $\phi$ defined on $S$ by
$$
\phi/_{\gamma_k}=w_\delta/\Omega_{i(k)}/\gamma_k,\; 1\leq k\leq K,
$$
satisfies, for $1\leq i\leq I$ and for any edge $\Gamma$ of $\Omega_i$ 
contained in $S$:
for all $\chi\in \mathbb{P}_{N_i-2}(\Gamma)$,
\begin{equation}
\int_{\Gamma}(w_\delta/_{\Omega_i}-\phi)(\tau)\chi(\tau) d\tau=0.\label{e3.3}
\end{equation}

Let $Y_\delta=X_\delta\times X_\delta$ be the discrete space of the discrete 
velocity. For the discrete pressure we consider the space:
$$
M_\delta=\{p_\delta\in L^2(\Omega): p_{\delta/\Omega_i}\in 
\mathbb{P}_{N_i-2}(\Omega_i)\text{ and }
\int_\Omega p_\delta(\mathbf{x})d\mathbf{x}=0\}.
$$
The space $M_\delta$ corresponds to the case where the pressure has no 
spurious modes (see \cite{BBCM,MPR}). We define the following norm on $Y_\delta$:
$$
\| \mathbf{w}_\delta\|=\Big(\sum_{i=1}^I\| \mathbf{w}_{\delta/\Omega_i}
\|^2_{[H^1(\Omega_i)]^2}\Big)^{1/2}.\\
$$
We are expanding the discrete space of the velocity $Y_{\delta }$. 
Let ${\kappa}_1=(\kappa_1^1,\kappa_1^2)$ be the first singular function of 
the velocity \cite{C3}. We denote by $Y_{\delta }^*$ the space
$$
Y_{\delta }^*= Y_\delta + \mathbb R{\kappa}_1.
$$
If $\mathbf{u}_\delta^*$ is in $Y_{\delta}^*$, there exists 
$\mathbf{u}_\delta=(u_{\delta 1},u_{\delta 2})$ in $Y_{\delta }$ and 
$\lambda_1$ in $\mathbb{R}$ such that
$$
\mathbf{u}_{\delta}^*=\mathbf{u}_{\delta }+\lambda_1 {\kappa}_1.
$$
Since $\kappa_1$ is in the space $[H^1(\Omega)]^2$, we consider the norm 
on the space $Y_{\delta}^*$ for all 
$\mathbf{u}_\delta^*= \mathbf{u}_{\delta }+\lambda_1 {\kappa}_1$ in 
$Y_{\delta}^*$:
$$
\| \mathbf{u}_\delta\|_*=\big(\| \mathbf{u}_\delta\|^2+
 | \lambda_1|^2 \| {\kappa}_1 \|^2\big)^{1/2}.
$$
Therefore, the discrete problem is defined as follows:
for a continuous data function $\mathbf{f}=(f_1,f_2)$ on $\bar \Omega$, find
$\mathbf{u}_\delta^*=(u_{\delta 1}^*,u_{\delta 2}^*)=(u_{\delta 1}+\lambda_1 \kappa_1^1,u_{\delta 2}+\lambda_1 \kappa_1^2)$ in $Y_{\delta
}^*$ and $p_\delta$ in $M_\delta$ such that for all $\mathbf{v}_\delta^*=(v_{\delta
1}^*,v_{\delta 2}^*)=(v_{\delta 1}+\mu_1 \kappa_1^1,v_{\delta 2}+\mu_1 \kappa_1^2)$ in $Y_{\delta}^*$ and for all $q_\delta$ in
$M_\delta$,
\begin{equation}
\begin{gathered}
 a^*_\delta(\mathbf{u}^*_\delta,\mathbf{v}^*_\delta)+
b_\delta^*(\mathbf{v}^*_\delta,p_\delta)
=(\mathbf{f},\mathbf{v}^*_\delta)_\delta\\
b^*_\delta (\mathbf{u}^*_\delta,q_\delta)=0,
\end{gathered} \label{e3.4}
\end{equation}
such that
$$
a^*_\delta(\mathbf{u}^*_\delta,\mathbf{v}^*_\delta)
=a^*_{1\delta}(u^*_{\delta1},v^*_{\delta1})
+ a^*_{2\delta}(u^*_{\delta2},v^*_{\delta2}),
$$
where the bilinear form $a^*_{k\delta}(\cdot,\cdot)$,
$k\in \{1,2\}$ is defined by
\begin{align*}
a_{k\delta}^*(u_{\delta k}^*,v_{\delta k}^*)
&= \sum_{i=1}^I\Bigl(
(\nabla u_{\delta k},\nabla v_{\delta k})_{N_i}
+\lambda_1\int_{\Omega_i} \nabla \kappa^k_1\nabla v_{\delta k} dx
+\mu_1 \int_{\Omega_i}\nabla u_{\delta k} \nabla \kappa^k_1 dx \\
&\quad +\lambda_1 \mu_1\int_{\Omega_i}(\nabla \kappa^k_1)^2dx\Bigr).
\end{align*}
Since we know $\kappa_1^k$, $k \in \{1,2\}$, we can compute exactly the value
of the integral $ {\int_{\Omega_i}}(\nabla \kappa_1^k)^2dx$.
To approximate the integral
$ {\int_{\Omega_i}} \nabla \kappa^k_1\nabla v_{\delta k} d\mathbf{x}$, and
$$
(\mathbf{f},\mathbf{v}^*_\delta)_\delta
=(f_1,v^*_{\delta 1})_\delta+(f_2,v^*_{\delta 2})_\delta
$$
such that
$$
(f_k,v^*_{\delta k})_\delta=(f_k,v_{\delta k})_\delta
+ \mu_1  {\int_{\Omega}}f_k(x,y)\kappa_1^k(x,y) d\mathbf{x}
$$
we use the algorithm presented in \cite{ABM}.

As $\operatorname{div}\kappa_1=0$, for 
$\mathbf{u}^*_{\delta}=\mathbf{u}_{\delta}+\lambda_1 {\kappa}_1$ in 
$Y^*_{\delta}$ and $q_{\delta}$ in $M_\delta$, we have
$$
b^*_\delta(\mathbf{u}^*_\delta,q_\delta)
=- \big (\sum_{i=1}^I(\operatorname{div}\mathbf{u}_\delta,q_\delta)_{N_i} 
+ \lambda_1 \int_{\Omega_i}\operatorname{div}{\kappa}_1 q_\delta dx\big )
=b_\delta(\mathbf{u}_\delta,q_\delta).
$$
Let
$$
V^*_\delta=\{\mathbf{v}^*_\delta \in Y^*_\delta,\;
 b^*_\delta(\mathbf{v}^*_\delta,q_\delta)=0,\; \forall q_\delta \in M_\delta\}
$$
be the kernel of the bilinear form $b^*_\delta(\cdot,\cdot)$.

The bilinear form $ a^*_\delta(\cdot,\cdot)$ is continuous on the space $Y_\delta^*$ by respecting the norm $\|.\|_*$. There exists a positive constant $\nu$ independent of $\delta$ such that for all $\mathbf{u}_\delta^*, \mathbf{v}_\delta^*$ in $Y_\delta^*$
$$
| a_\delta^*(\mathbf{u}_\delta^*, \mathbf{v}_\delta^*)|\leq \nu\|
\mathbf{u}_\delta^*\|_*\| \mathbf{v}_\delta^*\|_*.
$$
In order to prove that the problem \eqref{e3.4} is well posed, we need to define the following norm on the space $Y^*_\delta$:
$$
\| \mathbf{v}_\delta^*\|_{**} =\Bigl(\sum_{i=1}^I
{\| \mathbf{v}_\delta^*}_{/\Omega_i}\|_{[H^1(\Omega_i)]^2}^2
 \Bigr)^{1/2}\quad\mbox{for all $\mathbf{v}_\delta \in Y^*_\delta$}.
$$

The bilinear form $ a^*_\delta(\cdot,\cdot)$ is elliptic relative to the norm 
$\|\cdot\|_{**}$. There exists a positive constant $\beta$ independent 
of $\delta$ such that for all $\mathbf{u}_\delta^*$ in $Y_\delta^*$,
$$
 a_\delta^*(\mathbf{u}_\delta^*, \mathbf{u}_\delta^*)
\geq \beta\| \mathbf{u}_\delta^*\|_{**}^2.
$$

We have the following Inf-Sup condition on the bilinear form 
$b_\delta^*(\cdot,\cdot)$, which establishes the compatibility between the 
spaces $Y_\delta ^*$ and $M_\delta$: there exists a constant 
$\gamma_\delta=\inf(N_i^{-{1\over 2}}),\ 1\leq i\leq I$ such that \cite{BBCM}
$$
\forall t_\delta \in M_\delta , \quad
 \sup_{\mathbf{u}_\delta^*\in
Y_\delta^*}^{}{b_\delta^*(\mathbf{u}_\delta^*,t_\delta)\over\|
\mathbf{u}_\delta^*\|_*}
\geq \gamma_\delta\| t_\delta \|_{L^2(\Omega)}.
$$
Then we can conclude that  problem \eqref{e3.4} has a unique solution 
$(\mathbf{u}_\delta^*,p_\delta)$ in the space $Y_\delta^*\times M_\delta$ such that
$$
\| \mathbf{u}_\delta^* \|_{**} + \gamma_\delta\| p \|_{ L^2(\Omega)} \leq
C \| \mathbf{f} \|_{[L^{2}(\Omega)]^2},
$$
where C is a positive constant.

Since the equivalence constant for the two norms $\| \cdot \|_{*}$ and 
$\| \cdot \|_{**}$ depends on the discrete parameter $\delta$, we prefer 
to keep the norm $\| \cdot\|_{*}$ and proving an Inf-Sup condition on 
the bilinear form $ a^*_\delta(\cdot,\cdot)$. 
There exists a constant $\mu > 0$ independent of $\delta$ such that
 \cite{C1} for all $\mathbf{u}_\delta^* \in V_\delta^* $,
\begin{equation}
\sup_{\mathbf{v}_\delta^*\in
V_\delta^*}^{}{ a^*_\delta (\mathbf{u}_\delta^*,\mathbf{v}_\delta^*)\over\|
\mathbf{v}_\delta^*\|_*}
\geq \mu \| \mathbf{u}_\delta^* \|_{L^2(\Omega)}.\label{e3.5}
\end{equation}
From the Inf-Sup condition \eqref{e3.5} on the bilinear form
$ a^*_\delta(\cdot,\cdot)$ and the Strang's lemma, we conclude the
following error estimate
\begin{equation}
\begin{aligned}
\| \mathbf{u}-\mathbf{u}_\delta^*\|_*
&\leq  C\Big(  {\inf_{v_\delta^*\in V_\delta^*}}
\Big\{ \| \mathbf{u}-\mathbf{v}_\delta^*\|_*+ {\sup_{\mathbf{w}_\delta^*\in\
Y_\delta^*}}{( a- a_\delta^*)(\mathbf{v}_\delta^*,\mathbf{w}_\delta^*)\over\|
\mathbf{w}_\delta^*\|_*}\Big\}\\
&\quad + {\sup_{\mathbf{w}_\delta^*\in\ Y_\delta^*}}
 {{ {\int_{\Omega_i}} \mathbf{f} \mathbf{w}_\delta^* \,dx\,dy
 - (f ,\mathbf{w}_\delta^* )_\delta} \over\| \mathbf{w}_\delta^*\|_*}\\
&\quad + {\sup_{\mathbf{w}_\delta\in
V^*_\delta}}{{{ {\sum_{i=1}^I} {\sum_{i=1}^I}\int_{\gamma^{il}}
({-\nu{\partial \mathbf{u}\over\partial n}}+pn)[\mathbf{w}^*_\delta]}\over{\|
\mathbf{w}^*_\delta\|_*}}}\Big),
\end{aligned}\label{e3.6}
\end{equation}
where $\mathbf{u}$ is the solution of  \eqref{e2.1},
$\mathbf{u}_\delta^*$ is the solution of  \eqref{e3.4},
$[\mathbf{w}^*_\delta]$ is the jump of $\mathbf{w}^*_\delta$ through the
 edge $\gamma^{il}$ and $C$ is a positive constant independent of $\delta$.

If we decompose the solution as:
$$
\mathbf{u}=\mathbf{u}_r +\lambda_1{\kappa_1}+\lambda_2 {\kappa}_2,
$$
and we estimate each term of \eqref{e3.6}, the error between the 
continuous velocity $\mathbf{u}$ of  problem \eqref{e2.1} and the discrete 
velocity $\mathbf{u}^*_\delta $, of problem \eqref{e3.4} is
$$
\| \mathbf{u}-\mathbf{u}^*_\delta\|_* 
\leq C\sup\Big\{\sum_{i=1}^I
N_i^{-\sigma_i},\sum_{i=1}^I N_i^{-\rho_i}\Big\}\| \mathbf{f}\|_{[H^{s-2}
(\Omega)]^2}
$$
where $\mathbf{f}$ belongs to the space $[H^{s-2}(\Omega)]^2$, $s>0$, 
such that $\mathbf{f}/_{\Omega_i}$ in the space 
$[H^{\rho_i}(\Omega_i)]^2$, $\rho_i > 1$ and for all $\epsilon > 0$, 
$\sigma_i$ is given by
$$
\sigma_i=\begin{cases}
s-1 &\text{if } \bar \Omega_i \text{ contains no corner }\bar \Omega_i\\
\inf\{s-1, 2 z_2({\pi\over 2}) -\epsilon\}& \text{if $\bar \Omega_i$  
contains corners different of }\mathbf{c} \\
\inf\{s-1, 2 z_2(\alpha)-\epsilon\} &\text{if $\bar \Omega_i$ contains }\mathbf{c}
\end{cases}
$$
$z_2(\alpha)$ is the second real solution of the equation \eqref{e2.3} 
in the band $0<\operatorname{Re}(z)<1$ for the angle $\alpha$.


For a regular data $\mathbf{f}$ and $N= {\inf_{1\leq i\leq I}} N_i$, 
then for all $\epsilon$ positive we remark that
\begin{itemize}
\item  If $\alpha= \frac{3\pi}{2}$, then the order of convergence is
$N^{\epsilon-1,816}$;

\item  If $\alpha= {2\pi}$, then the order of convergence is $N^{\epsilon-1}$.
\end{itemize}
We conclude that we double the order of convergence when we compare with 
the case without the Strang and Fix algorithm \cite{BM1,C3}.

\section{Implementation and numerical results}

We present in this section the implementation of the Strang and Fix 
algorithm for the spectral element method of the Stokes problem. 
We consider just the first singular function on the velocity for 
$\alpha = \frac{3\pi}{2}$ and the case of the crack $\alpha = {2\pi}$. 
To simplify the implementation we neglect the singularity of the pressure.
 We start by choosing the basis for the discrete spaces of the 
velocity and the pressure. Then, we write the matrix system corresponding 
to the discrete problem \eqref{e3.4}. Afterward, we describe the iterative
 method for solving the matrix system. Finally, we present some numerical 
results showing that the Strang and Fix algorithm is an efficient method. 
We notice that the computation are performed with MATLAB code.

\subsection{Choice of basis}

To define the matrix system of the discrete problem, we have to choose 
a basis of the discrete spaces $Y_{\delta}^*=X_{\delta}^*\times X_{\delta}^*$
 and $M_{\delta}$. These basis are naturally defined through a local basis 
(on each sub-domain) and determine the structure of the matrix system. 
The choice of the resolution method is inherent in this structure.
Let $\varphi_j^N$ the Lagrange interpolating polynomials defined on $[-1,1]$ 
of degree less or equal to $N$ such that
$$
\varphi_j^N\in \mathbb P_N([-1,1]),\quad 
\varphi_j^N(\xi_l)=\delta_{lj},\quad 0\leq l,j\leq N
$$
where $\delta_{lj}$ is the Kronecker symbol. The polynomial $\varphi_j^N$ 
is defined as
$$
\varphi_j^N(\xi)={-1\over N (N+1)} {(1-\xi^2) L_N^{'}(\xi)\over(\xi-\xi_j)}
\quad \forall\xi\in[-1,1],
$$
where $L_N$ is the Legendre polynomials defined on $[-1,1]$ of degree less than
or equal to $N$. We denote $\xi_l^i=F^i(\xi_l)$ and the polynomial
 $\varphi_l^{N_i}$ satisfies $\varphi_l^{N_i} \circ F^i=\varphi_l^N$, 
then for all $v_\delta$ in the space $X_\delta$,
$$
v_\delta(x,y)_{/\Omega_i}=\sum_{l=0}^{N_i}\sum_{j=0}^{N_i}v^{l,j}_{N_i}
\varphi_l^{N_i}(x)
\varphi_j^{N_i}(y)
$$
where $v^{l,j}_{N_i}=v_\delta(\xi_l^i,\xi_j^i)_{/\Omega_i}$.

For each $v_\delta^*$ in $X_\delta^*$, there exists 
$(v_\delta, \lambda_1)$ in $X_\delta\times \mathbb R$ such that 
$v_\delta^*=v_\delta+\lambda_1 \kappa^k_1$, $k\in \{1,2\}$ then
\begin{equation}
v_\delta^*(x,y)_{/\Omega_i}=\sum_{l=0}^{N_i}\sum_{j=0}^{N_i}v^{l,j}_{N_i}
\varphi_l^{N_i}(x)
\varphi_j^{N_i}(y)+ \lambda_1 \kappa^k_{1/\Omega_i}.\label{e4.1}
\end{equation}

Now, we have to choose a basis for the space of the discrete pressure 
$M_\delta$. Let the polynomial $\bar \varphi^{N}_j$ in the space 
$\mathbb P_{N-2}(]-1,1[)$ defined by
$$
\bar \varphi^{N}_j(x) = {\varphi^N_j(x)\over (1-x^2)}.
$$
We recall the following equalities that will be used later. 
For $0\leq l,j\leq N$ we have
\begin{gather*}
\varphi^{N'}_l (\xi_j)={L_N(\xi_j)\over L_N(\xi_l)(\xi_j-\xi_l)},\quad 
\varphi^{N'}_l (\xi_l)=0\\
\bar \varphi^{N}_l (\xi_0)= {1\over 2} \varphi_l' (\xi_0), \quad
\bar \varphi^{N}_l (\xi_N)= {1\over 2} \varphi_l '(\xi_N).
\end{gather*}
The family of polynomials $\bar \varphi^N_l \bar \varphi^{N}_j$ for 
$1\leq l,j\leq N-1$
is a basis of the space $\mathbb P_{N-2}(\Omega)$. 
Let $\bar \varphi^{N_i}_l$ the polynomials such that 
$\bar \varphi^{N_i}_l\circ F^i=\bar \varphi^{N}_l$, then for all 
$q_\delta$ belongs to $M_\delta$
\begin{equation}
q_\delta(x,y)\slash_{\Omega_i}=\sum^{N_i-1}_{l=1} \ \sum^{N_i-1}_{j=1} \
q^{l\,j}_{N_i}\bar \varphi^{N_i}_l(x)\bar \varphi^{N_i}_j(y).\label{e4.2}
\end{equation}
where
$q^{l\,j}_{N_i}=q_\delta(\xi^{i}_l,\xi^{i}_j)(1-{\xi^{i}_l}^2)(1-{\xi^{i}_j}^2)$
(see \cite{AC}).

\subsection{Matrix system and algorithm of resolution}

In this section we formulate the discret problem \eqref{e3.4} as a matrix 
system using the basis of the velocity and the pressure spaces.

For $u^*_\delta$ in the space $Y^*_\delta $ we have 
$u^*_{\delta}=u_\delta+\lambda_1\kappa_1 $, where $\lambda_1$ 
is the first singular coefficient related to the first singular function 
$\kappa_1=(\kappa_1^1,\kappa_1^2)$. Although in the reality we have only 
one singular coefficient $\lambda_1$ in $\mathbb{R}$, we will suppose that we have 
two unknowns $\lambda_1^1$ and $\lambda_1^2$  to maintain the 
symmetry of the problem. The numerical results of $\lambda_1^1$ and $\lambda_1^2$ 
are approximatively equal to $\lambda_1$.

Having \eqref{e4.1} and \eqref{e4.2} we obtain the  matrix system
\begin{equation}
\begin{gathered}
A  U+ B^T P=F\\
B  U=0
\end{gathered}\label{e4.3}
\end{equation}
where $(U,P)$ are, respectively, the vector of admissible unknowns of
the velocity and the vector of admissible unknowns of the pressure.
The components of $U$ are the values of the solution on the nodes of the
sub-domains and the nodes on the corresponding boundaries. The components
of $P$ are the values of the pressure on the internal nodes of the sub-domains.

\subsubsection*{Construction of the matrix $A$}
We note that the matrix $A$ is written in the form
$$
A=\begin{pmatrix}
\tilde{A}_1 \quad 0\\
0 \quad \tilde{A}_2
\end{pmatrix}
$$
where the matrices $\tilde{A}_k $, for $k\in\{1,2\}$ are symmetric matrices 
such that:
\begin{itemize}
\item  the diagonal entries are 
$[\tilde A_k]_{i,i}=\Big(\nabla(\varphi^{N_i}_l \varphi^{N_i}_j) ;
\nabla(\varphi^{N_i}_p \varphi^{N_i}_q )\Big)_{N_i}$,  
$1\leq l,j\leq N_i-1$ and $1\leq p,q\leq N_i-1$, which represents the 
Laplace operator with the Neumann boundary condition on each sub-domain $\Omega_i$,
  $1\leq i\leq I$;

\item the entries on the last row and last column are
 $[\tilde A_k]_{i,{I+1}}=[\tilde A_k]_{{I+1},i}
={\int_{\Omega_i}\nabla\kappa^k_1\nabla(\varphi^{N_i}_p \varphi^{N_i}_q )}\,dx\,dy$ 
where $1\leq i\leq N_{\Sigma}$, and  $N_{\Sigma}$ is the number of
 rectangles contained  in ${\sum}$;
         
\item the bottom right entry is 
$[\tilde A_k]_{{I+1},{I+1}}={\int_\Omega(\nabla \kappa^k_1)^2 \,dx\,dy}$;
        
\item the other entries are zero.
    \end{itemize}


\subsubsection*{Construction of the matrix $B$}
The matrix $B$ is written as
$B=[B_1,\;B_2 ]$ where $ B_k= [\mathcal{B}, 0 ]$, $k\in\{1,2\}$,
$\mathcal{B}$ is composed by the matrices 
\[
\mathcal{B}^i=-\Big(\operatorname{div}(\varphi^{N_i}_l \varphi^{N_i}_j,
\varphi^{N_i}_l \varphi^{N_i}_j),(\bar \varphi^{N_i}_p 
\bar \varphi^{N_i}_q)\Big)_{N_i},
\] 
for $1\leq l,j\leq (N_i-1)$, $1 \leq p,q\leq N_i-1$ and 
$1\leq i\leq I$, which represents the gradient operator on each 
sub-domain $\Omega_i$.

In the following $i$ is omitted for simplicity,
\begin{align*}
&-\Big(\operatorname{div}(\varphi^{N}_l \varphi^{N}_j,\varphi^{N}_l \varphi^{N}_j),
(\bar \varphi^{N}_p \bar \varphi^{N}_q)\Big)_{N} \\
&= - \sum_{m=0}^N  \sum_{n=0}^N \Big({\varphi^N_{l}}'(\xi_m)\varphi^N_j(\xi_n)
  +\varphi^{N}_l(\xi_m){\varphi^N_j}'(\xi_n)\Big){\bar \varphi}^N_p(\xi_m){\bar
\varphi}^N_q (\xi_n) \rho_m \rho_n,
\end{align*}
thus for $p\neq l$ we have
\begin{align*}
\Big( {\varphi^N_l}'\,\varphi^N_j,\bar \varphi^N_p\bar \varphi^N_q \Big)_N 
&=  -{1\over 2} {\rho_j\over(1-\xi^2_j)L_N(\xi_l)}\Big( {\rho_0\over
L_N(\xi_p)(\xi_0-\xi_l)(\xi_0-\xi_p)} \\
&\quad +{2L_N(\xi_p)\rho_p\over (\xi_p-\xi_l)(1-{\xi_p}^2)}
 -{\rho_N\over(\xi_N-\xi_p)(\xi_N-\xi_l)L_N(\xi_p)}\Big)\delta_{jq}
\end{align*}
and for $p=l$,
\begin{align*}
&\Big(  {\varphi^N_l}' \varphi^N_j,\bar \varphi^N_p \bar
 \varphi^N_q \Big)_N \\
&=-{1\over 2} 
{\rho_j\over(1-\xi^2_j)L_N(\xi_l)L_N(\xi_p)}
\Big({\rho_0\over( \xi_0 - \xi_l)(\xi_0-\xi_p)} - {\rho_N\over(\xi_N-\xi_p)
(\xi_N-\xi_l)}\Big)\delta_{jq}.
\end{align*}
We note that if $l=j=0$ and $l=j=N$ are made in the same way, with some
 modifications. We denote by
$F=\big(F_1\; F_2\; 0\big) ^{\rm T}$
the second member where $F_k, k\in \{1,2\}$ are given by
$$
F_k= \begin{pmatrix} (\varphi^{N_1}_p\varphi^{N_1}_q,f_k)_{N_1} \\
 \dots \\
 (\varphi^{N_I}_p\varphi^{N_I}_q,f_k)_{N_I} \\
 {\int_{\Omega}f_k \kappa_1^k dx}\\
\end{pmatrix}.
$$

\subsubsection*{Solution of the matrix system}

System \eqref{e4.3} has false degrees of freedom, which are the values of the
 velocity on the boundary nodes of the sub-domains. These values are found 
by the action of the matching matrix $\widetilde Q$ on the mortar functions. 
They are linked by the integral matching condition \eqref{e3.3}. 
The calculation of this matrix is local for each pair edge-mortar.

We consider $\phi$ a mortar function such that:
$$
\phi_{/\gamma_k}=\sum_{j=0}^{N_{i(k)}}\phi_j^k \varphi_j^{N_{i(k)}}(s),\quad 
1\leq k \leq K, \quad s\in [-1,1].
$$
Now, we determine a basis of the test function space $\mathbb{P}_{N_i-2}(\Gamma)$ that identifies with
$\mathbb{P}_{N_i-2}(]-1,1[)$
$$
\chi_{/\Gamma}=\sum_{p=1}^{N_i-1}\chi_p \eta_p^{N_i-2}(s)
$$
where
$$
\eta_p^{N_i-2}(s)={(-1)^{N_k-p}L_{N_i}(s)\over (\xi_p-s)},\quad
p \in \{1,\dots ,N_i-1\},
$$
(more details are given in \cite{SS}). Hence, the integral matching 
condition \eqref{e3.3} is written as $w_\delta = \widetilde Q\phi$. 
Then, we present the matching global matrix $Q$.
$$
\underbrace{\begin{pmatrix}
w_{lj}^i/\text{internal} & \\
w_{lj}^i/\text{boundary} & \\
\lambda_1
 \end{pmatrix}}_{w^*_\delta}
=\underbrace{\begin{pmatrix}
I &0 &0\\
0 &\widetilde Q &0 \\
0 &0 &1
\end{pmatrix}}_Q
\underbrace{\begin{pmatrix}
w_{lj}^i/\text{internal} & \\
 \phi^k_j \\
\lambda_1\\
\end{pmatrix}}_{\tilde{w^*_\delta}},
$$
for $1\leq i\leq I$,  and $1\leq k\leq K $.

The matrix $Q$ is used to uncouple the system \eqref{e4.3}  to be solved locally.
 The transpose of the matrix $Q$ purges the vector of unknowns from the false 
degrees of freedom. So the resulting system is
\begin{equation}
\begin{gathered}
\mathbf{Q}^T  A  \mathbf{Q} \tilde{U}+\mathbf{Q}^T  B^T P
=\mathbf{Q}^T  F \\
B  \mathbf{Q}  \tilde{U} =0,
\end{gathered}\label{e4.4}
\end{equation}
where $\mathbf{Q}=\begin{pmatrix}Q &0\\ 0 &Q \end{pmatrix}$,
$\tilde{U}$ is composed with the values of the velocity on the internal
nodes and the mortar functions on the skeleton. To solve the system \eqref{e4.4}
 we apply the Uzawa algorithm which is appropriate for the Stokes problem
(see \cite{GR,T}). We consider
$$
\tilde{A}=\mathbf{Q}^T A \mathbf{Q},  \tilde{B}=B\mathbf{Q} \text{ and }
\tilde{F} =\mathbf{Q}^T F.
$$
We uncouple the two unknowns $\tilde{U}$ and $P$. The matrix $\tilde{A}$
is invertible since the bilinear form $a_\delta^*(\cdot,\cdot)$
is elliptic. The first equation of the system \eqref{e4.4} is written as
\begin{equation}
\tilde{U}=\tilde{A}^{-1}\big(\tilde{F}-\tilde{B}P\big).\label{e4.5}
\end{equation}
By inserting \eqref{e4.5} in the second equation of the system \eqref{e4.4}
 we obtain
\begin{equation}
\big(\tilde{B}\tilde{A}^{-1}\tilde{B}^T\big)P=\tilde{A}^{-1}\tilde F.\label{e4.6}
\end{equation}
To find the discrete pressure $P$, the linear system \eqref{e4.6} is solved
by the gradient conjugate method since the matrix
 $\tilde{B} \tilde{A}^{-1}  \tilde{B}^T$ is symmetric and positive define.
 We apply the same method on the system \eqref{e4.5}  to find the discrete
velocity $\tilde{U}$.

\noindent\textbf{Uzawa algorithm}
\begin{itemize}
\item Let $P_0$ arbitrary.

\item Initialization step:
$\tilde{A}\tilde{U}_0=\tilde{F}-\tilde{B}^T P_0$

\item Iterations: $m\geq0$
From $\tilde{U}_m$ and $P_m$:
\begin{gather*}
G_m=  -\tilde{B} \tilde{U}_m\\
\tilde{A} V_m =\tilde{B}^T G_m\\
\rho_m ={\| G_m\| ^2 \over (\tilde{B}^T G_m, V_m)}\\
P_{m+1} =P_m-\rho_m G_m\\
\tilde{U}_{m+1} =\tilde{U}_m+\rho_m V_m.
\end{gather*}
\end{itemize}

\subsection{Numerical results}

In this section we present some numerical tests which are in accordance 
with the theoretical results given in \cite{C3}. We provide the behavior 
of the error between the continuous solution and the discrete one. 
The convergence tests are done with analytical and singular solutions. 
The tests are established on two non convex domains which correspond 
respectively to the case where $\alpha=2\pi$ and $\alpha={{3\pi}\over 2}$ 
see Figure \ref{fig1}.

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.45\textwidth]{fig1a} % cracknnew.pdf
 \includegraphics[width=0.45\textwidth]{fig1b} % Lnnnew.pdf
\end{center}
\caption{Spectral mesh of the domain when $\alpha=2\pi$ (left) and 
$\alpha=\frac {3\pi}{2}$ (right).}\label{fig1}
\end{figure}

If we denote by $N$ the degree of the polynomial in 
$\bar\Omega_i$, $1 \leq i\leq I$ which contains the singular point $\mathbf{c}$.
 Then, the degree of the polynomial in the other rectangular sub-domains 
is chosen less than $N$ and is fixed for each numerical tests.

Next, we present some numerical cases  to show the efficiency of the 
method which consists in the use of the domain decomposition method 
associated to Strang and Fix algorithm.

We consider the viscosity $\nu=1$, and the pressure equal to $p(x,y)=xy$ 
for all the error curves test. 
The given solution for the velocity is the first singular function in the 
two situations:

(i) for $\alpha=3\pi/2$ the first singularity is constructed using the formula
\begin{equation}
\kappa_1(r,\theta)=\operatorname{curl}(\eta_1(r,\theta))
=\big({\cos(\theta)\over r}{\partial\eta_1\over\partial\theta}
+ \sin(\theta){\partial\eta_1\over\partial r},{\sin(\theta)\over r}
{\partial\eta_1\over\partial\theta}
- \cos(\theta){\partial\eta_1\over\partial r}\big) \label{e4.7}
\end{equation}
where
\begin{align*}
\eta_1(r,\theta)&=r^{1.544}\Big( 2.093\big(\cos(0.459\theta)
-\cos(1.544\theta)\big) \\
&\quad + 1.093\big(2.193\sin(0.459\theta)-0.647
\sin(1.544\theta)\big)\Big).
\end{align*}

(ii) for $\alpha=2\pi$ the first singularity is
\begin{equation}
\kappa_1(r,\theta)=r^{1/2}(3\sin\theta\sin {\theta\over 2}, 3(1-\cos\theta)\sin
{\theta\over 2}).\label{e4.8}
\end{equation}
Figure \ref{fig2} (resp. \ref{fig3}) presents the convergence curves of the 
error on the velocity for the solution issued from \eqref{e4.7} (resp. \eqref{e4.8}).
 We present the semi logarithmic scale (for $N$ varying from 5 to 60) and the
 Logarithmic scale for proving the slope. We remark that the obtained results 
for the error and slopes using Strang and Fix algorithm are better than 
those without.
\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig2a} % singmarche.jpg
\includegraphics[width=0.45\textwidth]{fig2b} % singmarchepente.jpg
\end{center}
\caption{Error curves for the solution defined from \eqref{e4.7}}\label{fig2}
\end{figure}

\begin{figure}[htbp]
\begin{center}
 \includegraphics[width=0.45\textwidth]{fig3a} % singcrack.jpg
 \includegraphics[width=0.45\textwidth]{fig3b} % singcrackpente.jpg
\end{center}
\caption{Error curves for the solution defined from \eqref{e4.8}.}\label{fig3}
\end{figure}

In the second test we choose an analytic velocity $u=\operatorname{curl}(\phi)$, 
$\phi$  is the stream function
\begin{equation}
\phi(x,y)=\sin(\pi x)^2\sin(\pi y)^2. \label{e4.9}
\end{equation}
Figure \ref{fig4} corresponds from left to right to the semi logarithmic
error curves with respect to $N$ for $\alpha=3\pi/2$ and $\alpha=2\pi$
in the case of the regular solution issued from \eqref{e4.9}. We obtain a good
convergence with or without Strang and Fix algorithm.
This is so because the Strang and Fix algorithm does not improve the
convergence for a regular function.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig4a} % regmarche.jpg
\includegraphics[width=0.45\textwidth]{fig4b} % regcrack.jpg
\end{center}
\caption{Error curves for the solution from \eqref{e4.9}, for 
$\alpha=3\pi/2$ and $\alpha=2\pi$.}\label{fig4}
\end{figure}

The given solution for the velocity is the second singular function in the 
two situations:
\\
(i) for $\alpha=\frac {3\pi}{2}$ the second singularity constructed using the 
formula
\begin{equation}
\kappa_2(r,\theta)=\operatorname{curl}(\eta_2(r,\theta))
=\big({\cos(\theta)\over r}{\partial\eta_2\over\partial\theta}
+ \sin(\theta){\partial\eta_2\over\partial r},{\sin(\theta)\over r}
{\partial\eta_2\over\partial\theta}
- \cos(\theta){\partial\eta_2\over\partial r}\big) \label{e4.10}
\end{equation}
where
\begin{align*}
\eta_2(r,\theta)&=r^{1.908}\Big(4.302\big(\cos(0.092\theta)
-\cos(1.908\theta)\big)  \\
&\quad - 1.815\big(10.869\sin(0.092\theta)-0.524\sin(1.908\theta)\big) \Big).
\end{align*}
(ii) for $\alpha=2\pi$ the second singularity is
\begin{equation}
\kappa_2(r,\theta)=r^{1/2}(2 \sin {\theta\over 2}
+\sin\theta\cos{\theta\over 2}, (1-\cos\theta)\cos {\theta\over 2}).\label{e4.11}
\end{equation}
We present in Figure \ref{fig5} from left to right the convergence curves
of the error corresponding to the second singular function defined
in \eqref{e4.10} for $\alpha=3\pi/2$ (respectively in \eqref{e4.11}
 for $\alpha=2\pi$). We remark that we do not obtain a good convergence with
 or without Strang and Fix algorithm. If we want to improve this convergence
we have to add the second singular function to the discrete space
$X^*_{\delta}$ which is difficult to implement numerically.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=5.5cm]{fig5a} % sing2marche.jpg
\includegraphics[width=5.5cm]{fig5b} % sing2crack.jpg
\end{center}
\caption{Error curves for the solution defined from \eqref{e4.10}, 
respectively from \eqref{e4.11}.} \label{fig5}
\end{figure}

Figure \ref{fig6} presents from left to right and top to bottom the values of 
the two  components of the velocity and of the pressure corresponding to the 
less regular data function $\mathbf{f}=(f_1,f_2)$, with
\begin{equation}
f_1=|x|^{0.5} \quad f_2=|y|^{0.5}\label{e4.12}
\end{equation}
obtained with $N=60$.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig6a}s % ux.jpg
\includegraphics[width=5.5cm]{fig6b}  % uy.jpg}
\includegraphics[width=0.45\textwidth]{fig6c} % pressure.jpg
\end{center}
\caption{The solution $(u_1,u_2,p)$ for the data function \eqref{e4.12} 
and the mesh $\alpha=3\pi/2$.} \label{fig6}
\end{figure}

\subsection*{Conclusion}
This paper dealt with the implementation of the mortar decomposition domain. 
This technique is applied to the spectral element method for approximating 
the standard velocity pressure formulation of the Stokes problem. 
We implemented Strang and Fix algorithm, which consists to enlarge the 
discrete space of the velocity by the first singular function.
 The obtained error curves confirm that Strang and Fix algorithm applied 
with the mortar spectral method is a good tool to improve the convergence 
order of the error. The extension of this method to the non linear Navier-Stokes 
equations and to the three axisymmetric domain is presently under consideration.

\subsection*{Acknowledgments} 
The authors would like to extend their sincere appreciation to the Deanship of 
Scientific Research at King Saud University
for funding this Research group No (RG-1435-026).

\begin{thebibliography}{00}

\bibitem {ABM} M. Amara, C. Bernardi, M. A. Moussaoui;
\emph{Handling corner singularities by mortar elements method}, 
Applicable Analysis, an International Journal, \textbf{46} (1992), 25-44.

\bibitem {AC} M. Azaiez, G. Coppoletta;
 \emph{Calcul de la pression dans le probl\`eme de Stokes par une 
m{\'e}thode spectrale de "quasi-collocation" {\`a} grille unique}, 
Annales Maghr\'{e}bines de l'Ing\'{e}nieur, \textbf{6} (1992), 41-61.

\bibitem {BR} I. Babu\v{s}ka, M. B. Rosenzweig; 
\emph{A finite element scheme for domains with corners}, 
Numer. Math. \textbf{20} (1972), 1-21.

\bibitem {BS} I. Babu\v{s}ka, M. Suri; 
\emph{The optimal convergence rate of the P-version of the finite element method}, 
SIAM J. Numer. Anal. \textbf{24}  (1987), 750-776.

\bibitem {BBCM} F. Ben Belgacem, C. Bernardi, N. Chorfi, Y. Maday;
 \emph{Inf-sup condition for the mortar spectral element discretization of 
the Stokes problem}, Numer. Math. \textbf{85}  (2000), 257-281.

\bibitem {BM1} C. Bernardi, Y. Maday;
\emph{Polynomial approximation of some singular functions}, 
Applicable Analysis, an International Journal, \textbf{42}  (1991), 1-32.

\bibitem {BM2} C. Bernardi, Y. Maday;
\emph{Spectral method}, P.G.Ciarlet, J.L.Lions eds., Handbook of Numerical Analysis,
 Volume V, North-Holland, Amsterdam, (1996).

\bibitem {BMP} C. Bernardi, Y. Maday, A. T. Patera;
\emph{A new nonconforming approch to domain decomposition: 
the mortar element method}, in Nonlinear Partial Differential Equations 
and their Applications, Coll{\`e}ge de France Seminar, H. Br{\'e}zis, 
J.-L. Lions, eds, (1991).

\bibitem {BRR} F. Brezzi, J. Rappaz, P. A. Raviart;
\emph{Finite dimensional approximation of nonlinear problems. 
Branches of nonsingular solution}, Numer. Math. \textbf{36} (1980), 1-25.

\bibitem {C1} N. Chorfi;
\emph{Handling geometric singularities by mortar spectral element method 
case of the Laplace equation}, Journal of Scientific Computing, Vol. \textbf{18}, 
No. 1 (2003), 25-48.

\bibitem {C2} N. Chorfi;
\emph{Geometric singularities of the Stokes problem}, Abstract and Applied Analysis,
Volume 2014, Article ID 491326, (2014).

\bibitem {C3} N. Chorfi;
\emph{Mortar spectral element discretisation of the Stokes problem in 
domain with corners}, Boundary Value Problems, DOI 10, 1186/s13661-015-0387-4, (2015).

\bibitem{FG} M. Fortin, R. Glowinski;
\emph{M\'{e}thode de Lagrangien augment\'{e}. Application \`a la r\'{e}solution
 num\'{e}rique de probl\`{e}mes aux limites}, vol. 9, Dunod, 1982.

\bibitem {GR} V. Girault, P. A. Raviart;
\emph{Finite Element Methods for the Navier-Stokes Equations, Theory and Algoritms}, 
Springer-Verlag, New York (1986).

\bibitem {G} P. Grisvard;
\emph{Elliptic Problems in Nonsmooth Domains}, Pitman, (1985).

\bibitem {K} V. A. Kondratiev;
\emph{Boundary value problems for elliptic equations in domain with conical
 or angular points}, Trans. Moscow Math. Soc. \textbf{16} (1967), 227-313.

\bibitem {MPR} Y. Maday, A. T. Patera, E. M. Ronquist;
\emph{Optimal Lagrange spectral element method for the Stokes semiperiodic problem}, 
SIAM J. Numer. Anal.

\bibitem{rad1} V. R\u{a}dulescu;
\emph{Nonlinear elliptic equations with variable exponent: old and new}, 
Nonlinear Analysis, Theory, Methods and Applications, 121 (2015), 336-369.

\bibitem{rad2} V. R\u{a}dulescu, D. Repov\v{s};
\emph{Partial Differential Equations with Variable Exponents: Variational
 Methods and Qualitative Analysis}, CRC Press, Taylor and Francis Group,
 Boca Raton FL, (2015).

\bibitem {SF} G. Strang, G. J. Fix;
\emph{An Analysis of the Finite Element Method}, Prentice-Hall, 
Englewood Cliffs (1973).

\bibitem {SS} A. H. Stroud, D. Secrest;
\emph{Gaussian Quadrature Formulas}, Prentice Hall, Englewood Cliffs, 
New Jersey, (1966).

\bibitem {T} R. Temam;
\emph{Theory and Numerical Analysis of the Navier-Stokes Equations}, 
North-Holland Publ. Co., Amsterdam, New York, Oxford, New York, (1977).

\end{thebibliography}

\end{document}