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

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2015 (2015), No. 157, pp. 1--13.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2015 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2015/157\hfil Approximation of the singularity coefficients]
{Approximation of the singularity coefficients of an
 elliptic equation by mortar spectral element method}

\author[N. Chorfi, M. Jleli \hfil EJDE-2015/157\hfilneg]
{Nejmeddine Chorfi, Mohamed Jleli}

\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}

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

\thanks{Submitted April 5, 2015. Published June 12, 2015.}
\subjclass[2010]{35J15, 78M22}
\keywords{Mortar spectral method; singularity; crack}

\begin{abstract}
 In a polygonal domain, the solution of a linear elliptic problem is written
 as a sum of a regular part and a linear combination of singular functions
 multiplied by appropriate coefficients. For computing the leading singularity
 coefficient we use the dual method which based on the first singular dual function.
 Our aim in this paper is the approximation of this leading singularity coefficient
 by spectral element method which relies on the mortar decomposition domain technics.
 We prove an optimal error estimate between the continuous and the discrete
 singularity coefficient. We present numerical experiments which are in perfect
 coherence with the analysis.
\end{abstract}

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

\section{Introduction}

If the data are smooth, the solution of an elliptic partial differential 
equation is not regular when the domain is polygonal. For a Dirichlet problem 
of the Laplace operator, we define some singular functions depending only on 
the geometry of the domain. The solution is written as the sum of a regular 
part and singular functions multiplied by appropriate coefficients \cite{K,G}. 
For approximating the leading singularity coefficients we use two algorithms. 
The first one is Strang and Fix algorithm \cite{SF}, which consists to add the 
leading singularity function to the discrete space see \cite{C}. 
For the second algorithm we apply the dual method. The numerical computation 
of the leading singularity coefficient has been performed by finite elements, 
see Amara and Moussaoui \cite{AA1,AA2}. This coefficient is physically 
meaningful in solid mechanics (crack propagation). We use the mortar spectral 
element method: the domain is decomposed in a union of finite number of 
disjoint rectangles, the discrete functions are polynomials of high degree 
on each rectangle and are enforced to satisfy a matching condition on the interfaces. 
This technique is nonconforming because the discrete functions are not continuous.
We refer to Bernardi, Maday and Patera \cite{BMP} for the introduction of the 
mortar spectral element method.

An outline of this article is as follows. In the second section, we give  
the dual singular function and the formula for finding the leading coefficient 
of the singularity. This  coefficient  does not depend on the data function 
or the geometry of the domain but it just depends on the solution. 
In section $3$, we present the discrete problem and the discrete leading 
singularity coefficient. The section $4$ is devoted to the estimation of the 
error and  we prove the optimality. Finally,  a results of a numerical test 
are given in Section $5$.

\section{Dual singular function and the coefficient of the singularity}

In the rest of the paper, we assume that our domain  ${\Omega}$ is a polygon 
of $\mathbb{R}^2$ such that there exists a finite number of open rectangles  
$\Omega_k,1 \leq k \leq K$, for which
\begin{equation}
\bar{\Omega}= \cup_{k=1}^{K}\bar \Omega_k\quad\text{and}\quad
\Omega_k\cap \Omega_l=\emptyset \quad \text{for } k\not = l.\label{2.1}
\end{equation}
We suppose  also that the intersection of each $\bar \Omega_k$ 
(for $1 \leq k \leq K$)  with the boundary $\partial \Omega$ is either 
empty or a corner or one or several entire edges of $\Omega_k$. 
We choose the coordinate axis parallel to the edge of the $\Omega_k$. 

Handling the singularity is a local process. Therefore, it is not restricted 
to suppose that $\Omega$ has a unique non-convex corner $\mathbf{a}$ with an
angle $\omega$ equal either to $3\pi/ 2$ or to $2\pi$ (case of the crack). 
We choose the origin of the coordinate axis at the point $\mathbf{a}$,
we introduce a system of polar coordinates $(r,\theta)$ where $r$ stands for 
the distance from $\mathbf{a}$ and $\theta$ is such that the line $\theta=0$
contains an edge of $\partial\Omega$. For more technical proof that will 
appear later we need to make the following conformity assumption: if the 
intersection of  $\bar \Omega_k$ and  $\bar \Omega_l, k \not= l$ contains 
$\mathbf{a}$, then it contains either  $\mathbf{a}$ or both an edge of $\Omega_k$ and $\Omega_l$.
Let $\Sigma$ be the open domain in $\Omega$ such that $\bar\Sigma$ is the union 
of the $\bar\Omega_k$ which contain $\mathbf{a}$.

The model equation  under consideration  is the following  Dirichlet problem 
for the Laplace operator:
\begin{equation}
\begin{gathered}
-\Delta u=f  \quad \text{in }\Omega\\
 u=0  \quad \text{on }\partial\Omega.
\end{gathered} \label{2.2}
\end{equation}
If the data $f$ belongs to $H^{s-2}(\Omega)$, then  the above problem admits 
a unique solution $u$ belongs to $H^{s}(\Omega)$. This solution is decomposed as:
\begin{equation}
u=u_r+\lambda S_1,\label{2.3}
\end{equation}
where the function $u_r$ belongs to $ H^{s}(\Omega)$ for $s<1+{2\pi\over{\omega}}$ 
such that
$$
\| u_r\|_{H^{s}(\Omega)}+\mid\lambda\mid<C \| f\|_{H^{s-2}(\Omega)}.
$$
Here $\lambda$ is the leading singularity coefficient and $S_1$ is the first 
singular function given by
$$
S_1(r,\theta)=\chi(r) r^{\pi\over\omega}\sin({\pi\theta/\omega}),
$$
with $\chi$ is a smooth cut-off function with supported on $\bar\Sigma$ and 
is equal to $1$ in a neighborhood of $\mathbf{a}$ \cite{MP}.

Since the Laplace operator with the homogeneous Dirichlet boundary conditions 
is self-adjoint operator,  we define the dual singular function associated with 
$ S_1$ by:
$$
S_1^{*}(r,\theta)=\chi(r)r^{-\pi/omega}\sin({\pi\theta/\omega}).
$$
 We remark that this function does not belong to $H^1(\Omega)$, however
 $\Delta S_1^*$ belongs to the space $L^2(\Omega)$. This allows us to 
 define the function $\varphi^*$ in $H^1(\Omega)$ solution of the following problem
\begin{equation}
\int_{\Omega} \nabla \varphi^* \nabla \psi \,dx\,dy
=\int_\Omega\Delta S_1^{*}\psi dx \quad\text{ for all  } \psi \in
 H^1_0(\Omega).\label{2.4}
\end{equation}
More details are given in  \cite{M} and  \cite{MP}.



\begin{proposition}  \label{prop2.1}
Let $u$, respectively $\varphi^*$ be the solution of the problem \eqref{2.2}, 
respectively \eqref{2.4}, then the coefficient of the singularity $\lambda$ 
satisfies
\begin{equation}
\lambda \pi=\int_{\Omega}f S_1^{*}\,dx
+\int_{\Omega} u\Delta S_1^{*}\,dx
=\int_{\Omega}f (S_1^{*}-\varphi^{*})\,dx.\label{2.5}
\end{equation}
\end{proposition}

\begin{proof} 
Let $C_1=\Omega\cap B(\mathbf{a},R)$ where $B(\mathbf{a},R)$ is the ball of center
$\mathbf{a}$ and the radius $R$. Since the cut off function  $\chi$ is equal
to $1$ in a neighborhood of $\mathbf{a}$, we can choose $R$ such that:
$\Delta S_1=\Delta S_1^*=0$ in $C_1$. Then,  using  \eqref{2.2} and \eqref{2.3}
 we have
$$
\int_{\Omega}f S_1^*\,dx=\int_{C_1}-\Delta u_r S_1^*\,dx
+\int_{\Omega\setminus C_1}-\Delta u S_1^*\,dx.
$$
Integrating by parts, this yields
$$
\int_{\Omega}f S_1^*\,dx+\int_{C_1}u \Delta S_1^*\,dx
=\int_0^{\omega}\Bigr(\partial_r(u-u_r)S_1^*-(u-u_r)\partial_r(S_1^*)
\Bigr)(R,\theta) r\,d\theta.
$$
Replacing $u-u_r$ with $\lambda S_1$, we obtain
$$
\int_{\Omega}f S_1^*\,dx+\int_{\Omega}u \Delta S_1^*\,dx
=\lambda \int_0^{\omega}\Bigr(\partial_r(S_1)S_1^*
-S_1\partial_r(S_1^*)\Bigr)(R,\theta)=\pi\lambda.
$$
To get  the second equality,  we replace $\psi$ by $u$ in the problem \eqref{2.4} 
and we  integrate by parts.
\end{proof}

\section{Approximation of the singularity coefficient}

In this section we recall some basic notion concerning the spectral element
method and the mortar matching condition. Since the discretization is essentially 
a Galerkin method relying on the variational formulation \eqref{2.4}, we need 
to define the discrete space and give the quadrature formula which is used to 
compute the integrals of polynomials \cite{SS}.

The discretization parameter is a $K$-tuple of integers $(N_1,\dots,N_K)$ 
larger than or equal to $2$ denoted by $\delta$. For any nonnegative integer $n$
and for $1\leq k\leq K$, we denote by $\mathbb{P}_n(\Omega_k)$ the space of polynomials
on $\Omega_k$ such that their degree with respect to each variable $x$ and $y$ 
is less than or equal to $n$. The restriction of discrete functions to 
$\Omega_k$ will belong to $\mathbb{P}_{N_k}(\Omega_k)$.

Let us recall the Gauss-Lobatto quadrature formula: for any positive integer $n$, 
there exists a unique set of $(n+1)$ nodes $\xi_j$ in $[-1,1]$ and of 
$(n+1)$ positive weights $\rho_j$ ( for $ 0\leq j\leq n$ ) such that the 
following equality holds for any polynomial $\phi$ with degree less than or 
equal to $2n-1$:
$$
\int_{-1}^1\phi(z)dz =\sum_{j=0}^n\phi(\xi_j)\rho_j.
$$
If $T^k$ denotes an affine mapping from $(-1,1)^2$ onto $\Omega_k$, 
we define a bilinear form on continuous functions $u$ and $v$ on $\bar\Omega_k$  
as follow:
$$
(u,v)_{N_k}={\mid\Omega_k\mid\over 4}
\sum_{i=0}^{N_k}\sum_{j=0}^{N_k}(uoT^k)(\xi_i,\xi_j)(voT^k)(\xi_i,\xi_j)\rho_i\rho_j,
$$
where $\mid\Omega_k\mid$ is  the area of $\Omega_k$.

We need some more notations to enforce the matching condition. 
Let $\nu$ be the set of all the corners of the $\Omega_k$ for $1\leq k\leq K$. 
We denote by $\Gamma^{k,j}$ for $1\leq j\leq 4$ the edges of $\Omega_k$ and 
$\gamma^{kl}=\bar\Omega_k\cap\bar\Omega_l$ for $ k\not = l$. 
We make the further assumption that the boundary $\partial \Omega$ consists 
of entire edges of the $\Omega_k$. Next we introduce the skeleton of the 
decomposition:
$$
S=(\cup_{k=1}^K \partial \Omega_k) \setminus \partial\Omega,
$$
and we assume that it is a disjoint union of mortars  
$(\gamma_m)\,, 1\leq m \leq M$ ($M$ is positive integer)
$$
S=\cup_{m=1}^M\gamma_m\text{ and }\gamma_m\cap\gamma_{m'}=
\emptyset\quad \text{for } m\not = m',
$$
where each mortar $\gamma_m$ is an entire edge of one rectangle $\Omega_k$, 
denoted by $\Omega_{k(m)}$.

For any nonnegative integer $n$ and for any segment $\gamma$, we denote by 
$\mathbb{P}_n(\gamma)$
the space of polynomials with degree less or equal to $n$ on $\gamma$. 
The mortar space $W_\delta$ is then defined by
$$
W_\delta=\{\phi\in L^2(S), \text{ such that for all  } m, \; 1 \leq m \leq M,\;
\phi/_{\gamma_m}\in \mathbb{P}_{N_{k(m)}}(\gamma_m)\}.
$$
The space $X_\delta$ is then defined as in the standard mortar method \cite{BMP}. 
It is the space of function $v_\delta$ in $L^2(\Omega)$ such that
\begin{itemize}
\item The restriction of $v_\delta$ to $\Omega_k,$ for  $1\leq k\leq K$,  belongs to $\mathbb{P}_{N_k}(\Omega_k)$.
\item The function $v_\delta$ vanishes on $\partial\Omega$.
\item The mortar function $\varphi$ defined on  $S$ by
$$
\varphi/_{\gamma_m}=v_\delta/\Omega_{k(m)}/\gamma_m,\text{ for } 1\leq m\leq M,
$$
satisfies that  for $1\leq k\leq K$ and for any edge $\Gamma$ of $\Omega_k$ 
contained in $S$,
$$
\forall \psi\in \mathbb{P}_{N_k-2}(\Gamma),\;
\int_{\Gamma}(v_\delta/_{\Omega_k}-\varphi)(\tau)\psi(\tau) d\tau=0.
$$
\end{itemize}
Later we will need to define the space
$$
X^-_\delta=\{v_\delta\in X_\delta\text{ such that }v_{\delta/\Omega_k}
\in \mathbb{P}_{N_k-1}(\Omega_k)\text{ for } 1\leq k\leq K\}.
$$

\begin{remark} \rm
In the general case of nonconforming decomposition, the space $X_\delta$ 
is not embedded in $H^1(\Omega)$ since the functions in $X_\delta$ are not 
necessarily continuous through the interfaces. To ensure the continuity in 
a neighbourhood of $\mathbf{a}$, we assume the following conformity assumption:
\begin{quote}
Conformity assumption: If $\mathbf{a}$ is the extremity of an edge 
$\Gamma^{k(m),j(m)}$, this edge necessarily coincides with the edge 
of an another rectangle $\Omega_l$.
\end{quote}
 Then we suppose that $N_{k(m)}\leq N_l$.
\end{remark}

Now, we define the  scalar product on $\Omega$ as follows:
$$
(\varphi,\psi)_\delta=\sum_{k=1}^K(\varphi,\psi)_{N_k}\quad
\text{ for all }\varphi,\psi\in{ \mathcal C}^0(\bar\Omega_k).
$$
Our goal now is to give a more accurate approximation than the Strang 
and Fix algorithm for calculating the singularity coefficient $\lambda$. 
In fact, let $X^*_{\delta}= X_{\delta}+\mathbb{R} S_1$ be the augmented discrete space.
We consider the following two discrete problems:
\begin{itemize}
\item Find a function $u_\delta^*=u_\delta + \lambda S_1$ in $X_\delta^*$ such that
\begin{equation}
 a_\delta^*(u_\delta^*,v_\delta^*)
=\sum_{k=0}^K\int_{\Omega_k}f/_{\Omega_k} v_\delta^*/_{\Omega_k}\,dx\,dy\quad
\text{for all  } v_\delta^*=v_\delta + \mu S_1\in X_\delta^* .\label{3.1}
\end{equation}

\item Find $\varphi^*_{\delta}$ in $X^*_{\delta}$ such that
\begin{equation}
a^*_{\delta}(\varphi^*_{\delta},\psi^*_{\delta})=(\Delta
S^*_1;\psi^*_{\delta})_\delta \quad\text{ for all } 
\psi^*_{\delta}\in X^*_{\delta},\label{3.2}
\end{equation}
where
\begin{align*}
a_\delta^*(u_\delta^*,v_\delta^*)
&=\sum_{k=1}^K\Big((\nabla u_k,\nabla v_k)_{N_k}
+\lambda  \int_{\Omega_k}   \nabla S_1\nabla v_k \,dx\,dy  \\
&\quad + \mu  \int_{\Omega_k} \nabla u_k \nabla S_1 \,dx\,dy 
+\lambda \mu  \int_{\Omega_k} (\nabla S_1)^2\,dx\,dy\Big).
\end{align*}
\end{itemize}
We refer to \cite{C} for the numerical analysis and the implementation
 by the mortar spectral element method of this problems.

\begin{proposition}   \label{prop3.2}
Let  $u,\varphi^*, u^*_{\delta}$ and $\varphi^*_{\delta}$ be the respective 
solutions of the problems  \eqref{2.2},\eqref{2.4}, \eqref{3.1} and \eqref{3.2}, 
then
\begin{equation}
\pi \lambda^*_{\delta}=\int_{\Omega}f
S_1^{*}\,dx+\int_{\Omega} u^*_{\delta}\Delta S_1^{*}\,dx)= \int_{\Omega}f (
S_1^{*}-\varphi^{*}_{\delta})\,dx\label{3.3}
\end{equation}
and
\begin{equation}
\begin{aligned}
 \pi(\lambda- \lambda^*_{\delta})
&=  \sum^{K}_{k=1}\int_{\Omega_k}\nabla(u-u^*_{\delta})_{|\Omega_k}
\nabla(\varphi^*-\varphi^*_{\delta})/_{\Omega_k}\,dx\\
&\quad +\sum_{1\leq k<l\leq K}\int_{\gamma^{kl}}({\partial u\over\partial
n_k})(\varphi^*_{\delta}/_{\Omega_k}-\varphi^*_{\delta}/_{\Omega_l})-
({\partial\varphi^*\over\partial
n_k})(u^*_{\delta}/_{\Omega_k}-u^*_{\delta}/_{\Omega_l})\,d\tau .
\end{aligned} \label{3.4}
\end{equation}
\end{proposition}


\begin{proof}
To  obtain \eqref{3.3},  we proceed in the same way as in Proposition \ref{prop2.1}.
 Now,  using \eqref{2.5} and \eqref{3.3}, we obtain
$$
\pi(\lambda- \lambda^*_{\delta})=\int_{\Omega}f (
\varphi^*_{\delta}-\varphi^{*})\,dx=\sum^{K}_{k=1}\int_{\Omega_k}
\Delta u_{/\Omega_k}
(\varphi^*-\varphi^*_{\delta})_{/\Omega_k}\,dx.
$$
Integrating by part  yields
\begin{equation}
\pi(\lambda- \lambda^*_{\delta})=\sum^{K}_{k=1}\int_{\Omega_k}
\nabla u\nabla(\varphi^*-\varphi^*_{\delta})\,dx+
\sum^{}_{1\leq k<l\leq K}\int_{\gamma^{kl}}({\partial u\over\partial
n_k})(\varphi^*_{\delta_{/\Omega_k}}-\varphi^*_{\delta_{/\Omega_l}})\,d\tau.
\label{3.5}
\end{equation}
Let $\varphi^*_{\delta}=\varphi_{\delta}+\mu S_1$ and
$u^*_{\delta}=u_{\delta}+\lambda S_1$ be in $X^*_{\delta}$. Since
\begin{align*}
a^*_{\delta}(\varphi^*_{\delta},u^*_{\delta})
&=\sum^{K}_{k=1}(\nabla\varphi_{\delta_{/\Omega_k}},\nabla
u_{\delta_{/\Omega_k}})_{N_k}
+\lambda\int_{\Omega_k}\nabla S_1\nabla u_{\delta}\,dx\\
&\quad +\mu\int_{\Omega_k}\nabla \varphi_{\delta} \nabla S_1\,dx
+\lambda\mu\int_{\Omega_k}(\nabla S_1)^2 \,dx.
\end{align*}
If $u_{\delta}$ is in $X^-_{\delta}$ then
$$
\sum^{K}_{k=1}(\nabla\varphi_{\delta}/_{\Omega_k},\nabla
u_{\delta}/_{\Omega_k})_{N_k}=\sum^{K}_{k=1}\int_{\Omega_k}\nabla
\varphi_{\delta_{/\Omega_k}}
\nabla u_{\delta_{/\Omega_k}}\,dx.
$$
Using \eqref{3.2}, we deduce that
\begin{equation}
a^*_{\delta}(\varphi^*_{\delta},u^*_{\delta})
=\sum^{K}_{k=1}\int_{\Omega_k}\nabla
\varphi^*_{\delta}\nabla u^*_{\delta}\,dx
= \sum^{K}_{k=1}\int_{\Omega_k}\Delta S_1^* u^*_{\delta}\,dx.\label{3.6}
\end{equation}
From \eqref{3.6}, $\Delta\varphi^*=\Delta S^*_1$ in the
distributions sense and since $\varphi^*=0$ on $\partial\Omega$, this yields
$$
a^*_{\delta}(\varphi^*_{\delta},u^*_{\delta})=\sum^{K}_{k=1}\int_{\Omega_k}\Delta
\varphi^*_{/\Omega_k}
 u^*_{\delta_{/\Omega_k}}\,dx,
$$
hence, integration by part yields
\begin{equation}
a^*_{\delta}(\varphi^*_{\delta},u^*_{\delta})=\sum^{K}_{k=1}\int_{\Omega_k}\nabla
\varphi^*_{/\Omega_k}
 \nabla u^*_{\delta_{/\Omega_k}}\,dx-
\sum^{}_{1\leq k<l\leq K}\int_{\gamma^{kl}}({\partial \varphi^*\over\partial
n_k})(u^*_{\delta}/_{\Omega_k}-u^*_{\delta}/_{\Omega_l})\,d\tau.
\label{3.7}
\end{equation}
Thanks to \eqref{3.6} and \eqref{3.7}, we conclude that
$$
\sum^{K}_{k=1}\int_{\Omega_k}\nabla
{(\varphi^*-\varphi^*_\delta)_{/\Omega_k}}
 \nabla u^*_{\delta_{/\Omega_k}}\,dx
=\sum^{}_{1\leq k<l\leq K}\int_{\gamma^{kl}}({\partial \varphi^*\over\partial
n_k})(u^*_{\delta}/_{\Omega_k}-u^*_{\delta}/_{\Omega_l})\,d\tau.
$$
Adding this equality to \eqref{3.6}, we obtain de desired result.
\end{proof}

\section{Error estimation}

In the following theorem, we give our result  concerning  the error on 
the approximation of the singularity coefficient.

\begin{theorem} \label{thm4.1}
Let  $s\geq 0$ and $f\in H^{s-2}(\Omega)$. The error between  $\lambda$ and
$\lambda^*_{\delta}$ satisfies the following estimation: For $\varepsilon > 0$
$$
|\lambda-\lambda^*_{\delta}|
\leq C N^{-1}\Big(\sum^{K}_{k=1}N_k^{-\sigma_k}\Big)\| f\|_{H^{s-2}(\Omega)},
$$
where $N=\inf_{1\leq k\leq K}^{} N_k$ and
$$
\sigma_k =\begin{cases}  
s-1 &\text{if $\bar{\Omega}^k$  contains no corners of $\bar{\Omega}$},\\
\inf\{s-1,8-\varepsilon\} &\text{if $\bar{\Omega}^k$ contains corners
different from $\mathbf{a}$},\\
\inf\{s-1,{4\pi\over\omega}-\varepsilon\} &\text{if $\bar{\Omega}^k$ contains 
$\mathbf{a}$}.
\end{cases}
$$
\end{theorem}

\begin{proof}
To estimate $| \lambda- \lambda^*_{\delta}|$, we have to estimate each term 
of the inequality \eqref{3.4}. For the first term, using Cauchy-Schwarz 
and  Poincar{\'e}-Friedrichs inequalities we deduce that
\begin{equation}
\big| \sum^{K}_{k=1}\int_{\Omega_k}\nabla(u-u^*_{\delta})
\nabla(\varphi^*-\varphi^*_{\delta})\,dx \big|
\leq C\| u-u^*_{\delta}\|_*\|\varphi^*-\varphi^*_{\delta}\|_*
\label{4.1}.
\end{equation}
Since $u$ (respectively $u^*_{\delta}$) is the solution of the continuous 
problem \eqref{2.2} (respectively the discrete problem \eqref{3.1}).  
As the same  for $\varphi$ and $\varphi^*_{\delta}$, are respectively the 
solutions of the problems \eqref{2.4} and \eqref{3.2} with second member equal 
to $\Delta S^*_1$ in $L^2(\Omega)$, we conclude by 
\cite[result (5.16)]{C}.

From the continuity of $S_1$, we deduce that for any function
$\varphi^*_{\delta}$ in $X^*_{\delta}$ the jump
$\varphi^*_{\delta_{/\Omega_k}}-\varphi^*_{\delta_{/\Omega_l}}$ is equal to
$\varphi_{\delta_{/\Omega_k}}-\varphi_{\delta_{/\Omega_l}}$ which vanishes on
$\Sigma$ due to the conformity assumption.

We note also that $u=u_r$ on $\Omega\setminus\bar{\Sigma}$.  Hence
$$
\int_{\gamma^{kl}}({\partial u\over\partial
n_k})(\varphi^*_{\delta_{/\Omega_k}}-\varphi^*_{\delta_{/\Omega_l}})\,d\tau
= \int_{\gamma^{kl}}({\partial u_r\over\partial
n_k})(\varphi^*_{\delta_{/\Omega_k}}-\phi)\,d\tau
-\int_{\gamma^{kl}}({\partial u_r\over\partial
n_k})(\varphi^*_{\delta_{/\Omega_l}}-\phi)\,d\tau,
$$
where  $\phi$ is the mortar function associated to $\varphi_{\delta}$. 
So, the estimation of this quantity can be evaluated as in 
\cite[result (5.24)]{C}. If $\Gamma^k$ is not a mortar then
\begin{equation}
\begin{aligned}
&\big|\sum^{}_{1\leq k<l\leq K}\int_{\gamma^{kl}}({\partial u\over\partial
n_k})(\varphi^*_{\delta_{/\Omega_k}}-\varphi^*_{\delta_{/\Omega_l}})\,d\tau\big| \\
&\leq C \sum^{K}_{ k=1}\sum^{4}_{j=1}\inf_{\psi_{k,j}\in
{\mathbb{P}}_{N_{k^{-2}}}(\Gamma^{kj})}^{}\| {\partial u_r\over\partial
n_k}-\psi_{kj}\|_{(H^{1/2}(\Gamma^{kj}))'}
\|\varphi-\varphi^*\|_*,
\end{aligned}\label{4.2}
\end{equation}
by taking $\psi_{k,j}=\Pi_{N_k-2}({\partial u_r\over\partial n_k})$, 
where $\Pi_{N_k-2}$ is the orthogonal projection operator from $L^2(\Gamma^k)$ in
${\mathbb{P}}_{N_{k^{-2}}}(\Gamma^k)$.
Now, let us consider
$$
u_r=\widetilde u_r+\widetilde \lambda S_2,
$$
where $S_2$ is the second singular function \cite{C}, its proved in the same 
way that
\begin{align*}
&\big|\sum^{}_{1\leq k<l\leq K}\int_{\gamma^{kl}}({\partial u\over\partial
n_k})(\varphi^*_{\delta_{/\Omega_k}}
 -\varphi^*_{\delta_{/\Omega_l}})\,d\tau\big| \\
&\leq C \sum^{K}_{ k=1}\sum^{4}_{j=1}\inf_{\psi_{k,j}\in
{\mathbb{P}}_{N_{k^{-2}}}(\Gamma^{kj})}^{}\| {\partial \widetilde u_r\over\partial
n_k}-\psi_{kj}\|_{(H^{1/2}(\Gamma^{kj}))'} \|\varphi-\varphi^*\|_*.
\end{align*}
Using the fact that  $\widetilde u_r$ belongs to $H^s(\Omega)$ for $s$ between 
$3/2$  and $1+{3\pi\over\omega}$, we take 
$\psi_{kj}= \Pi_{N_k-2}({\partial\widetilde
u_r\over\partial n_k})$ in $\mathbb{P}_{N_k-2}(\Gamma^{kj})$ 
to deduce  the desired estimation.

Now, to get the  estimation of the third term, we note that $\Delta S_1^*$ 
belongs to $L^2(\Omega)$. Since the function $\varphi^*$ is the sum of 
$\widetilde\varphi_r^*$ in $H^2(\Omega)$ and a linear combination of  
$S_1$ and $S_2$ then, for the same reasons as above we have
\begin{equation}
\begin{aligned}
&\big|\sum^{}_{1\leq k<l\leq K}\int_{\gamma^{kl}}({\partial \varphi^*\over\partial
n})(u^*_{\delta_{/\Omega_k}}-u^*_{\delta_{/\Omega_l}})\, d\tau\big|\\
&\leq C (\inf_{v_\delta\in {X}_\delta}\| u_r-v_\delta\|_*)
\sum^{K}_{ k=1}\sum^{4}_{j=1}\inf_{{\mathcal X}_{k,j}\in
{\mathbb{P}}_{N_{k^{-2}}}(\Gamma^{kj})}\| {\partial \widetilde\varphi_r\over\partial
n_k}-{\mathcal X}_{kj}\|_{(H^{1/2}(\Gamma^{kj}))'}.
\end{aligned}\label{4.3}
\end{equation}
Since ${\partial \widetilde\varphi_r\over\partial n_k}$ belongs to
$H^{1/2}(\Gamma^{kj})$, we  can choose  
${\mathcal X}_{k,j}=\Pi_{N_k-2}({\partial\widetilde
\varphi_r\over\partial n_k})$ in $\mathbb{P}_{N_k-2}(\Gamma^{kj})$ 
and consider the fact that
$$
\inf_{v_\delta\in {X}_\delta}^{}\| u_r-v_\delta\|_*\leq
\inf_{v_\delta\in {X}_\delta}^{}\| \widetilde
u_r-v_\delta\|_*+\mid\widetilde\lambda\mid\inf_{v_\delta\in
{X}_\delta}^{}\| S_2-v_\delta\|_*.
$$
 We complete the proof by combining \eqref{4.1}, \eqref{4.2} and \eqref{4.3}.
\end{proof}


\begin{remark} \rm
Note that the convergence order is as $N^{\epsilon-3}$ in the case of the crack. 
However, in the case $\omega={3\pi\over2}$, the convergence order is as 
$N^{\epsilon-{11\over3}}$. This proves the high accuracy of the method.
\end{remark}

\section{Implementation and numerical results}

To write the algebraic system of the problem \eqref{3.1}, it is necessary 
to choose a basis of the discrete space $X_{\delta}^*$. This basis is naturally
 defined by the local basis (on each sub-domain).
Let $h_j^N$ be the Lagrange polynomials
$$
h_j^N\in \mathbb{P}([-1,1]),\quad h_j^N(\xi)=\delta_{ij}\quad\text{for }
 0\leq i,j\leq N.
$$
Here $\xi_i$ are  the ($N+1$) Gauss-Lobatto nodes. Then for all 
$v_\delta$ in $X_\delta$
$$
v_\delta(x,y)_{/\Omega_k}=\sum_{i=0}^{N_k}
\sum_{j=0}^{N_k}v^{i,j}_{N_k}h_i^{N_k}(x) h_j^{N_k}(y),
$$
where $v^{i,j}_{N_k}=v_\delta(\xi_i^k,\xi_j^k)_{/\Omega_k}$,
the  $\xi_i^k$, respectively  $h_i^{N_k}$  deduced from $\xi_i$, respectively
from $h_i^{N}$  by translation and dilatation.

Now,  for $v_\delta^*$ in $X_\delta^*$ there exists a function 
$(v_\delta, \lambda)$ belongs to $X_\delta\times \mathbb{R}$
such that $v_\delta^*=v_\delta+\lambda S_1$ and
$$
v_\delta^*(x,y)_{/\Omega_k}=\sum_{i=0}^{N_k}
\sum_{j=0}^{N_k}v^{i,j}_{N_k}h_i^{N_k}(x)
h_j^{N_k}(y)\quad+\quad \lambda S_{1/\Omega_k}.
$$
Note that the internal degrees of freedom $v^{i,j}_{N_k}$ for 
$1\leq i,j \leq(N_k-1)$ are free. However, the boundary degrees of freedom 
found by the transformation of the mortar functions, they are linked by the 
following integral matching condition
\begin{equation}
 \int_{\Gamma^{k,l}}
(v_{\delta/\Omega^k}-\phi)(\tau)\psi (\tau) d\tau=0\quad \text{for all }
 \psi\in \mathbb{P}(\Gamma^{k,l}).\label{5.1}
\end{equation}

Let $ Q$  be  the matrix reflecting the transmission condition on the interfaces 
of the sub-domains. $Q^T$ its  transpose  purges the vector of unknowns from 
the false degrees of freedom. The calculation of this matrix is local for each 
pair edge-mortar.

Now, let $\phi$ be a mortar function such that
$$
\phi_{/\gamma_m}=\sum_{i=0}^{N_{k(m)}}\phi_i^mh_i^{N_{k(m)}}(s),\quad 
1\leq m \leq M,
$$
here $s$ is a local coordinate  defined on the interval $[-1,1]$.

At this stage, we need to determine a basis of the space of test functions 
$\mathbb{P}_{N_k-2}(\Gamma^{k,l})$ identified with
$\mathbb{P}_{N_k-2}(]-1,1[)$. Then
$$
\psi_{/\Gamma^{k,l}}=\sum_{q=1}^{N_k-1}\beta_q \eta_q^{N_k-2}( \tilde s)
$$
where
$$
\eta_q^{N_k-2}(z)={(-1)^{N_k-q}L_{N_k}(z)\over (\xi_q-z)},\quad z \in ]-1,1[,\;
q \in \{1,\dots,N_k-1\},
$$
with $\tilde s$ is a local variable in $(-1,1)$ and $L_N$ is a Legendre 
polynomial of order $N$ ( more details are given in \cite{SS}). 
Hence, the integral matching condition \eqref{5.1} is given as $v_\delta= Q\phi$.

The matrix $Q$ is not built and the  computation  of the vector $v_\delta$ 
on the spectral nodes of the edge $\Gamma^{k,l}$, $1\leq k\leq K$, $1\leq l\leq 4$ 
depends only on the value of $\phi/\gamma^m$.  Then, we introduce the matching 
global matrix
$$
\underbrace{\begin{pmatrix}
v_{ij}^k/\text{internal}   \\
v_{ij}^k/\text{boundary}   \\
\lambda 
 \end{pmatrix}}_{v^*_\delta}
=\underbrace{\begin{pmatrix}
I    &0  &0\\\
0   &\widetilde Q &0 \\\
0 &0 &1\\
\end{pmatrix}}_ Q
\underbrace{\begin{pmatrix}
v_{ij}^k/\text{internal}   \\
 \phi^m_j  \\\
\lambda
\end{pmatrix}} _{\tilde v^*_\delta},
$$
for $1\leq k\leq K$, and $1\leq  m\leq M $.

After giving the basis of the space $X_\delta^*$ and the matrix $Q$, 
we can write the discrete problem \eqref{3.1} in matrix form
\begin{equation}
\widetilde A  u_\delta^*=F \label{5.2}
\end{equation}
where $\widetilde A$ takes the form
$$
\left(\begin{smallmatrix} 
(\nabla(h_i h_j) ;\nabla(h_p h_q ) )_{N_1}
&0               &.       &.       &0         &{\int_{\Omega_1}    \nabla
S_1\nabla(h_p h_q )}\\
 0       &         &        &       &      &.       \\\
.&      &       &       &      &.\\
.&      &       &      &      &      {\int_{\Omega_{N_\Delta}}    \nabla
S_1\nabla(h_p h_q )}\\\
 .&       &       &       &      &0 \\
.&        &       &       &      &. \\
.&        &       &       &      &. \\
 0        &.      &     . &0      & (\nabla(h_i h_j) ;\nabla(h_p h_q ))_{N_k} &0 \\
 {\int_{\Omega_1}    \nabla
S_1\nabla(h_i h_j )}       &.      & {\int_{\Omega_{N_\Delta}}    \nabla
S_1\nabla(h_i h_j)}      &0      &.       &{\int_\Omega(\nabla S_1)^2}
\end{smallmatrix}\right)
$$
The $k$-th block $(\nabla(h_i h_j) ;\nabla(h_p h_q ))_{N_k}$ for 
$1\leq i,j\leq N_k-1$ and $1\leq p,q\leq N_k-1$, represents the Neumann 
Laplace operator on the sub-domain $\Omega_k$, and $F$ is the second member 
given by
$$
F= \begin{pmatrix} 
(h_ph_q,f)_{N_1} \cr
  \vdots \\
      (h_ph_q,f)_{N_k} \\
    {\int_{\Omega}f S_1 dx}
\end{pmatrix}.
$$
Finally, $u_\delta^*$ is  the vector of admissible unknown. 
It is constitute by the values of the solution at a collocation nodes 
in sub-domains and their respective boundary.

Note that  we do not solve the system \eqref{5.2} because it has false degrees 
of freedom. However, the global system that we solve is the following
\begin{equation}
Q^T\widetilde A Q\widetilde u_\delta^*=Q^T F,\label{5.3}
\end{equation}
where $\widetilde u_\delta^*$ is the vector composed from the unknown on 
the internal nodes and the value of the mortar function $\phi$  on the 
skeleton $S$. The matrix $A=Q^T\widetilde AQ$ is symmetric and positive defined. 
Then, we  will use the gradient conjugate method to solve the problem \eqref{5.3}.

For this, let $u_0$ arbitrary, $r_0=Q^TF-A u_0$, $q_0=r_0$ and
\begin{gather*}
\alpha_n={(r_n,r_n)\over(q_n,Aq_n)}\,, \quad
u_{n+1}=u_n +\alpha_n q_n\,,\\
r_{n+1}=r_n-\alpha_n Aq_n\,,\quad
\beta_n={(r_{n+1},r_{n+1})\over(r_n,r_n)}\,,\\
q_{n+1}=r_{n+1}+\beta_n q_n\,.
\end{gather*}
 Although the resolution by the gradient algorithm is global, the calculation 
is made locally. Indeed, the product matrix-vector which is the most expensive 
induces the following elementary products
\begin{equation}
\widetilde A_{ijmn}^k u_m^k+C_{ij}^k\lambda
=\delta_{jn}\rho_q\rho_nD_{qj}u_{mn}^k +
\delta_{im}\rho_m\rho_q D_{qi}D_{qn} u_{mn}^k+C_{ij}^k\lambda,
\label{5.4}
\end{equation}
where
\begin{gather*}
\widetilde A_{ijmn}^k
=(\nabla(h_ih_j),\nabla(h_mh_n))_{N_k},\\
C_{ij}^k=\int_{\Omega_k}\nabla S_1\nabla (h_ih_j)dx, \quad
D_{ij}={dh_j\over dz}(\xi_i)
\end{gather*}
is the monodimensional derivative operator and $\delta_{ij}$ is the Kronecker 
symbol.


The local matrices  $\widetilde A_{ijmn}^k$ are full, thus the coast is 
as $O(N^4)$ operations and $O(N^4)$ memory space. The tensorisation which 
consists to write the product matrix-vector as:
$$
\widetilde A_{ijmn}^k u_m^k=[\delta_{jn}D_{qj}[\rho_q\rho_n[D_{qm}u_{mn}^k]]]+
[\delta_{in}D_{qi}[\rho_q\rho_m[D_{qn}u_{mn}^k]]]
$$
reduces the cost of operations to $O(N^3)$ and the memory space to $O(N^2)$ 
by sub-domain.

We present now some numerical results on the approximation of the solution 
of the problem \eqref{3.1} and the  singularity coefficient by applying 
the dual method. It consists of varying the parameter of discretization $N$ 
and the problem data.

The test cases are done in a neighborhood of the singular corner $\mathbf{a}$, 
i.e. four sub-domains in the case of the crack and three ones in the case
of $\omega = 3 \pi/ 2$  (see Figure \ref{fig1}).

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig1a} % cracknnew.pdf
\includegraphics[width=0.45\textwidth]{fig1b} % Lnnnew.pdf
\end{center}
\caption{The domain when $\omega=2\pi$ and $\omega=3\pi/2$.}
\label{fig1}
\end{figure}

We present below some numerical results related to the calculation of the discrete 
solution of the problem \eqref{3.1} and the leading singularity
coefficient by the dual method. In the following examples, $\lambda^*_{\delta} $ 
denotes the discrete singularity coefficient.


\noindent{\bf Example 1:}
$u(x,y)=\sin(\pi x^2)\sin(\pi y^2)$ and  $\omega={3\pi/2}$.
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
N & 5 &  10 &  15 & 20 & 30\\
\hline
$\lambda^*_{\delta}$ &$8.0\,10^{-4}$&$8.127\, 10^{-7}$&$-2.731\, 10^{-13}$&$-5.021\,10^{-14}$&$7.481\,10^{-14}$ \\
\hline
\end{tabular}
\end{center}
\bigskip

\noindent{\bf Example 2:} $u(r,\theta)=r^{1/2}\cos({\theta/2})$ 
and $\omega={2\pi}$.
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
N & 5 & 15 & 20 & 30\\
\hline
$\lambda^*_{\delta}$ &0.9995. &0.9999. &1. & 1. \\
\hline
\end{tabular}
\end{center}
\bigskip


\noindent{\bf Example 3:} $u(r,\theta)=r^{2/3}\cos({2\theta/3})$ and 
$\omega={3\pi/2}$.
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
N & 5 & 10 & 15 & 20 & 30\\
\hline
$\lambda^*_{\delta}$ &0.9987. &0.9996. &0.9999. & 1. &1. \\
\hline
\end{tabular}
\end{center}
\bigskip

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig2a} % err_crack.jpg
\includegraphics[width=0.45\textwidth]{fig2b} % err_L.jpg
\end{center}
\caption{The error on the solution and the singularity coefficient.}
\label{fig2}
\end{figure}



Figure \ref{fig2}  presents the error curves on the solution of the 
problem \eqref{3.1} (curves in blue) and the curves of error on the singularity 
coefficient (curves in red) in both $\omega=2\pi$ and $\omega=3\pi/2$. 
The continuous solutions are equal to the first singular function which corresponds 
to a singularity coefficient equal to $1$  (Example $2$  and Example $3$). 
Error curves are made in logarithmic scale which permits the computing of the 
convergence order corresponding to the slope of the curve. We remark that 
the convergence order on the singularity coefficient  is better than the one 
of the solution. This order is equal to $2.8976$ for the crack and to $1.9779$ 
for the L-domain. However in the case of the solution, this order is equal to
 $1.9997$ for the crack and to $0.9994$ for the L-domain.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig3a} % iso_crack.jpg
\includegraphics[width=0.45\textwidth]{fig3b} % isonew_L.jpg
\end{center}
\caption{The discrete solution $\omega=2\pi$ and $\omega=3\pi/2$.}
\label{fig3}
\end{figure}


Let $\Gamma_0=\{(r,\theta) \text{ such that } \theta=0  \text{ and  }\theta =\omega\}$.
Figure $3$ presents the iso-values of the discrete solution in the case when 
$\omega=\frac {3\pi}{2}$  for the following Dirichlet problem
\begin{gather*}
-\Delta u =0\quad \text{in }\Omega\\
 u=x\quad \text{on }{\partial\Omega}/{\bar\Gamma_0}\\
 u=0\quad \text{on }{\Gamma_0}
\end{gather*}
and in the case when $\omega=2\pi$ the discrete solution corresponding 
to the problem:
\begin{gather*}
-\Delta u =1\quad \text{in }\Omega\cr
 u=x^2\quad \text{on }\partial\Omega
\end{gather*}
We consider in the Example $4$ the calculation of the leading singularity 
coefficient in the case of the crack for the following Neumann problem:
\begin{gather*}
-\Delta u =f\quad \text{in }\Omega\\
{{\partial u}\over{\partial n}}=g\quad \text{on }\partial\Omega/{\bar\Gamma_0}\\
{{\partial u}\over{\partial n}}=0\quad \text{on }{\Gamma_0}.
\end{gather*}
Under the following condition assuring the existence of the solution
$$
\int_{\Omega} f\, dx+\int_{\partial\Omega} g\, d\tau=0.
$$


\noindent{\bf Example $4$:} $f=0,\ g=x$ on $\partial\Omega/{\Gamma_0}$, 
$g=0$ on ${\Gamma_0}$ and $\omega={2\pi}$.
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
N & 10 & 15 & 20 & 30 & 40\\
\hline
$\lambda^*_{\delta}$ &0.2680. &0.2701. &0.2752. & 0.2787. &0.2787. \\
\hline
\end{tabular}
\end{center}
\bigskip

For  Example $4$,  if we compare these results with those found in the case 
of a discretization with finite element method (see \cite{AA2}), 
we obtain a better precision. This is due to the accuracy of the spectral method.


\subsection*{Conclusion} 
In this work we dealt with the approximation of the leading singularity 
coefficient by mortar spectral element method. The results obtained using 
the dual method are better than those obtained by Strang and Fix algorithm
 (see \cite{C}). This confirm  the theory since the dual method gives us an 
optimal error estimate and bring to light the efficiency of spectral 
discretization of this type of problem.


\subsection*{Acknowledgments}
The authors would like to extend their sincere appreciation to
the Deanship of Scientific Research at King Saud University for
its funding of this research through the Research Group
Project No RGP-1435-034.

\begin{thebibliography}{99}

\bibitem {AA1} M. Amara, M. A. Moussaoui;
\emph{Approximation de coefficient de singularit{\'e}s}, C.R.A.S. Paris 
{\bf313} S{\'e}rie I, pp 335-338, (1991).

\bibitem {AA2} M. Amara, M. A. Moussaoui;
\emph{Approximation  of solution and singularity coefficients for an elliptic
 equation in a plane polygonal domain}, preprint, E.N.S. lyon, (1989).

\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 {C}  N. Chorfi;
\emph{Handling geometric singularities by mortar spectral element method: 
case of the Laplace equation}, Journal of Scientific Computing, Vol {\bf 18}, 
N1, 25-48. (2003).

\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. {\bf16},  pp
227-313, (1967).

\bibitem {MP}  V. G.  Maz'ja, B. A. Plamenvskij;
\emph{On the coefficients in the asymptotic
of solution of elliptic boundary-value problems near conical point},  Dokl.
Akad. Nauk SSSR {\bf219} (1974) Sov. Math. Dolk. {\bf19}, 1570-1574 (1974).

\bibitem {M}  M. Moussaoui;
\emph{Sur l'approximation des solutions du probl\` eme de
Dirichlet dans un ouvert avec coin}, in singularities and constructive, 
Grisvard-Wendland and Whiteman {\'e}d., Lecture notes in math., n 1121 
Springer Verlag (1985).

\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).

\end{thebibliography}

\end{document}