A Group of Immersed Finite Element Spaces \\ For Elliptic Interface Problems \thanks{This research was partially supported by the National Science Foundation (Grant Number DMS 1016313)
Ruchi Guo and Tao Lin
(ruchi91@vt.edu and tlin@vt.edu)
Department of Mathematics
Virginia Tech
We present a unified framework for developing and analyzing immersed finite element (IFE) spaces for solving typical elliptic interface problems with interface independent meshes. This framework allows us to construct a group of new IFE spaces with either linear, or bilinear, or the rotated-$Q_1$ polynomials. Functions in these IFE spaces are locally piecewise polynomials defined according to the sub-elements formed by the interface itself instead of its line approximation. We show that the unisolvence for these IFE spaces follows from the invertibility of the Sherman-Morrison matrix. A group of estimates and identities are established for the interface geometry and shape functions that are applicable to all of these IFE spaces. These fundamental preparations enable us to develop a unified multipoint Taylor expansion procedure for proving that these IFE spaces have the expected optimal approximation capability according to the involved polynomials.
A_Group_of_IFE_Space.tex
\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage[thmmarks,amsmath,amsthm]{ntheorem}

\usepackage{amsfonts, psfrag, graphicx, amssymb,enumitem,indentfirst,float,geometry,color,enumerate,graphicx,subfigure,algorithm,algpseudocode,pifont,hyperref}

\allowdisplaybreaks[4]
\geometry{left=2cm,right=2cm,top=2.5cm,bottom=2.5cm}

\numberwithin{equation}{section}
\numberwithin{figure}{section}

\usepackage{lipsum}
\newcommand\blfootnote[1]{%
  \begingroup
  \renewcommand\thefootnote{}\footnote{#1}%
  \addtocounter{footnote}{-1}%
  \endgroup
}


\newtheorem{hyp}{Hypothesis}
%\newtheorem{algorithm}{Algorithm}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{definition}{Definition}[section]
\newtheorem{thm}{Theorem}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{rem}{Remark}[section]

\newcommand{\commentout}[1]{{}} % for large block comments

%\newtheorem{thm}{Theorem}
%\newtheorem{rem}{Remark}
%\newtheorem{exmp}{Example}
%\newtheorem{lemma}{Lemma}
%\numberwithin{thm}{section}
%\numberwithin{exmp}{section}
%\numberwithin{lemma}{section}
%\numberwithin{rem}{section}

\newcommand{\abs}[1]{\left|#1\right|}
\newcommand{\norm}[1]{\left\|#1\right\|}

\newcommand{\bfA}{{\bf A}}
\newcommand{\bfa}{{\bf a}}
\newcommand{\bfB}{{\bf B}}
\newcommand{\bfb}{{\bf b}}
\newcommand{\bfC}{{\bf C}}
\newcommand{\bfc}{{\bf c}}
\newcommand{\bfD}{{\bf D}}
\newcommand{\bfd}{{\bf d}}
\newcommand{\bfE}{{\bf E}}
\newcommand{\bfe}{{\bf e}}
\newcommand{\bfF}{{\bf F}}
\newcommand{\bff}{{\bf f}}
\newcommand{\bfG}{{\bf G}}
\newcommand{\bfg}{{\bf g}}
\newcommand{\bfH}{{\bf H}}
\newcommand{\bfh}{{\bf h}}
\newcommand{\bfI}{{\bf I}}
\newcommand{\bfi}{{\bf i}}
\newcommand{\bfj}{{\bf j}}
\newcommand{\bfJ}{{\bf J}}
\newcommand{\bfk}{{\bf k}}
\newcommand{\bfK}{{\bf K}}
\newcommand{\bfl}{{\bf l}}
\newcommand{\bfL}{{\bf L}}
\newcommand{\bfM}{{\bf M}}
\newcommand{\bfm}{{\bf m}}
\newcommand{\bfn}{{\bf n}}
\newcommand{\bfN}{{\bf N}}
\newcommand{\bfP}{{\bf P}}
\newcommand{\bfp}{{\bf p}}
\newcommand{\bfQ}{{\bf Q}}
\newcommand{\bfq}{{\bf q}}
\newcommand{\bfR}{{\bf R}}
\newcommand{\bfr}{{\bf r}}
\newcommand{\bfs}{{\bf s}}
\newcommand{\bfS}{{\bf S}}
\newcommand{\bft}{{\bf t}}
\newcommand{\bfT}{{\bf T}}
\newcommand{\bfu}{{\bf u}}
\newcommand{\bfV}{{\bf V}}
\newcommand{\bfv}{{\bf v}}
\newcommand{\bfw}{{\bf w}}
\newcommand{\bfW}{{\bf W}}
\newcommand{\bfx}{{\bf x}}
\newcommand{\bfX}{{\bf X}}
\newcommand{\x}{{\bf x}}
\newcommand{\bfy}{{\bf y}}
\newcommand{\bfY}{{\bf Y}}
\newcommand{\bfz}{{\bf z}}
\newcommand{\bfZ}{{\bf Z}}



\newcommand{\bfbeta}{{\boldsymbol \beta}}
\newcommand{\bfdelta}{\boldsymbol{\delta}}
\newcommand{\bfepsilon}{\boldsymbol{\epsilon}}
\newcommand{\bfgamma}{\boldsymbol{\gamma}}
\newcommand{\bfGamma}{\boldsymbol{\Gamma}}
\newcommand{\bfLambda}{\boldsymbol{\Lambda}}
\newcommand{\bfphi}{\boldsymbol{\phi}}
\newcommand{\bfPhi}{\boldsymbol{\Phi}}
\newcommand{\bfpsi}{\boldsymbol{\psi}}
\newcommand{\bfPsi}{\boldsymbol{\Psi}}
\newcommand{\bfsigma}{\boldsymbol{\sigma}}
\newcommand{\bftau}{\boldsymbol{\tau}}
\newcommand{\bftheta}{\boldsymbol{\theta}}
\newcommand{\bfTheta}{\boldsymbol{\Theta}}
\newcommand{\bfxi}{\boldsymbol{\xi}}

%\input{tlin_symbol}


\begin{document}
\title{A Group of Immersed Finite Element Spaces \\
For Elliptic Interface Problems \thanks{This research was partially supported by the National Science Foundation (Grant Number DMS 1016313)}
\blfootnote{Keywords: Interface problems, discontinuous coefficients, finite element spaces, Cartesian mesh.}
\blfootnote{2010 Mathematics Subject Classification: 35R05, 65N30, 97N50.}}
\author{Ruchi Guo and Tao Lin\\
(ruchi91@vt.edu and tlin@vt.edu) \\
Department of Mathematics \\
Virginia Tech }
\date{}
\maketitle

\begin{abstract}
We present a unified framework for developing and analyzing immersed finite element (IFE) spaces for solving typical elliptic interface problems
 with interface independent meshes. This framework allows us to construct a group of new IFE spaces with either linear, or bilinear, or the rotated-$Q_1$ polynomials. Functions in these IFE spaces are locally piecewise polynomials defined according to the sub-elements formed by the interface itself instead of its line approximation. We show that the unisolvence for these IFE spaces follows from the invertibility of the Sherman-Morrison matrix. A group of estimates and identities are established for the interface geometry and shape functions that are applicable to all of these IFE spaces. These fundamental preparations enable us to develop a unified multipoint Taylor expansion procedure for proving that these IFE spaces have the expected optimal approximation capability according to the involved polynomials.
\end{abstract}


%\blfootnote{Keywords: Interface problems, discontinuous coefficients, finite element spaces, Cartesian mesh.}
%\blfootnote{2010 Mathematics Subject Classification: 35R05, 65N30, 97N50.}

\section{Introduction}

This article presents a unified framework for developing and analyzing a group of immersed finite element (IFE) spaces that use interface independent meshes (such as highly structured \textrm{Cartesian} meshes) to solve interface problems of the typical second order elliptic partial differential equations:
\begin{align}
\label{eq1_1}
 -\nabla\cdot(\beta\nabla u)=f, \;\;\;\; & \text{in} \; \Omega^-  \cup \Omega^+, \\
 u=g, \;\;\;\; &\text{on} \; \partial\Omega,
\end{align}
where, without loss of generality, the domain $\Omega\subseteq \mathbb{R}^2$ is separated by an interface curve $\Gamma$ into two subdomains $\Omega^+$ and $\Omega^-$, the diffusion coefficient $\beta(X)$ is discontinuous such that
\begin{equation*}
\beta(X)=
\left\{\begin{array}{cc}
\beta^- & \text{if} \; X\in \Omega^- ,\\
\beta^+ & \text{if} \; X\in \Omega^+ ,
\end{array}\right.
\end{equation*}
where $\beta^{\pm}$ are positive constants. In addition, the solution $u$ is assumed to satisfy the jump conditions:
\begin{align}
[u]_{\Gamma} &= 0, \label{eq1_3} \\
\big[\beta \nabla u\cdot \mathbf{n}\big]_{\Gamma} &= 0, \label{eq1_4}
%[\beta \nabla u\cdot \mathbf{n}]_{\Gamma} &=& 0, \label{eq1_4}
\end{align}
%\begin{eqnarray}
%[u]_{\Gamma} &=& 0, \label{eq1_3} \\
%\big[\beta \nabla u\cdot \mathbf{n}\big]_{\Gamma} &=& 0, \label{eq1_4}
%%[\beta \nabla u\cdot \mathbf{n}]_{\Gamma} &=& 0, \label{eq1_4}
%\end{eqnarray}
where $\mathbf{ n}$ is the unit normal vector to the interface $\Gamma$, and for every piece-wise function $v$ defined as
\begin{equation*}
%\label{ }
v=
\left\{\begin{array}{cc}
v^-(X) & \text{if} \; X\in \Omega^- , \\
v^+(X) & \text{if} \; X\in \Omega^+ ,
\end{array}\right.
\end{equation*}
we adopt the notation $[v]|_{\Gamma}=v^+|_{\Gamma}-v^-|_{\Gamma}$.

\par
It is well known that the standard finite element method can be applied to interface problems provided that the mesh is formed according to the interface, see \cite{1970Babuska,1996BrambleKing,1998ChenZou} and references therein. Many efforts have been made to develop alternative finite element
methods based on unfitted meshes for solving interface problems. The advantages of using unfitted meshes are discussed in
\cite{2010Bordas,2014MoumnassiFigueredo,2011MohammedSalim,2013NadalAlbelda}. A variety of finite element methods that can use interface independent meshes to solve interface problems have been reported in the literature, see
\cite{2003BabuskaBanerjeeOsborn,
1998BabuskaZhang,
1987BarrettElliott,
2009BastianEngwer,
2001DolbowMoesBelytschko,2009EfendievHou,2009GuyomarchLeeJeon,2015GuzmanSanchezSarkisMathComp,1996BabuskaMelenk} for a few examples. In particular,
instead of modifying the shape functions on interface elements which is an approach to be discussed in this article, methods in \cite{2015BurmanClaus,2016BurmanGuzmanSanchezSarkisHcontrast,2002HansboHansbo,2016WangXiaoXu} employ standard finite element functions but use the Nitsche's penalty along the interface in the finite element schemes.

\commentout{
{\color{red}In the literature, alternative finite element methods have been developed for using unfitted meshes. For example, in stead of modifying the shape functions on interface elements, methods in \cite{2015BurmanClaus,2016BurmanGuzmanSanchezSarkisHcontrast,2002HansboHansbo,2016WangXiaoXu}
employ the Nitsche's penalty along the interface in the finite element schemes. Other finite element methods that can use interface independent meshes to solve interface problems can be found in
\cite{2003BabuskaBanerjeeOsborn,
1998BabuskaZhang,
1987BarrettElliott,
2009BastianEngwer,
2001DolbowMoesBelytschko,2009EfendievHou,2009GuyomarchLeeJeon,2015GuzmanSanchezSarkisMathComp,1996BabuskaMelenk}. The unfitted mesh is also advantageous when the boundary geometry of the computation is complicated or moving during the computation \cite{2014MoumnassiFigueredo,2011MohammedSalim,2013NadalAlbelda}. In addition, the authors in \cite{2010Bordas} review different approaches aiming at reducing the mesh burden in computation.}
}

This article focuses on the immersed finite element (IFE) methods, whose basic idea was introduced in \cite{1998Li}, for those applications where it is preferable to solve interface problems with a mesh independent of the interface,
for example, the Particle-In-Cell method for plasma particle simulations \cite{2005KafafyWangLin,2007KafafyWang,2002LinWang}, the problems with moving interfaces
\cite{2013HeLinLinZhang,2013LinLinZhang1}, and the electroencephalography forward problem \cite{2010VallaghePapadopoulo}. IFE methods for interface problems of other types of partial differential equations can be found in \cite{2015AdjeridChaabaneLin,2013HeLinLinZhang,2012HouLiWangWang,2005LiYang,2011LinLinSunWang,2013LinLinZhang1,2013LinSheenZhang,2016KihyoMoon_thesis,2003YangLiLi}.

The IFE spaces developed in this article are extended from the IFE spaces constructed with linear polynomials \cite{2007GongLiLi,2010KwakWeeChang,2004LiLinLinRogers,2003LiLinWu}, bilinear polynomials
\cite{2009HeTHESIS,2008HeLinLin,2001LinLinRogersRyan}, and rotate-$Q_1$ polynomials \cite{2013LinSheenZhang,2013ZhangTHESIS} using the standard Lagrange type local degrees of freedom imposed either
at element vertices as usual or at midpoints of element edges in the Crouzeix-Raviart way \cite{1973CrouzeixRaviart}. We note that, the local linear IFE space on each triangular interface element constructed here with the Lagrange local degrees of freedom imposed at vertices is very similar to the one recently introduced in \cite{2015GuzmanSanchezSarkisP1}.
The IFE spaces in this article are new because, locally on each interface element, they are Hsieh-Clough-Tocher type macro finite element functions \cite{2001Braess,1966CloughTocher} defined with sub-elements formed by the interface curve itself in contrast to those IFE spaces in the literature defined with sub-elements formed by a straight line approximating the interface curve.

Our research presented here is motivated by two issues. The first issue concerns the general $O(h^2)$ order accuracy for a line to approximate a curve which is a fundamental ingredient for the optimal approximation capability of those IFE spaces in the literature. We hope the study of IFE spaces based on curve sub-elements can shed light on the development of higher degree IFE spaces for which the $O(h^2)$ order is not sufficient. For examples, those techniques incorporating the exact geometry for constructing basis functions \cite{2012LianBordas,2015NguyenAnitescu,2008SevillaSoniaHuerta} may be considered. The second issue is the attempt to unify the fragmented framework for developing and analyzing the IFE spaces in the literature. For IFE spaces based on different meshes, different polynomials, and different local degrees of freedom, we show that their unisolvence, i.e.,
the existence and uniqueness of IFE shape functions, can be established through a uniform procedure related with the invertibility of the Sherman-Morrison matrix. We have derived a group of identities for the interface geometry and shape functions that are applicable to all of these IFE spaces, and this enables us to derive error estimates for the interpolation in these new IFE spaces in a general unified multipoint Taylor expansion approach in which, IFE functions defined according to the given interface actually simplify the analysis because we only need to apply the same arguments to two sub-elements formed by the interface while the analysis for the IFE spaces in the literature has to use a different set of arguments to handle the sub-elements sandwiched between the interface curve and its approximate line.
Also, inspired by \cite{2015GuzmanSanchezSarkisP1}, we have made an effort to show how the error bounds explicitly depend on the maximum curvature of the interface curve and the ratio between $\beta^-$ and $\beta^+$, which are two important problem dependent characteristics effecting the approximation capability of IFE spaces. We note that the dependence of constants in the error bounds on the ratio between $\beta^-$ and $\beta^+$ in our article is similar to the one discussed in \cite{2010ChuGrahamHou}, and we think this coincidence follows from the fact that we analyze the approximation capability of
finite element spaces with Lagrange type degrees of freedom.

%{\color{red} And in the final error estimates, we specify the dependence of the constants on the curvature of the curve and $\beta^{\pm}$ for
% all the types of finite elements considered}.


This article consists of 5 additional sections. In the next section we describe common notations and some basic assumptions used in this article. In Section 3, we derive estimates and identities associated with the interface and the jump conditions in an element. From these estimates, we can see how their bounds explicitly depend on
curvature of the interface and the ratio between $\beta^-$ and $\beta^+$, and how the mesh size $h$ is subject to the interface curvature.
%{\color{red}And we give a criteria for smallness of the mesh size $h$ based on the curvature.}
In Section 4 we present generalized multipoint Taylor expansions for piecewise $C^2$ functions in an interface element. Estimates for the remainders in these expansions are derived in terms of pertinent Sobolev norms. In Section 5, first, we establish the unisolvence of immersed finite element functions constructed with linear, bilinear, Crouzeix-Raviart and rotated-$Q_1$ polynomials, i.e., we show the standard Lagrange local degrees of freedom imposed at the nodes of an interface element can uniquely determine an IFE function that satisfies the interface jump conditions in a suitable approximate sense. Then, we show that the IFE shape functions have several desirable properties such as the partition of unity and the critical identities in Theorem \ref{thm4_3}.
Finally, with a unified analysis, we show that the IFE spaces have the expected optimal approximation capability. In Section 6, we demonstrate features of these IFE spaces by numerical examples.


%many others. However it is time consuming to generate the body-fitting mesh for arbitrary curve interface, especially for moving interface problems. So it is advantageous to develop numerical methods for non-body-fitting meshes, such as the Cartesian meshes. A typical example is the Particle-In-Cell method for plasma particle simulations \cite{2007KafafyWang,2005KafafyWangLin,2002LinWang,2002LinWang} developed by Lin and Wang, which is very efficient for particle tracking. And in 1977, Peskin introduced \emph{the immersed boundary method} \cite{1977Peskin,2002Peskin,1999AlexandrePeskinBerger} based on finite difference method. There have been many other finite difference schemes based on Cartesian meshes proposed, such as \emph{the immersed interface method} \cite{2003DumettKeener,1994LevequeLi,1994LevequeLi,2006LiIto} and so on.
%
%\par One of the most important numerical methods using finite elements formulation for interface problems is the immersed finite element (IFE) method which was first introduced by Li \cite{1998Li}. It allows interfaces to cut the interior of elements and thus, can be applied on non-body-fitting meshes. So there are two types of elements, interface elements and non-interface elements. On the interface elements, the IFE basis functions are constructed to approximate the jump conditions \eqref{eq1_3} and \eqref{eq1_4}. And the standard basis functions are used on those non-interface elements. An arbitrary $p$-th order accuracy was proved by Adjerid and Lin \cite{2009AdjeridLin} for one dimensional case. For two dimensional cases, various articles have been published over the past decade, such as \cite{2012Chou,2010ChouKwakWee,2007GongLiLi,2009HeTHESIS,2009HeLinLin,2011HeLinLin,2012HeLinLin,2013HeLinLinZhang,2004LiLinLinRogers,2003LiLinWu,2001LinLinRogersRyan}.
%
%\par In this article, we construct the two dimensional nonconforming IFE spaces on rectangular meshes. Rannacher and Turek first introduced the two basic nonconforming finite elements in \cite{1992RanacherTurek} using midpoint-value and integral-value degrees of freedom respectively. And Ciarlet's triplet definition \cite{1978Ciarlet} differentiates the two nonconforming elements as $(T,\Pi_T,\Sigma^P_T)$ and $(T,\Pi_T,\Sigma^I_T)$. Refer to \cite{2001Braess,1994BrennerScott,2005Chen} for more details. Compared with the conforming IFE methods, nonconforming IFE methods have better numerical performance around the interface \cite{2013ZhangTHESIS} and are more reliable for some complicated problems from solid and fluid mechanics. In \cite{2013ZhangTHESIS}, Zhang developed nonconforming IFE spaces using both the $(T,\Pi_T,\Sigma^P_T)$ and $(T,\Pi_T,\Sigma^I_T)$  which are piecewise polynomials partitioned by the line connecting the two intersection points of the curve with boundary of the element. Actually all the previous formulation of two dimensional IFE spaces are based on partitioning by the line. Due to the geometric limits, only the second order accuracy can be achieved. Therefore, the main purpose of this article is to develop new nonconforming IFE basis functions partitioned by the curve and to investigate the approximation capability. In our formulation, instead of using weak flux continuity on the line in \cite{2013ZhangTHESIS}, we use the exact flux continuity on an arbitrary point on the curve to construct the basis functions. The size of the linear system from the formulation of basis functions is only one half of the degrees of freedom. Besides our analysis of interpolation error can be simply applied to other IFE spaces, such as bilinear spaces on rectangular meshes and linear spaces on triangular meshes. We hope the methods developed here could motivate the future study of higher order accuracy.


\section{Preliminaries}
Throughout the article, $\Omega\subset\mathbb{R}^2$ denotes a bounded domain as a union of finitely many rectangles. The interface curve $\Gamma$ separates $\Omega$ into two subdomains $\Omega^+$ and $\Omega^-$ such that $\overline{\Omega} = \overline{\Omega^+} \cup \overline{\Omega^-} \cup \Gamma$. For every measurable subset $\tilde \Omega \subseteq \Omega$, let $W^{k,p}(\tilde \Omega)$ be the standard Sobolev spaces on $\tilde \Omega$ associated with the norm $\|\cdot\|_{k,p,\tilde \Omega}$ and the semi-norm $|v|_{k,p,\tilde \Omega}=\|D^{\alpha}v\|_{0,p,\tilde \Omega}$, for $|\alpha|=k$. The corresponding Hilbert space is $H^k(\tilde \Omega)=W^{k,2}(\tilde \Omega)$. When
$\tilde \Omega^s = \tilde \Omega \cap \Omega^s \not = \emptyset, s = \pm$, we let
%Due to the jump conditions, the solution $u$ has lower global regularity than the solution to the usual elliptic PDE, i.e $u\notin W^{k,p}(\Omega)$, for $k\geqslant 2$. However we note in each %subdomain, $u|_{\Omega^{s}}\in H^1(\Omega^{s})$, $s=+,-$, so the following spaces are important:
\begin{equation*}
%\label{ }
PH^k_{int}(\tilde \Omega)=\{ u: u|_{\tilde{\Omega}^s}\in H^k(\tilde \Omega^s),\; s=\pm; \; [u]_\Gamma=0 \; \text{and} \; [\beta\nabla u\cdot \mathbf{n}]_\Gamma=0\; \text{on}\; \Gamma \cap \tilde \Omega \},
\end{equation*}
\begin{equation*}
%\label{ }
PC^k_{int}(\tilde \Omega)=\{ u: u|_{\tilde \Omega^s}\in C^k(\tilde \Omega^s),\; s=\pm; \; [u]_\Gamma=0 \; \text{and} \; [\beta\nabla u\cdot \mathbf{n}]_{\Gamma}=0\; \text{on}\; \Gamma \cap \tilde \Omega \}.
\end{equation*}
The norms and semi-norms to be used on $PH^k_{int}(\tilde \Omega)$ are
\begin{equation*}
\|\cdot\|^2_{k,\tilde \Omega}=\|\cdot\|^2_{k,\tilde \Omega^+}+\|\cdot\|^2_{k,\tilde \Omega^-}, \;\;\;\;\; |\cdot|^2_{k,\tilde \Omega}=|\cdot|^2_{k,\tilde \Omega^+}+|\cdot|^2_{k,\tilde \Omega^-},
\end{equation*}
\begin{equation*}
\|\cdot\|_{k,\infty,\tilde \Omega}=\max(\|\cdot\|_{k,\infty,\tilde \Omega^+} \;,\; \|\cdot\|_{k,\infty,\tilde \Omega^-}), \;\;\;\;\; |\cdot|_{k,\infty,\tilde \Omega}=\max(|\cdot|_{k,\infty,\tilde \Omega^+} \;,\; |\cdot|_{k,\infty,\tilde \Omega^-}).
\end{equation*}

Let $\mathcal{T}_h$ be a Cartesian triangular or rectangular mesh of the domain $\Omega$ with the maximum length of edge $h$. An element $T\in \mathcal{T}_h$ is called an interface element provided the interior of $T$ intersects with the interface $\Gamma$; otherwise, we name it a non-interface element. We let $\mathcal{T}^i_h$ and $\mathcal{T}^n_h$ be the set of interface elements and non-interface elements, respectively. Similarly, $\mathcal{E}^i_h$ and $\mathcal{E}^n_h$ are sets of interface edges and non-interface edges, respectively. In addition, as in \cite{2009HeLinLin}, we assume that $\mathcal{T}_h$ satisfies the following hypotheses when the mesh size $h$ is small enough:

\begin{itemize}[leftmargin=30pt]
  \item [\textbf{(H1)}] The interface $\Gamma$ cannot intersect an edge of any element at more than two points unless the edge is part of $\Gamma$.
  \item [\textbf{(H2)}] If $\Gamma$ intersects the boundary of an element at two points, these intersection points must be on different edges of this element.
  \item [\textbf{(H3)}] The interface $\Gamma$ is a piecewise $C^2$ function, and the mesh $\mathcal{T}_h$ is formed such that the subset of $\Gamma$ in every interface element $T\in\mathcal{T}^i_h$ is $C^2$.
  \item [\textbf{(H4)}] The interface $\Gamma$ is smooth enough so that $PC^2_{int}(T)$ is dense in $PH^2_{int}(T)$ for every interface element $T\in\mathcal{T}^i_h$.
\end{itemize}

On an element $T \in \mathcal{T}_h$, we consider the local finite element space $(T, \Pi_T, \Sigma_T^P)$ with
\begin{eqnarray}
\Pi_T &=& \begin{cases}
\textrm{Span} \{ 1,x,y\}, & \text{for $P_1$ or Crouzeix-Raviart (C-R) finite element functions}, \\
\textrm{Span} \{ 1,x,y,xy\}, & \text{for $Q_1$ (bilinear) finite element functions}, \\
\textrm{Span} \{ 1,x,y,x^2-y^2 \}, & \text{for rotated-$Q_1$ finite element functions},
\end{cases} \label{eq4_1}\\
\Sigma_T^P &=& \{ \psi^P_T(M_i): i \in \mathcal{I}, \; \forall \psi^P_T \in \Pi_T \}, \label{eq4_2}
\end{eqnarray}
where $\mathcal{I} = \{1, 2, \cdots, DOF(T)\}$, $DOF(T) = 3$ or $4$ depending on whether $T$ is triangular or rectangular, $M_i, i\in \mathcal{I}$ are the local nodes to determine
shape functions on $T$, and the super script $P$ is to emphasize the Lagrange type degrees of freedom imposed at the points $M_i$s. For $P_1$ and $Q_1$ finite elements, $M_i = A_i, i\in \mathcal{I}$ where $A_i$'s are vertices of $T$. For C-R and rotated-$Q_1$ finite elements,
$M_i$ is the midpoint of the $i$-th edge of $T$ for $i\in \mathcal{I}$. It is well known \cite{1994BrennerScott,1978Ciarlet, 1973CrouzeixRaviart,1992RanacherTurek} that
$(T, \Pi_T, \Sigma_T^P)$ has a set of shape functions $\psi_i^P(X), i\in \mathcal{I}$ such that
\begin{equation}
\label{eqn_psi_iT_values_bounds}
\psi^P_{i,T}(M_j)=\delta_{ij},~~\norm{\psi^P_{i,T}}_{\infty, T} \leq C,~~\norm{\nabla \psi^P_{i,T}}_{\infty, T} \leq Ch^{-1}, \;\; i,j\in \mathcal{I},
\end{equation}
where $\delta_{ij}$ is the \textit{Kronecker} delta function.

%Besides for every interface element $T\in\mathcal{T}^i_h$, denote $T^+=\Omega^+\cap T$, $T^-=\Omega^-\cap T$ and thus, we can define
%\begin{equation*}
%%\label{ }
%PH^k_{int}(T)=\{ u: u|_{T^s}\in H^k(T^s),\; s=+\; or\; -; \; [u]=0 \; and \; [\beta\nabla u\cdot \mathbf{ n}_{\Gamma}]=0\; on\; \Gamma\cap T \},
%\end{equation*}
%\begin{equation*}
%%\label{ }
%PC^k_{int}(T)=\{ u: u|_{T^s}\in C^k(T^s),\; s=+\; or\; -; \; [u]=0 \; and \; [\beta\nabla u\cdot \mathbf{ n}_{\Gamma}]=0\; on\; \Gamma\cap T \},
%\end{equation*}
%with the norms and semi-norms
%\begin{equation*}
%%\label{ }
%\| \cdot \|^2_{k,T}=\| \cdot \|^2_{k,T^+}+\| \cdot \|^2_{k,T^-} , \;\;\;\;\; | \cdot |^2_{k,T}=| \cdot |^2_{k,T^+}+| \cdot |^2_{k,T^-},
%\end{equation*}
%\begin{equation*}
%\|\cdot\|_{k,\infty,T}=\max(\|\cdot\|_{k,\infty,T^+} \;,\; \|\cdot\|_{k,\infty,T^-}), \;\;\;\;\; |\cdot|_{k,\infty,T}=\max(|\cdot|_{k,\infty,T^+} \;,\; |\cdot|_{k,\infty,T^-}).
%\end{equation*}
Throughout this article, without loss of generality, we assume that $\beta^+\geqslant\beta^-$ and let $\rho=\beta^-/\beta^+\leqslant1$. In addition, on any $T\in \mathcal{T}^i_h$, we use $D$, $E$ to denote the intersection points of $\Gamma$ and $\partial T$ and let $l$ be the line connecting $DE$.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% section2

\section{Geometric Properties of the Interface}

In this section, we discuss geometric properties on interface elements that are useful for developing and analyzing IFE spaces. Let $T$ be an interface element.
As illustrated in Figure \ref{fig2_2}, for a point $\widetilde{X}$ on $\Gamma\cap T$, we let
$\mathbf{n}(\widetilde{X})=(\tilde{n}_x(\widetilde{X}), \tilde{n}_y(\widetilde{X}))$ be the normal of $\Gamma$ at $\widetilde{X}$, and we let $\widetilde{X}_{\bot}\in l$ be the orthogonal projection of $\widetilde{X} \in \Gamma\cap T$ onto the line $l$. Also, for the line $l$, we let $\bar{\mathbf{n}}=(\bar{n}_x,\bar{n}_y)$ be its unit normal vector and, consequently, $\bar{\mathbf{t}}=(\bar{n}_y, -\bar{n}_x)$ is the vector tangential to $l$. Without loss of generality, we assume the orientation of all the normal vectors are from $T^-$ to $T^+$. In addition, we let $\kappa$ be the maximum curvature of the curve $\Gamma$.


%figure3_1
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.4\textwidth}
\centering
\includegraphics[width=2.5in]{interf_elem_3normal_vectors.pdf}
\caption{ The local system }
\label{fig2_2}
\end{minipage}
\begin{minipage}[t]{0.5\textwidth}
\centering
\includegraphics[width=2.5in]{local_coordinates_interf_elem.pdf}
\caption{The geometry of an interface element }
\label{fig2_1}
\end{minipage}
\end{figure}

For $T \in \mathcal{T}_h^i$, without loss of generality, we can introduce a local coordinate system such that
the point on $\Gamma \cap T$ can expressed as $(\xi,\eta) = (\xi, \omega(\xi))$ for a suitable function $\omega$ such that, in this local system, $D$ is origin and its $x$-axis is aligned with $l$, as shown in Figure \ref{fig2_1}. We start from the following lemma which extends
similar results in \cite{2004LiLinLinRogers}:

\begin{lemma}
\label{lem_geo_curv}
Given any $\epsilon\in(0,\frac{\sqrt{2}}{2})$, assume $h\kappa\leqslant\epsilon$, then for any interface element $T\in\mathcal{T}^i_h$, there hold
\begin{eqnarray}
&&|w(\xi)|\leqslant 2(1-2\epsilon^2)^{-3/2}\kappa h^2, \label{eq2_1} \\
&&|w'(\xi)|\leqslant \sqrt{2}(1-2\epsilon^2)^{-3/2}\kappa h. \label{eq2_2}
\end{eqnarray}
\end{lemma}
\begin{proof}
In the local system, let $\xi_E$ be the coordinate of the point $E$. And by the Mean Value Theorem, there is some $\xi_0\in[0,\xi_E]$ such that $\omega'(\xi_0)=0$. Consider a function $g$ as well as its derivative
$$
g(\xi)=\frac{\omega'(\xi)}{\sqrt{1+(\omega'(\xi))^2}}, ~~~ \textrm{and}, ~~~ g'(\xi)=\frac{\omega''(\xi)}{(1+(\omega'(\xi))^2)^{3/2}}.
$$
Note that $g'(\xi)$ is the curvature of $\Gamma$ at $\xi$; hence, we have $|g'(\xi)| \leq \kappa$. Then, by $g(\xi_0)=0$, we have $|g(\xi)|=\abs{\int^{\xi}_{\xi_0}g'(s)ds}\leqslant\int^{\xi}_{\xi_0}|g'(s)|ds\leqslant\sqrt{2}\kappa h=\sqrt{2}\epsilon$, which implies $|\omega'(\xi)|\leqslant\frac{\sqrt{2}\epsilon}{\sqrt{1-2\epsilon^2}}$. By the definition of $\kappa$, we have $|\omega''(\xi)|\leqslant (1-2\epsilon^2)^{-3/2} \kappa$.

Now using the Taylor expansion for $\omega$ around $D$ leads to
$
\omega(\xi)=\omega'(0)\xi+\frac{1}{2}\omega''(\bar{\xi})\xi^2
$
 for some $\bar{\xi}\in[0,\xi_E]$. Note that $\omega(\xi_E)=0$ shows $\omega'(0)=-\frac{1}{2}\omega''(\bar{\xi}_E))\xi_E$. Thus we have $\omega(\xi)=-\frac{1}{2}\omega''(\bar{\xi}_E))\xi_E\xi+\frac{1}{2}\omega''(\bar{\xi})\xi^2$ and therefore, $|\omega(\xi)|\leqslant 2\|\omega''\|_{\infty}h^2\leqslant2(1-2\epsilon^2)^{-3/2}\kappa h^2$, which yields \eqref{eq2_1}. And using the Taylor expansion again for $\omega'$ around $\xi_0$, we have $\omega'(\xi)=\omega''(\tilde{\xi})(\xi-\xi_0)$ for some $\tilde{\xi}$ between $\xi_0$ and $\xi$. Finally we obtain $|\omega'(\xi)|\leqslant \sqrt{2}\|\omega''\|h\leqslant\sqrt{2}(1-2\epsilon^2)^{-3/2}h$.
\end{proof}
%
%\begin{lemma}
%\label{lem2_1}
%There exists constants $C$ only depending on $f$, such that
%\begin{eqnarray}
%&&|f(\xi)|\leqslant Ch^2, \label{eq2_1} \\
%&&|f'(\xi)|\leqslant Ch. \label{eq2_2}
%\end{eqnarray}
%\end{lemma}

%
%\begin{lemma}
%\label{lem2_2}
%For every $u\in PC^2_{int}(T)$, we have
%\begin{equation}
%\label{eq2_5}
%\nabla u^+(\widetilde{X}) = M^-(\widetilde{X}) \nabla u^-(\widetilde{X}), \;\;\;\;\;
%\nabla u^-(\widetilde{X})  = M^+(\widetilde{X}) \nabla u^+(\widetilde{X}).
%\end{equation}
%\end{lemma}

We note that the argument in the lemma above is similar to Assumption 3.14 in \cite{2010ChuGrahamHou} with a minor difference that
a local polar coordinate system is used on the interface element in \cite{2010ChuGrahamHou}.
%Another difference is that we explicitly specify the dependence of constants on the curvature.
The following lemmas provide estimates about various geometric quantities defined at points on $\Gamma \cap T$.

\begin{lemma}
\label{lem2_3}
Given any $\epsilon\in(0,\frac{\sqrt{2}}{2})$, assume $h\kappa\leqslant\epsilon$, then for any interface element $T\in\mathcal{T}^i_h$ and any point $\widetilde{X}\in\Gamma\cap T$, the following inequality holds:
\begin{equation}
\label{eq2_6}
    \| \widetilde{X}-\widetilde{X}_{\bot}\| \leqslant 2(1-2\epsilon^2)^{-3/2}\kappa h^2,
\end{equation}
and for any $\widetilde{X}_1$, $\widetilde{X}_2 \in \Gamma \cap T$, we have
\begin{subequations} \label{eq2_7}
\begin{align}
    & \| \mathbf{ n}(\widetilde{X}_1)-\mathbf{ n}(\widetilde{X}_2)\|\leqslant \sqrt{2}(1+(1-2\epsilon^2)^{-3/2})\kappa h,  \label{eq2_7a} \\
    & \mathbf{ n}(\widetilde{X}_1)\cdot\mathbf{ n}(\widetilde{X}_2) \geqslant 1- 2(1+(1-2\epsilon^2)^{-3/2})^2\kappa^2h^2. \label{eq2_7b}
\end{align}
\end{subequations}
\end{lemma}
\begin{proof}
Estimate \eqref{eq2_6} directly follows from \eqref{eq2_1}. For \eqref{eq2_7a}, we assume $\widetilde{X}_1 = (\xi_1, w(\xi_2))$ and $\widetilde{X}_2 = (\xi_2, w(\xi_2))$ in the local system, respectively. Then we have
\begin{equation*}
%\label{ }
\mathbf{ n}(\widetilde{X}_1)=\frac{1}{\sqrt{1+(w'(\xi_1))^2}}
\left(\begin{array}{c}
-w'(\xi_1) \\
1
\end{array}\right),
 \;\;\;
\mathbf{ n}(\widetilde{X}_2)=\frac{1}{\sqrt{1+(w'(\xi_2))^2}}
\left(\begin{array}{c}
-w'(\xi_2) \\
1
\end{array}\right).
\end{equation*}
By the calculation in Lemma \ref{lem_geo_curv} and Mean Value Theorem, there is some $\bar{\xi}\in[0,\xi_E]$ such that
\begin{equation*}
\begin{split}
%\label{ }
\abs{\frac{w'(\xi_1)}{\sqrt{1+(w'(\xi_1))^2}} - \frac{w'(\xi_2)}{\sqrt{1+(w'(\xi_2))^2}}}
=\frac{|\omega''(\bar{\xi})|}{(1+(\omega'(\bar{\xi}))^{2})^{3/2}}|\xi_1-\xi_2|\leqslant \sqrt{2}\kappa h,
%\leqslant  \abs{\frac{w'(\xi_1)}{\sqrt{1+(w'(\xi_1))^2}}} + \abs{\frac{w'(\xi_2)}{\sqrt{1+(w'(\xi_2))^2}}}
%\leqslant 2\sqrt{2}\kappa_Th,
\end{split}
\end{equation*}
and $\tilde{\xi}\in[0,\xi_E]$ such that
\begin{equation*}
\begin{split}
%\label{ }
\abs{\frac{1}{\sqrt{1+(w'(\xi_1))^2}} - \frac{1}{\sqrt{1+(w'(\xi_2))^2}}} = \frac{|\omega'(\tilde{\xi})| |\omega''(\tilde{\xi})|}{ (1+(\omega'(\tilde{\xi}))^{2})^{3/2} } |\xi_1-\xi_2| % \leqslant 2(1-2\epsilon^2)^{-3/2}\kappa^2 h^2
\leqslant 2(1-2\epsilon^2)^{-3/2}\epsilon \kappa h.
% \leqslant & \abs{\sqrt{1+(w'(\xi_2))^2} - \sqrt{1+(w'(\xi_1))^2}} \\
%= & \frac{|(w'(\xi_2))^2-(w'(\xi_1))^2|}{ \sqrt{1+(w'(\xi_2))^2} + \sqrt{1+(w'(\xi_1))^2} } \leqslant \kappa_T^2h^2.
\end{split}
\end{equation*}
Then \eqref{eq2_7a} follows by applying these estimates in the local coordinate forms of $\mathbf{ n}(\widetilde{X}_1)$ and $\mathbf{ n}(\widetilde{X}_2)$.
Furthermore, by \eqref{eq2_7a} and
\begin{equation*}
\begin{split}
%\label{ }
% 2(1+(1-2\epsilon^2)^{-3/2})^2h^2\geqslant
\| \mathbf{ n}(\widetilde{X}_1)-\mathbf{ n}(\widetilde{X}_2)\|^2
=\| \mathbf{ n}(\widetilde{X}_1)\|^2+\|\mathbf{ n}(\widetilde{X}_2)\|^2-2\mathbf{ n}(\widetilde{X}_1)\cdot\mathbf{ n}(\widetilde{X}_2)
=2-2\mathbf{ n}(\widetilde{X}_1)\cdot\mathbf{ n}(\widetilde{X}_2),
\end{split}
\end{equation*}
we have \eqref{eq2_7b}.
\end{proof}

\begin{rem}\label{rem2_1}
Note that there exists a point $\widetilde{X}_1 = (\xi_1, w(\xi_1))\in \Gamma\cap T$ such that $w'(\xi_1)=0$ which means $\mathbf{n}(\widetilde{X}_1) = \bar{\mathbf{n}}$. Hence, by Lemma \ref{lem2_3}, we have the following estimates for an arbitrary point $\widetilde{X} \in \Gamma \cap T$:
\begin{subequations}
\begin{align}
%\label{eq2_8}
    & \| \mathbf{ n}(\widetilde{X})-\bar{\mathbf{ n}}\|\leqslant  \sqrt{2}(1+(1-2\epsilon^2)^{-3/2})\kappa h,  \label{eq2_8a}\\
    & \mathbf{ n}(\widetilde{X})\cdot \bar{\mathbf{ n}} \geqslant 1- 2(1+(1-2\epsilon^2)^{-3/2})^2\kappa^2 h^2. \label{eq2_8b}
\end{align}
\end{subequations}
\end{rem}

The two lemmas above have suggested a criteria about how small $h$ should be according to the maximum curvature $\kappa$ of $\Gamma$ so that the related analysis is valid.
Therefore, for all discussions from now on, we further assume that
\begin{itemize}[leftmargin=30pt]
\item
$h$ is sufficiently small such that for some fixed parameter $\epsilon\in (0,\sqrt{2}/2)$ and $\bar{\kappa} \in (0, 1]$ of one's own choice, there holds
\begin{equation}
\label{h_small_assump}
h<\min{ \left\{ \frac{\sqrt{\bar{\kappa}}}{\sqrt{2}(1+(1-2\epsilon^2)^{-3/2})\kappa}, \; \frac{\epsilon}{\kappa} \right\} }.
\end{equation}
\end{itemize}
Obviously $\epsilon$ is the proportion by which we should choose the mesh size $h$ according to the interface curvature $\kappa$. Also, by \eqref{h_small_assump}
and \eqref{eq2_8b}, we have
\begin{eqnarray}
\mathbf{ n}(\widetilde{X})\cdot \bar{\mathbf{ n}} \geq 1 - \bar{\kappa} \label{n_nbar_angle_lower_bnd}
\end{eqnarray}
which shows how much the angle between the normal of $\Gamma\cap T$ and $\bar{\mathbf{ n}}$ can vary in an interface element $T \in \mathcal{T}^i_h$, a larger value of
$\bar{\kappa}\in (0, 1]$ allows $\mathbf{ n}(\widetilde{X})$ to vary more from $\bar{\mathbf{ n}}$ up to, but not equal to, $90$ degree. Therefore, we will call $\bar{\kappa}$ the angle allowance.

%Lemma \ref{lem2_3} shows that \eqref{eq2_6} and \eqref{eq2_7} hold for $h$ satisfying the assumption \eqref{h_small_assump}.
In the rest of this article, all the generic constants $C$ are assumed to possibly depend only on the parameter $\epsilon$ and $\bar{\kappa}$, but they are independent of the interface location, $\beta^{\pm}$, and the curvature $\kappa$.

We now consider some matrices associated with the normal of interface $\Gamma$ and the normal of $l$. First, for any $\widetilde{X}\in\Gamma\cap T$, we use the normal $\mathbf{n}(\widetilde{X})$ to form two matrices:
\begin{equation*}
%\label{ }
N^s(\widetilde{X})=\left(\begin{array}{cc}
\tilde{n}_y(\widetilde{X}) & -\tilde{n}_x(\widetilde{X}) \\
\beta^s \tilde{n}_x(\widetilde{X}) & \beta^s \tilde{n}_y(\widetilde{X})
\end{array}\right), ~~s = \pm.
%\;\;\; and \;\;\;
%N^+(\widetilde{X})=\left(\begin{array}{cc}
%\tilde{n}_y(\widetilde{X}) & -\tilde{n}_x(\widetilde{X}) \\
%\beta^+ n_x(\widetilde{X}) & \beta^+n_y(\widetilde{X})
%\end{array}\right).
\end{equation*}
Since $Det(N^s(\widetilde{X}))=\beta^s, s = \pm$, these matrices are nonsingular; therefore, we can define another two matrices at the
point $\widetilde{X}\in\Gamma\cap T$:
\begin{equation}
\label{eq2_3}
M^-(\widetilde{X})=\left(N^+(\widetilde{X})\right)^{-1}N^-(\widetilde{X})=
\left(\begin{array}{cc}
\tilde{n}^2_y(\widetilde{X})+\rho \tilde{n}^2_x(\widetilde{X}) & (\rho-1)\tilde{n}_x(\widetilde{X})\tilde{n}_y(\widetilde{X}) \\
(\rho-1)\tilde{n}_x(\widetilde{X})\tilde{n}_y(\widetilde{X}) & \tilde{n}^2_x(\widetilde{X})+\rho \tilde{n}^2_y(\widetilde{X})
\end{array}\right),
\end{equation}
\begin{equation}
\label{eq2_4}
M^+(\widetilde{X})=\left(N^-(\widetilde{X})\right)^{-1}N^+(\widetilde{X})=
\left(\begin{array}{cc}
\tilde{n}^2_y(\widetilde{X})+1/\rho \tilde{n}^2_x(\widetilde{X}) & (1/\rho-1)\tilde{n}_x(\widetilde{X})\tilde{n}_y(\widetilde{X}) \\
(1/\rho-1)\tilde{n}_x(\widetilde{X})\tilde{n}_y(\widetilde{X}) & \tilde{n}^2_x(\widetilde{X})+1/\rho \tilde{n}^2_y(\widetilde{X})
\end{array}\right).
\end{equation}
For matrices $M^-(\widetilde{X})$ and $M^+(\widetilde{X})$, we recall from \cite{2004LiLinLinRogers} the following results
\begin{equation}
\label{eq2_5}
\nabla u^+(\widetilde{X}) = M^-(\widetilde{X}) \nabla u^-(\widetilde{X}), \;\;\;\;\;
\nabla u^-(\widetilde{X})  = M^+(\widetilde{X}) \nabla u^+(\widetilde{X}), ~\forall \widetilde{X}\cap T \in \Gamma, ~~\forall u\in PC^2_{int}(T).
\end{equation}

%\begin{lemma}
%\label{lem2_2}
%For every $u\in PC^2_{int}(T)$, we have
%\begin{equation}
%\label{eq2_5}
%\nabla u^+(\widetilde{X}) = M^-(\widetilde{X}) \nabla u^-(\widetilde{X}), \;\;\;\;\;
%\nabla u^-(\widetilde{X})  = M^+(\widetilde{X}) \nabla u^+(\widetilde{X}).
%\end{equation}
%\end{lemma}


%Now take $\bar{\mathbf{ t}}=(\bar{n}_y, -\bar{n}_x)^T$ to be the tangential vector of $l$.
In addition, for $\widetilde{X}\in\Gamma\cap T$, we can use the normal vectors $\mathbf{n}(\widetilde{X})$ and $\bar{\mathbf{n}}$ to form the following matrices:
\begin{equation*}
%\label{ }
\overline{N}^s(\widetilde{X})=\left(\begin{array}{cc}
\bar{n}_y & -\bar{n}_x \\
\beta^s\tilde{n}_x(\widetilde{X}) & \beta^s\tilde{n}_y(\widetilde{X})
\end{array}\right), ~s = \pm.
%\;\;\; and \;\;\;
%\overline{N}^+(\widetilde{X})=\left(\begin{array}{cc}
%\bar{n}_y & -\bar{n}_x \\
%\beta^+\tilde{n}_x(\widetilde{X}) & \beta^+\tilde{n}_y(\widetilde{X})
%\end{array}\right).
\end{equation*}
By Remark \ref{rem2_1}, we have
\begin{equation*}
\begin{split}
%\label{ }
Det(\overline{N}^{s}(\widetilde{X}))= \beta^{s}\mathbf{ n}(\widetilde{X})\cdot \bar{\mathbf{ n}}\geqslant \beta^{s} (1- \bar{\kappa}), \;\;\; s=\pm ,
\end{split}
\end{equation*}
which means $\overline{N}^{s}(\widetilde{X})$ are non-singular when $h$ is small enough; hence, we can use them to form
\begin{equation}
\label{eq2_9}
\overline{M}^+(\widetilde{X})=(\overline{N}^-(\widetilde{X}))^{-1}\overline{N}^+(\widetilde{X}), \;\;\;\; \overline{M}^-(\widetilde{X})=(\overline{N}^+(\widetilde{X}))^{-1}\overline{N}^-(\widetilde{X}).
\end{equation}


% Then $\overline{M}^{\pm}$ has the following property:

\begin{lemma}
\label{lem2_6}
For the mesh $\mathcal{T}_h$ with $h$ sufficiently small, there exists a constant $C$ independent of interface location, $\beta^{\pm}$, and $\kappa$, such that, for two arbitrary points $\widetilde{X}_i, i = 1, 2$ on $\Gamma \cap T$, we have
\begin{equation}
\label{bounds_for_Mbar}
\| \overline{M}^-(\widetilde{X}_1)\| \leqslant C, ~~~ \| \overline{M}^+(\widetilde{X}_1)\| \leqslant \frac{C}{\rho}, ~~~ \| M^-(\widetilde{X}_2)\| \leqslant C, ~~~ \| M^+(\widetilde{X}_2)\| \leqslant \frac{C}{\rho},
\end{equation}
and
\begin{equation}
\label{eq2_13}
\| \overline{M}^-(\widetilde{X}_1)-M^-(\widetilde{X}_2)\| \leqslant C\kappa h, ~~~~ \| \overline{M}^+(\widetilde{X}_1)-M^+(\widetilde{X}_2)\| \leqslant \frac{C\kappa}{\rho}h.
\end{equation}
\end{lemma}
\begin{proof}
\eqref{bounds_for_Mbar} can be verified directly. We only prove \eqref{eq2_13} for the $-$ case and the arguments for the $+$ case are similar. For simplicity, we denote $\mathbf{ n}(\widetilde{X}_i)=(\tilde{n}_{ix}, \tilde{n}_{iy})$, $i=1,2$. Then by direct calculations, we have
\begin{equation*}
\begin{split}
\overline{M}^-(\widetilde{X}_1)-M^-(\widetilde{X}_2) =&
\left(\begin{array}{cc}
\tilde{n}_{1y}\bar{n}_{y}-\tilde{n}^2_{2y}+\rho(\tilde{n}_{1x}\bar{n}_{x}-\tilde{n}^2_{2x}) & (\rho-1)(\tilde{n}_{1y}\bar{n}_x-\tilde{n}_{2x}\tilde{n}_{2y}) \\
(\rho-1)(\tilde{n}_{1x}\bar{n}_y-\tilde{n}_{2x}\tilde{n}_{2y})  & \tilde{n}_{1x}\bar{n}_{x}-n^2_{2x}+\rho(\tilde{n}_{1y}\bar{n}_{y}-\tilde{n}^2_{2y})
\end{array}\right)\\
&+\frac{1-\bar{\mathbf{ n}}\cdot\mathbf{ n}(\widetilde{X}_1)}{\bar{\mathbf{ n}}\cdot\mathbf{ n}(\widetilde{X}_1)}
\left(\begin{array}{cc}
\tilde{n}_{1y}\bar{n}_{y}+\rho\tilde{n}_{1x}\bar{n}_{x} & (\rho-1)\tilde{n}_{1y}\bar{n}_x \\
(\rho-1)\tilde{n}_{1x}\bar{n}_y  & \tilde{n}_{1x}\bar{n}_{x}+\rho\tilde{n}_{1y}\bar{n}_{y}
\end{array}\right).
\end{split}
\end{equation*}
By the triangular inequality, \eqref{eq2_7a}, \eqref{eq2_8b}, \eqref{h_small_assump}, and $\rho\leqslant1$, we can verify that $\| \overline{M}^-(\widetilde{X}_1)-M^-(\widetilde{X}_2)\| \leqslant C\kappa h$ for a constant $C$ independent of interface location, $\beta^{\pm}$, and $\kappa$.
%By the triangular inequality, we have
%\begin{eqnarray}
%&&\| \overline{M}^-(\widetilde{X}_1) - M^-(\widetilde{X}_2) \| = \|(\overline{N}^+(\widetilde{X}_1))^{-1}\overline{N}^-(\widetilde{X}_1) - (N^+(\widetilde{X}_2))^{-1}N^-(\widetilde{X}_2)\| \nonumber \\
%&\leq& \|(\overline{N}^+(\widetilde{X}_1))^{-1}\| \, \| \overline{N}^-(\widetilde{X}_1) - N^-(\widetilde{X}_2)\| + \|(\overline{N}^+(\widetilde{X}_1))^{-1} - (N^+(\widetilde{X}_2))^{-1}\| \, \| N^-(\widetilde{X}_2) \|. \label{lem2_6_1}
%\end{eqnarray}
%By Lemma \ref{lem2_3} and Remark \ref{rem2_1}, we have
%\begin{equation}
%\label{lem2_6_2}
%\| \overline{N}^-(\widetilde{X}_1) - N^-(\widetilde{X}_2) \|
%=\norm{\left(\begin{array}{cc}
%\bar{n}_y-\tilde{n}_y(\widetilde{X}_2) & -\bar{n}_x + \tilde{n}_x(\widetilde{X}_2) \\
%\beta^-(\tilde{n}_x(\widetilde{X}_1)-\tilde{n}_x(\widetilde{X}_2)) & \beta^-(\tilde{n}_y(\widetilde{X}_1)-\tilde{n}_y(\widetilde{X}_2))
%\end{array}\right)} \leqslant Ch.
%\end{equation}
%In addition, we note that
%\begin{equation*}
%\begin{split}
%%\label{ }
%&(\overline{N}^+(\widetilde{X}_1))^{-1} - (N^+(\widetilde{X}_2))^{-1}\\
%=&\frac{1}{\beta^+\mathbf{ n}(\widetilde{X}_1)\cdot \bar{\mathbf{ n}}}
%\left(\begin{array}{cc}
%\beta^+\tilde{n}_y(\widetilde{X}_1) & \bar{n}_x \\
%-\beta^+\tilde{n}_x(\widetilde{X}_1) & \bar{n}_y
%\end{array}\right)-
%\frac{1}{\beta^+}
%\left(\begin{array}{cc}
%\beta^+\tilde{n}_y(\widetilde{X}_2) & \tilde{n}_x(\widetilde{X}_2) \\
%-\beta^+\tilde{n}_x(\widetilde{X}_2) & \tilde{n}_y(\widetilde{X}_2)
%\end{array}\right)\\
%=&\frac{1}{\beta^+\mathbf{ n}(\widetilde{X}_1)\cdot \bar{\mathbf{ n}}}
%\left(\begin{array}{cc}
%\beta^+(\tilde{n}_y(\widetilde{X}_1)-\tilde{n}_y(\widetilde{X}_2)) & (\bar{n}_x-\tilde{n}_x(\widetilde{X}_2)) \\
%-\beta^+(\tilde{n}_x(\widetilde{X}_1)-\tilde{n}_x(\widetilde{X}_2)) & (\bar{n}_y-\tilde{n}_y(\widetilde{X}_2))
%\end{array}\right) \\
%& \hspace{0.5in} + \frac{1-\mathbf{ n}(\widetilde{X}_1)\cdot \bar{\mathbf{ n}}}{\beta^+\mathbf{ n}(\widetilde{X}_1)\cdot \bar{\mathbf{ n}}}
%\left(\begin{array}{cc}
%\beta^+\tilde{n}_y(\widetilde{X}_2) & \tilde{n}_x(\widetilde{X}_2) \\
%-\beta^+\tilde{n}_x(\widetilde{X}_2) & \tilde{n}_y(\widetilde{X}_2)
%\end{array}\right).
%\end{split}
%\end{equation*}
%Then, applying Lemma \ref{lem2_3} and Remark \ref{rem2_1} to the identity above leads to
%\begin{eqnarray}
%\norm{(\overline{N}^+(\widetilde{X}_1))^{-1} - (N^+(\widetilde{X}_2))^{-1}} \leqslant Ch. \label{lem2_6_3}
%\end{eqnarray}
%Finally, applying estimates \eqref{lem2_6_2} and \eqref{lem2_6_3} to \eqref{lem2_6_1} yields the desired estimate \eqref{eq2_13}.
\end{proof}

The following lemmas provide a group of identities on interface elements.
\begin{lemma}
\label{lem2_4}
For the mesh $\mathcal{T}_h$ with $h$ sufficiently small, the following results hold for all $\widetilde{X} \in \Gamma$:
\begin{itemize}
\item
$\overline{M}^-(\widetilde{X})$ and $\overline{M}^+(\widetilde{X})$ are inverse matrices to each other, i.e.,
\begin{equation}
\label{eq2_10}
\overline{M}^-(\widetilde{X})\overline{M}^+(\widetilde{X})=I,~~\overline{M}^+(\widetilde{X})\overline{M}^-(\widetilde{X}) = I.
\end{equation}

\item
Matrix $(\overline{M}^-(\widetilde{X}))^T$ has two eigenvalues $1$ and $\rho$ with the corresponding eigenvectors $\bar{\mathbf{ t}}$ and $\mathbf{ n}(\widetilde{X})$, i.e.,
\begin{equation}
\label{eq2_12}
\left( \overline{M}^-(\widetilde{X}) \right)^T \bar{\mathbf{ t}}=\bar{\mathbf{ t}},~~\left( \overline{M}^- (\widetilde{X})\right)^T \mathbf{ n}(\widetilde{X})=\rho\mathbf{ n}(\widetilde{X}).
\end{equation}

\item
Similarly, matrix $(\overline{M}^+(\widetilde{X}))^T$ has two eigenvalues $1$ and $1/\rho$ with the corresponding eigenvectors $\bar{\mathbf{ t}}$ and $\mathbf{ n}(\widetilde{X})$, respectively, i.e.,
\begin{equation}
\label{eq2_11}
\left( \overline{M}^+(\widetilde{X}) \right)^T \bar{\mathbf{ t}}=\bar{\mathbf{ t}},~~\left( \overline{M}^+(\widetilde{X}) \right)^T \mathbf{ n}(\widetilde{X})=\frac{1}{\rho}\mathbf{ n}(\widetilde{X}).
\end{equation}

\end{itemize}
\end{lemma}
\begin{proof} First it is easy to see that $\overline{M}^-(\widetilde{X})\overline{M}^+(\widetilde{X})=(\overline{N}^-(\widetilde{X}))^{-1}\overline{N}^+(\widetilde{X})(\overline{N}^+(\widetilde{X}))^{-1}\overline{N}^-(\widetilde{X})=I$. Next by direct calculation, we have
\begin{equation*}
\begin{split}
%\label{ }
\overline{M}^-(\widetilde{X})
&=\frac{1}{\bar{n}_x\tilde{n}_x(\widetilde{X})+\bar{n}_y\tilde{n}_y(\widetilde{X})}
\left(\begin{array}{cc}
\bar{n}_y\tilde{n}_y(\widetilde{X})+\rho\bar{n}_x\tilde{n}_x(\widetilde{X}) & -\bar{n}_x\tilde{n}_y(\widetilde{X})+\rho\bar{n}_x\tilde{n}_y(\widetilde{X}) \\
-\bar{n}_y\tilde{n}_x(\widetilde{X})+\rho\bar{n}_y\tilde{n}_x(\widetilde{X}) & \bar{n}_x\tilde{n}_x(\widetilde{X})+\rho\bar{n}_y\tilde{n}_y(\widetilde{X})
\end{array}\right)
\end{split}
\end{equation*}
from which we can easily verify that $(\overline{M}^-(\widetilde{X}))^{T}\bar{\mathbf{ t}}=\bar{\mathbf{ t}}$ and
$ (\overline{M}^-(\widetilde{X}))^{T}\mathbf{ n}(\widetilde{X})=\rho \; \mathbf{ n}(\widetilde{X})$.
The results about $(\overline{M}^+(\widetilde{X}))^{T}$ follow from the fact $(\overline{M}^-(\widetilde{X}))^{T}(\overline{M}^+(\widetilde{X}))^{T}=I$.
\end{proof}

\begin{lemma}
\label{lem2_5}
Let $\mathcal{T}_h$ be a mesh with $h$ sufficiently small. Let $P\in \Omega$ and $\widetilde{X}$ be an arbitrary point on $\Gamma \cap T$. Then the following vectors are independent of $\overline{X}\in l$:
\begin{equation*}\label{}
(\overline{M}^+(\widetilde{X}) - I)^T(P - \overline{X}) \;\;\; and \;\;\; (\overline{M}^-(\widetilde{X}) - I)^T(P - \overline{X}).
\end{equation*}
\end{lemma}
\begin{proof}
For two arbitrary points $\overline{X}_i \in l, i = 1, 2$, $\overline{X}_1 - \overline{X}_2$ is a scalar multiple of $\bar{\mathbf{ t}}$. Hence, by Lemma \ref{lem2_4},
\begin{eqnarray*}
(\overline{M}^-(\widetilde{X}) - I)^T(P - \overline{X}_1) - (\overline{M}^-(\widetilde{X}) - I)^T(P - \overline{X}_2) = (\overline{M}^-(\widetilde{X}) - I)^T(\overline{X}_1 - \overline{X}_2) = 0
\end{eqnarray*}
which leads to $(\overline{M}^-(\widetilde{X}) - I)^T(P - \overline{X}_1) = (\overline{M}^-(\widetilde{X}) - I)^T(P - \overline{X}_2)$. Therefore
$(\overline{M}^-(\widetilde{X}) - I)^T(P - \overline{X})$ does not change when $\overline{X}\in l$ varies. The result for  $(\overline{M}^+(\widetilde{X}) - I)^T(P - \overline{X})$ can be proven similarly.
\end{proof}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%section3
\section{Multipoint Taylor Expansions on Interface Elements}
In this section, extending those in \cite{2009HeTHESIS,2008HeLinLin,2004LiLinLinRogers,1982Xu,2013ZhangTHESIS}, we derive multipoint Taylor expansions in more general formats for a function $u\in PC^2_{int}(T)$ over an arbitrary interface element $T\in\mathcal{T}^i_h$, in which $u(M_i), i \in \mathcal{I}$ is described in terms of $u$ and its derivatives at $x \in T^s, s = \pm$.  We also estimate the remainders in these expansions.
And as in \cite{2008HeLinLin}, we call a point $X\in T$ an obscure point if one of the lines $\overline{XM_i}, 1 \leq i \leq DOF(T)$ can intersect $\Gamma$ more than once. To facilitate a clear expository presentation of main ideas in our analysis, we carry out error estimation only for interface elements without any obscure points.
For the case containing obscure points, we can use a first order expansion for $u$ and use the argument that the measure of obscure points is bounded by $\mathcal{O}(h^3)$.


%Our discussion here focuses on $T^-$, and all the results and related arguments are extendable to $T^+$.

First, we partition $\mathcal{I}$ into two index sets: $\mathcal{I}^+=\{ i: M_i\in T^+ \}$ and $\mathcal{I}^-=\{ i: M_i\in T^- \}$ according to the locations of
$M_i, i \in \mathcal{I}$. For every $X\in T$, we let $Y_i(t,X)=tM_i+(1-t)X$. When
$X$ and $M_i$ are on different sides of $\Gamma$, we let $\tilde{t}_i=\tilde{t}_i(X)\in [0,1]$ such that $\widetilde{Y}_i=Y_i(\tilde{t}_i,X)$ is on the curve $\Gamma\cap T$, see Figure \ref{fig3_1} for an illustration in which the rotated-$Q_1$ finite elements are considered. When $M_i$ and $X$ are on the same side of $\Gamma$, by the standard second order Taylor expansion of $u\in PC^2_{int}(T)$, we have
\begin{eqnarray}
u^s(M_i)&=&u^s(X)+ \nabla u^s(X)\cdot(M_i-X) + R_i^s(X),~ i\in \mathcal{I}^s, s = \pm, ~\forall X \in T^s, \label{eq3_1}\\
\text{with~~} R_i^s(X) &=&\int_0^1(1-t)\frac{d^2}{dt^2}u^s(Y_i(t,X))dt, ~i\in \mathcal{I}^s, ~\forall X \in T^s. \label{eq3_2}
\end{eqnarray}

%graph
\begin{figure}[H]
\centering
\includegraphics[width=2.3in]{interf_elem_nonconf_midpoint_mult_expansion.pdf}
\caption{ The expansion of $u$ in a rectangular interface element}
\label{fig3_1}
\end{figure}

In the following discussion, we denote $s = \pm, s' = \mp$, i.e., $s$ and $s'$ take opposite signs whenever a formula have them both.
When, $M_i$ and $X$ are on different sides of $\Gamma$, the expansions in \cite{2009HeTHESIS,2008HeLinLin,2004LiLinLinRogers,2013ZhangTHESIS} can be
generalized to the following format for $u\in PC^2_{int}(T)$:
\begin{equation}
\begin{split}
\label{eq3_3}
u^{s'}(M_i)=&u^{s}(X)+\nabla u^{s}(X)\cdot (M_i-X)+ \left(\left( M^{s}(\widetilde{Y}_i)-I \right)\nabla u^{s}(X)\right) \cdot (M_i-\widetilde{Y}_i)\\
& + R_i^{s}(X), ~i\in \mathcal{I}^{s'}, ~\forall X \in T^{s},
%&+(1-\tilde{t}_i)\int_{0}^{\tilde{t}_i}\frac{d}{dt} \left( (M^-(\widetilde{Y}_i)-I)\nabla u^-(Y_i(t,X))\cdot (M_i-X) \right)dt\\
%&+\int_{0}^{\tilde{t}_i}(1-t)\frac{d^2u^-}{dt^2}(Y_i(t,X))dt+\int_{\tilde{t}_i}^{1}(1-t)\frac{d^2u^+}{dt^2}(Y_i(t,X))dt,
\end{split}
\end{equation}
with
\begin{equation}
\begin{cases} \label{eq3_4}
&R_i^s(X)=R_{i1}^s(X)+R_{i2}^s(X)+R_{i3}^s(X), ~i \in \mathcal{I}^{s'}, X \in T^s, \\
&R_{i1}^s(X)=\int_{0}^{\tilde{t}_i}(1-t)\frac{d^2u^s}{dt^2}(Y_i(t,X))dt, \;\;\; R_{i2}^s(X)=\int_{\tilde{t}_i}^{1}(1-t)\frac{d^2u^{s'}}{dt^2}(Y_i(t,X))dt, ~i \in \mathcal{I}^{s'}, X \in T^s,\\
&R_{i3}^s(X)=(1-\tilde{t}_i)\int_{0}^{\tilde{t}_i}\frac{d}{dt} \left( (M^s(\widetilde{Y}_i)-I)\nabla u^s(Y_i(t,X))\cdot (M_i-X) \right)dt, ~i \in \mathcal{I}^{s'}, X \in T^s,
\end{cases}
\end{equation}
where $M^s(\widetilde{Y}_i)$ are from \eqref{eq2_3} and \eqref{eq2_4}. We proceed to estimate remainders in \eqref{eq3_1} and \eqref{eq3_3}.

%%If $i\in \mathcal{I}^+$, the expansion of $u^+(M_i)$ has to be constructed piecewisely in the following lemma.
%\begin{lemma}
%\label{lem3_1}
%For $u\in PC^2_{int}(T)$ we have the following expansion:
%%\begin{subequations}
%%\label{eq3_3_total}
%\begin{equation}
%\begin{split}
%\label{eq3_3}
%u^{s'}(M_i)=&u^{s}(X)+\nabla u^{s}(X)\cdot (M_i-X)+ \left(\left( M^{s}(\widetilde{Y}_i)-I \right)\nabla u^{s}(X)\right) \cdot (M_i-\widetilde{Y}_i)\\
%& + R_i^{s}(X), ~i\in \mathcal{I}^{s'}, ~\forall X \in T^{s},
%%&+(1-\tilde{t}_i)\int_{0}^{\tilde{t}_i}\frac{d}{dt} \left( (M^-(\widetilde{Y}_i)-I)\nabla u^-(Y_i(t,X))\cdot (M_i-X) \right)dt\\
%%&+\int_{0}^{\tilde{t}_i}(1-t)\frac{d^2u^-}{dt^2}(Y_i(t,X))dt+\int_{\tilde{t}_i}^{1}(1-t)\frac{d^2u^+}{dt^2}(Y_i(t,X))dt,
%\end{split}
%\end{equation}
%%\begin{equation}
%%\begin{split}
%%\label{eq3_3b}
%%u^-(M_i)=&u^+(X)+\nabla u^+(X)\cdot (M_i-X)+ \left(\left( M^+(\widetilde{Y}_i)-I \right)\nabla u^+(X)\right) \cdot (M_i-\widetilde{Y}_i)\\
%%& + R_i^+(X), ~i\in \mathcal{I}^-, ~\forall X \in T^+,
%%%&+(1-\tilde{t}_i)\int_{0}^{\tilde{t}_i}\frac{d}{dt} \left( (M^-(\widetilde{Y}_i)-I)\nabla u^-(Y_i(t,X))\cdot (M_i-X) \right)dt\\
%%%&+\int_{0}^{\tilde{t}_i}(1-t)\frac{d^2u^-}{dt^2}(Y_i(t,X))dt+\int_{\tilde{t}_i}^{1}(1-t)\frac{d^2u^+}{dt^2}(Y_i(t,X))dt,
%%\end{split}
%%\end{equation}
%%\end{subequations}
%with
%\begin{equation}
%\begin{cases} \label{eq3_4}
%&R_i^s(X)=R_{i1}^s(X)+R_{i2}^s(X)+R_{i3}^s(X), ~i \in \mathcal{I}^{s'}, X \in T^s, \\
%&R_{i1}^s(X)=\int_{0}^{\tilde{t}_i}(1-t)\frac{d^2u^s}{dt^2}(Y_i(t,X))dt, \;\;\; R_{i2}^s(X)=\int_{\tilde{t}_i}^{1}(1-t)\frac{d^2u^{s'}}{dt^2}(Y_i(t,X))dt, ~i \in \mathcal{I}^{s'}, X \in T^s,\\
%&R_{i3}^s(X)=(1-\tilde{t}_i)\int_{0}^{\tilde{t}_i}\frac{d}{dt} \left( (M^s(\widetilde{Y}_i)-I)\nabla u^s(Y_i(t,X))\cdot (M_i-X) \right)dt, ~i \in \mathcal{I}^{s'}, X \in T^s,
%\end{cases}
%\end{equation}
%where $M^s(\widetilde{Y}_i)$ are from \eqref{eq2_3} and \eqref{eq2_4}.
%\end{lemma}
%\begin{proof}
%dfasdfasdfsadfaf?????
%\end{proof}



\begin{lemma}
\label{lem3_2}
Assume $u\in PC^2_{int}(T)$. Then there exist constants $C>0$ independent of the interface location and $\beta^{\pm}$ such that
\begin{equation}
\label{eq3_5}
   \int_{T^s}(1-t)^2|u_{d_1d_2}(Y_i(t,X))|^2dX \leqslant C|u|^2_{2,T},~s = \pm, ~\forall t \in [0, 1],
\end{equation}
where $d_1, d_2 =x$ or $y$.
\end{lemma}
\begin{proof}
Let $M_i = (x_i, y_i), i \in \mathcal{I}$ and let $\xi=t x_i+(1-t)x$ and $\eta=t y_i+(1-t)y$.
For each fixed $t \in [0,1]$, define
\begin{equation}
%\label{ }
T^s(t)=\{ tM_i+(1-t)X ~|~ X\in T^s \}.
\end{equation}
Since $T^s(t) \subseteq T, s = \pm$, we have
\begin{eqnarray*}
 \int_{T^s}(1-t)^2|u_{d_1d_2}(Y(t,X))|^2dX
= \int_{T^s(t)}(1-t)^2 u^2_{d_1d_2}(\xi,\eta)(1-t)^{-2}d\xi d\eta \leqslant |u|^2_{2,T}
\end{eqnarray*}
which leads to \eqref{eq3_5}. \end{proof}
%
%If $i\in \mathcal{I}^-$, $T^-(t)$ must be the subset of $T^-$. However if $i\in \mathcal{I}^+$, $T^-(t)$ might not be the subset of $T^-$, since $\tilde{t}_i=\tilde{t}_i(X)$ depends on $X$. But $T^-(t)$ must be the subset of $T$. So we need to estimate the integral on the whole element $T$. Actually it is not hard to see that


By a direct calculation we have
\begin{equation}
\begin{split}
\label{eq3_7}
\frac{d^2}{dt^2}u^s(Y_i(t,X))=(M_i-X)^TH_u^s(Y_i(t,X))(M_i-X), ~s = \pm,
\end{split}
\end{equation}
where
\begin{equation*}
%\label{ }
H_u^s(Y_i(t,X))=
\left(\begin{array}{cc}
u_{xx}^s(Y_i(t,x)) & u_{xy}^s(Y_i(t,x)) \\
u_{yx}^s(Y_i(t,x)) & u_{yy}^s(Y_i(t,x))
\end{array}\right), ~s = \pm,
\end{equation*}
is the Hessian matrix of $u^s$. We are now ready to derive bounds for the remainders in the following lemmas.

\begin{lemma}
\label{lem3_3}
Assume $u\in PC^2_{int}(T)$, there exist constants $C>0$ independent of the location of the interface and $\beta^{\pm}$ such that
\begin{eqnarray}
&&\| R_i^s \|_{0,T^s} \leqslant Ch^2 |u|_{2,T}, \;\;\; i\in \mathcal{I}^s, ~s = \pm, \label{eq3_8} \\
&&\| R_{i1}^s \|_{0,T^s} \leqslant Ch^2 |u|_{2,T}, \;\;\; \| R_{i2}^s \|_{0,T^s}  \leqslant Ch^2 |u|_{2,T}, \;\;\; i \in \mathcal{I}^{s'}, ~s = \pm. \label{eq3_9}
\end{eqnarray}
\end{lemma}
\begin{proof}
According to \eqref{eq3_2} and \eqref{eq3_7}, for $s = \pm$, we have
\begin{equation*}
\begin{split}
%\label{ }
\| R_i^s \|_{0,T^s}=& \left(  \int_{T^s} \left( \int_{0}^{1}(1-t)(M_i-X)H_u^s(Y_i(t,X))(M_i-X)dt \right)^2 dX  \right)^{\frac{1}{2}}\\
%\leqslant& Ch^2 \left( \int_{T^-} \left( \int_0^1(1-t)(u_{xx}(Y_i(t,X))+u_{xy}(Y_i(t,X))+u_{yy}(Y_i(t,X)))dt\right)^2dX \right)^{\frac{1}{2}}\\
%\leqslant& Ch^2 \int_{0}^{1} \left( \int_{T^-}(1-t)^2|u_{xx}(Y_i(t,X))+u_{xy}(Y_i(t,X))+u_{yy}(Y_i(t,X))|^2dX\right)^{\frac{1}{2}}dt \\
\leqslant& Ch^2  \int_{0}^{1} \left( \int_{T^s}(1-t)^2(|u_{xx}^s(Y_i(t,X))|^2+|u_{xy}^s(Y_i(t,X))|^2+|u_{yy}^s(Y_i(t,X))|^2)dX\right)^{\frac{1}{2}}dt.
\end{split}
\end{equation*}
Then \eqref{eq3_8} follows from Lemma \ref{lem3_2}. Estimate \eqref{eq3_9} can be derived similarly.
\end{proof}

%\begin{lemma}
%\label{lem3_4}
%Assume $u\in PC^2_{int}(T)$, then there exists constants $C>0$ independent of the interface location such that
%\begin{equation}
%\label{eq3_9}
%   \| R_{i1}^s \|_{0,T^s} \leqslant Ch^2 |u|_{2,T}, \;\;\; \| R_{i2}^s \|_{0,T^s}  \leqslant Ch^2 |u|_{2,T}, \;\;\; i \in \mathcal{I}^{s'}.
%\end{equation}
%\end{lemma}
%\begin{proof}
%The argument is similar to Lemma \ref{lem3_3}.
%\end{proof}

Moreover, note that for $i\in\mathcal{I}^{s'}, X \in T^s$ and $t\in [0, \tilde{t}_i(X)]$, it can be verified that
\begin{equation}
\begin{split}
\label{eq3_10}
\frac{d}{dt}\left( (M^{s}(\widetilde{Y}_i)-I)\nabla u^s(Y_i(t,X))\cdot(M_i-X) \right)= (M_i-X)^T H_u^s(Y_i(t,X))(M^{s}(\widetilde{Y}_i)-I)^T(M_i-X).
\end{split}
\end{equation}

\begin{lemma}
\label{lem3_5}
Assume $u\in PC^2_{int}(T)$, there exist constants $C>0$ independent of the interface location and $\beta^{\pm}$ such that
\begin{equation}
\label{eq3_11}
\| R_{i3}^- \|_{0,T^-}\leqslant Ch^2 |u|_{2,T}, \;\;\; i\in \mathcal{I}^{+}, ~~~~ \| R_{i3}^+ \|_{0,T^+}\leqslant \frac{C}{\rho}h^2 |u|_{2,T}, \;\;\; i\in \mathcal{I}^{-}.
\end{equation}
\end{lemma}
\begin{proof}
We only provide the proof of $\| R_{i3}^- \|_{0,T^-}$ and the argument for $\| R_{i3}^+ \|_{0,T^+}$ is similar. According to \eqref{eq3_10} and the fact $0\leqslant1-\tilde{t}_i(X)\leqslant 1-t$ for any $t\in [0,\tilde{t}_i(X)]$, we have
\begin{equation*}
\begin{split}
%\label{ }
\|R^-_{i3}\|_{0,T^-}
=&\left( \int_{T^-} (1-\tilde{t}_i(X))^2 \left( \int_{0}^{\tilde{t}_i(X)}  (M_i-X)^T H_u^-(Y_i(t,X))(M^{-}(\widetilde{Y}_i)-I)^T(M_i-X)dt \right)^2 dX \right)^{\frac{1}{2}}\\
\leqslant& Ch^2 \left( \int_{T^-} \left( \int_0^1|1-t| \; \big(\abs{u^-_{xx}(Y_i(t,X))} + \abs{u^-_{xy}(Y_i(t,X))} + \abs{u^-_{yy}(Y_i(t,X))}\big)dt\right)^2dX \right)^{\frac{1}{2}}\\
\leqslant& Ch^2 \int_{0}^{1} \left( \int_{T^-}(1-t)^2\big(|u^-_{xx}(Y_i(t,X))|^2+|u^-_{xy}(Y_i(t,X))|^2+|u^-_{yy}(Y_i(t,X))|^2\big)dX\right)^{\frac{1}{2}}dt,
\end{split}
\end{equation*}
where we use the fact $\rho\leqslant1$. Then applying the estimates in Lemma \ref{lem3_2} to the above leads to \eqref{eq3_11}.
\end{proof}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%section 5
\section{IFE Spaces and Their Properties}
\label{sec:IFE Spaces and Their Properties}

In this section, we discuss IFE spaces constructed from the related finite elements $(T,\Pi_T,\Sigma^P_T)$ for $T \in \mathcal{T}_h$
described in \eqref{eq4_1} and \eqref{eq4_2}. We will first address the unisolvence of the immersed finite elements on interface elements. We will then present a few fundamental properties of IFE functions. Moreover, we will show that these IFE spaces have the optimal approximation capability according to the polynomials used to constructed them.

\subsection{Local IFE spaces}

%In this section, we are going to derive the IFE shape functions with the the degrees of freedom imposed at mid points in interface elements and study the approximation capabilities for the related IFE space. Recall the nonconforming $Q_1$ finite elements $(T,\Pi_T,\Sigma^P_T)$ from \cite{1978Ciarlet}:
%%\begin{equation}
%%\label{eq4_1}
%%\Pi_T=\textrm{Span} \{ 1,x,y,x^2-y^2 \},
%%\end{equation}
%%\begin{equation}
%%\label{eq4_2}
%%\Sigma_T^P= \{ \psi^P_T(M_i): i \in \mathcal{I}, \; \forall \psi^P_T \in \Pi_T \}.
%%\end{equation}
%For any $\psi^P_T\in \Pi_T$, denote the coefficient of $x^2-y^2$ by $d(\psi^P_T)$. And on each element $T\in\mathcal{T}_h$, the basis functions for $(T,\Pi_T,\Sigma^P_T)$
%are denoted by $\psi^P_{i,T}$, $i\in\mathcal{I}$, which have the following properties
%%\begin{equation}
%%\label{eqn_psi_iT_values_bounds}
%%\psi^P_{i,T}(M_j)=\delta_{ij},~~\norm{\psi^P_{i,T}}_{\infty, T} \leq C,~~\norm{\nabla \psi^P_{i,T}}_{\infty, T} \leq Ch^{-1}, \;\; i,j\in \mathcal{I},
%%\end{equation}
%%where $\delta_{ij}$ is the \textit{Kronecker} delta function.


First, on each element $T \in \mathcal{T}_h$, the standard finite element $(T,\Pi_T,\Sigma^P_T)$ leads to the following local finite element space:
\begin{equation}
\label{eq4_3}
S^{P}_h(T)= \textrm{Span} \{ \psi^P_{i,T}:\; i\in\mathcal{I} \},
\end{equation}
where $\psi^P_{i,T}, i\in\mathcal{I}$ are the shape functions satisfying \eqref{eqn_psi_iT_values_bounds}. This local finite element space is then naturally used as the local
IFE space on every non-interface element $T \in \mathcal{T}_h^n$. Therefore, our effort here focuses on the local IFE space on interface elements. We will discuss the unisolvence, i.e., we will show that the local degrees of freedom $\Sigma^P_T$ can uniquely determine an IFE function with a suitable set of interface jump conditions. The unisolvence guarantees the existence and uniqueness of IFE shape functions that can span the local IFE space on interface elements.

Let $T\in \mathcal{T}^i_h$ be a typical interface element with vertices $A_i, i \in \mathcal{I}$. Without loss of generality, we assume
\begin{eqnarray}\label{eq4_4}
\begin{array}{ll}
A_1=(0,0),~~A_2=(h,0),~~A_3=(h,h), & \text{($T$ is a triangular element)}, \\
A_1=(0,0),~~A_2=(h,0),~~A_3=(h,h),~~A_4=(0,h), & \text{($T$ is a rectangular element)},
\end{array}
\end{eqnarray}
and the edges of $T$ are denoted as
\begin{eqnarray}\label{eq4_5}
\begin{array}{ll}
b_1=\overline{A_1A_2},~~b_2=\overline{A_2A_3},~~b_3=\overline{A_3A_1}, & \text{($T$ is a triangular element)}, \\
b_1=\overline{A_1A_2},~~b_2=\overline{A_2A_3},~~b_3=\overline{A_3A_4},~~b_4=\overline{A_4A_1}, & \text{($T$ is a rectangular element)}.
\end{array}
\end{eqnarray}

%we consider the following typical element $T$
%\begin{equation}
%\label{eq4_4}
%A_1=(0,0), \;\;\; A_2=(h,0), \;\;\; A_3=(h,h), \;\;\; A_4=(0,h),
%\end{equation}
%with the four edges:
%\begin{equation}
%\label{eq4_5}
%b_1=\overline{A_1A_2}, \;\;\; b_2=\overline{A_2A_3}, \;\;\; b_3=\overline{A_3A_4}, \;\;\; b_4=\overline{A_4A_1}. \;\;\;
%\end{equation}

On each interface element $T$, we consider IFE functions in the following piecewise polynomial format:
\begin{equation}
\label{eq4_10}
\phi^P_T(X) =
\left\{
\begin{aligned}
\phi^{P\;-}_T(X)=\phi^{-}_T(X)\in \Pi_T \;\;\;\;\; & \text{if} \;\; X\in T^-, \\
\phi^{P\;+}_T(X)=\phi^{+}_T(X)\in \Pi_T \;\;\;\;\; & \text{if} \;\; X\in T^+,
\end{aligned}
\right.
\end{equation}
such that it can satisfy the jump conditions \eqref{eq1_3} and \eqref{eq1_4} in an approximate sense as follow:
\begin{eqnarray}
%&&\phi^{-}_T\in \Pi_{T}, \;\;\; and \;\;\;  \phi^{+}_T\in \Pi_{T}, \label{eq4_6} \\
&&\begin{cases}
\phi^{-}_T|_{l}=\phi^{+}_T|_{l}, &\text{($T$ is a triangular element)}, \\
\phi^{-}_T|_{l}=\phi^{+}_T|_{l}, ~~d(\phi^{-}_T)=d(\phi^{+}_T), &\text{($T$ is a rectangular element)},
\end{cases} \label{eq4_8} \\
&&\beta^-\nabla\phi^{-}_T(F)\cdot\mathbf{ n}(F)=\beta^+\nabla\phi^{+}_T(F)\cdot\mathbf{ n}(F), \label{eq4_9}
\end{eqnarray}
where $d(p)$ denotes the coefficient in the second degree term for $p \in \Pi_T$ and $F$ is an arbitrary point on $\Gamma \cap T$.
For an IFE function $\phi^P_T$ such that
\begin{equation}
\label{eq4_12}
\phi^P_T(M_i)=v_i, \;\;\;\; i\in \mathcal{I},
\end{equation}
we can first expand $\phi^P_T$ on the sub-element with more degrees of freedom, i.e., on $T^+$ with the assumption that
$\abs{\mathcal{I}^+} \geq \abs{\mathcal{I}^-}$ without loss of generality, and the condition \eqref{eq4_8} then implies that
\begin{equation}
\label{eq4_13}
\phi^P_T(X) =
\begin{cases}
 \phi^{P\;-}_T(X)  = \phi^{P\;+}_T(X)+c_0L(X) & \text{if} \;\; X\in T^-, \\
 \phi^{P\;+}_T(X)  = \sum_{i\in\mathcal{I}^-}c_i\psi^P_{i,T}(X)+\sum_{i\in\mathcal{I}^+}v_i\psi^P_{i,T}(X)& \text{if} \;\; X\in T^+,
\end{cases}
\end{equation}
where the function
\begin{equation}
\label{eq4_11}
L(X)=\bar{\mathbf{ n}}\cdot(X-D)
\end{equation}
is such that $L(X)=0$ is the equation of the line $l$ and $\nabla L(X) = \bar{\mathbf{n}}$.

Recall from Remark \ref{rem2_1},
$
\nabla L(F)\cdot\mathbf{ n}(F)=\bar{\mathbf{n}}\cdot\mathbf{ n}(F)\geqslant 1-\bar{\kappa} >0,
$
when $h$ is small enough; hence
\begin{eqnarray}
\label{k}
\mu = \left( \frac{1}{\rho}-1 \right)\frac{ 1}{\bar{\mathbf{n}}\cdot\mathbf{ n}(F)}
\end{eqnarray}
is well defined, and, by $\rho\leqslant1$, we have
\begin{eqnarray}
\label{bound_for_k}
0\leqslant \mu  \leqslant \left( \frac{1}{\rho}-1 \right) \frac{1}{1 - \bar{\kappa}}.
\end{eqnarray}
By condition \eqref{eq4_9}, we then have
\begin{equation}
\label{eqn4_18}
c_0=
\mu \left( \sum_{i\in\mathcal{I}^-}c_i\nabla\psi^P_{i,T}(F)\cdot\mathbf{ n}(F)+\sum_{i\in\mathcal{I}^+}v_i\nabla\psi^P_{i,T}(F)\cdot\mathbf{ n}(F) \right).
\end{equation}
Putting this formula for $c_0$ in formula \eqref{eq4_13} for $\phi_T^P(X)$ and setting
$\phi^{P\;-}_T(M_j)=v_j$ for $j \in \mathcal{I}^-$ leads to the following linear system for $c_i, i \in \mathcal{I}^-$:
\begin{equation}
\begin{split}
\label{eqn4_20}
&\sum_{i\in\mathcal{I}^-} \left( \psi^P_{i,T}(M_j) + \mu \nabla\psi^P_{i,T}(F)\cdot\mathbf{ n}(F) L(M_j) \right)c_i \\
=& v_j-\sum_{i\in\mathcal{I}^+} \left( \psi^P_{i,T}(M_j) + \mu \nabla\psi^P_{i,T}(F)\cdot\mathbf{ n}(F) L(M_j) \right)v_i,~~j \in \mathcal{I}^-.
\end{split}
\end{equation}
Since $\psi^P_{i,T}(M_j)=\delta_{ij}$, for $i,j\in\mathcal{I}^-$, we can write the linear system \eqref{eqn4_20} in the following matrix form:
\begin{equation}
\label{eqn4_15}
(I + \mu \,\bfdelta \bfgamma^T )\bfc = \bfb,
\end{equation}
where $\mathbf{c}=(c_i)_{i\in \mathcal{I}^-}$,
\begin{equation}
\label{eqn4_16}
\bfgamma = \left( \nabla\psi^P_{i,T}(F)\cdot\mathbf{ n}(F) \right)_{i\in\mathcal{I}^-}, \;\;\; \bfdelta=\left( L(M_i) \right)_{i\in\mathcal{I}^-}
\end{equation}
and
\begin{equation}
\label{eqn4_17}
\bfb =
\left(
v_i- \mu L(M_i) \sum_{j\in\mathcal{I}^+} \nabla\psi^P_{j,T}(F)\cdot\mathbf{ n}(F)v_j
 \right)_{i\in\mathcal{I}^-}
\end{equation}
are all column vectors. We proceed to show that $\phi_T^P(X)$, i.e., its coefficients $c_0, \bfc$ are uniquely determined. We need the following two lemmas. Let
$\bar{\bfgamma}=\left( \nabla\psi_{i,T}(F_{\bot})\cdot \bar{\bfn} \right)_{i\in\mathcal{I}^-}$.
\begin{lemma}
\label{lemma4_3}
For all the interface elements, we have
$
\bar{\bfgamma}^T\bfdelta \in [0,1]
$.
And for the linear and bilinear $\psi_{i, T}, i\in\mathcal{I}$, there holds
\begin{equation}
\label{eq_Ldelgam_1}
|L(M_i)|<2h \, \sqrt{\bar{\bfgamma}^T\bfdelta}, ~~~ \forall i\in\mathcal{I}^-.
\end{equation}
Furthermore, for the bilinear $\psi_{i, T}, i\in\mathcal{I}$, if $F$ is chosen to be such that $F_{\bot}$ is the mid-point of the line $\overline{DE}$, then
\begin{equation}
\label{eq_Ldelgam_2}
|L(M_i)|<2h \, \bar{\bfgamma}^T\bfdelta, ~~~ \forall i\in\mathcal{I}^-.
\end{equation}
\end{lemma}
\begin{proof}
We only give the proof for the case in which $\Pi_T$ is the rotated-$Q_1$ polynomial space, the interface element $T$ is such that $\mathcal{I}^-=\{ 1 \}$ and $\mathcal{I}^+=\{2,3,4\}$ with $D=(hd, 0)$ and $E=(0,he)$ for some $d\in[1/2,1]$ and $e\in[0,1/2]$. Similar arguments apply to all other cases. First, $\bar{\mathbf{ n}}=(e,d)/\sqrt{d^2+e^2}$. Hence, $F_{\bot}=(td,e(h-t))$ for some $t\in (0,h)$. By direct calculation, we have
\begin{equation}
\begin{split}
\label{}
 \bar{\bfgamma}^T \bfdelta= \frac{1}{h(d^2+e^2)}\left[ h(e-2d)-2(2t-h)de \right] (\frac{1}{2}-d)e,
\end{split} \label{eq_lemma5.1_1}
\end{equation}
Note that
\begin{equation*}
\label{}
\bar{\bfgamma}^T \bfdelta = \frac{e-2d+2de}{d^2+e^2}(\frac{1}{2}-d)e\in [0,1] \;\;\text{if} \;\; t=0 \;\;\; \text{and} \;\;\;
\bar{\bfgamma}^T \bfdelta = \frac{e-2d-2de}{d^2+e^2}(\frac{1}{2}-d)e\in [0,1] \;\;\text{if} \;\; t=h,
\end{equation*}
which leads to $\bar{\bfgamma}^T\bfdelta \in [0,1]$ because $\bar{\bfgamma}^T \bfdelta$ is a linear function of $t$ according to \eqref{eq_lemma5.1_1}. And \eqref{eq_Ldelgam_1} and \eqref{eq_Ldelgam_2} follow from similar calculation.
\end{proof}

%\begin{lemma}
%\label{lemma_Ldelgam}
%For linear and bilinear elements, given all interface elements, there exist constants $C$ independent of interface location, such that
%\begin{equation}
%\label{eq_Ldelgam}
%|L(M_i)|<Ch \, \bar{\bfdelta}\bfgamma, ~~~ \forall i\in\mathcal{I}^-.
%\end{equation}
%\end{lemma}
%\begin{proof}
%%For simplicity, we only give the proof for the case discussed in Lemma \ref{lemma4_3}. Lemma \ref{lemma4_3} shows that $\bar{\bfdelta}\bfgamma$ is monotony and linear with respect to $t$. So we have
%%$$
%%\frac{L(M_1)}{\bar{\bfdelta}\bfgamma}=
%%$$
%The inequalities can be verified directly by similar arguments as the Lemma \ref{lemma4_3}.
%\end{proof}

\begin{lemma}
\label{lemma4_4}
For $h$ small enough, we have
\begin{equation}
\label{eqn4_22}
%1+k\,\bfgamma^T\bfdelta \geq % \min{ \left(1, \frac{1}{\rho}-Ch^2 \right) }-Ch.
%\min{ \left(1, \frac{1}{\rho} \right) }-Ch
1+\mu \,\bfgamma^T\bfdelta % \geqslant 1 + \mu \bar{\bfgamma}^T\bfdelta - 4\sqrt{2}(1+(1-2\epsilon^2)^{-3/2})\kappa \mu h
 % \geqslant 1- \frac{4\sqrt{\bar{\kappa}}}{1-\bar{\kappa}} \left( \frac{1}{\rho}-1 \right).
  \geqslant 1- 4\mu\sqrt{\bar{\kappa}}.
\end{equation}
And for the linear and Crouzeix-Raviart $\psi_{i, T}, i\in\mathcal{I}$, if $F$ is chosen such that $\mathbf{n}(F)=\bar{\mathbf{ n}}$, then
\begin{equation}
\label{eq_lower_bound_2}
1+\mu\bfgamma^T\bfdelta=1+\mu\bar{\bfgamma}^T\bfdelta\geqslant1.
\end{equation}
In addition, for the bilinear $\psi_{i, T}, i\in\mathcal{I}$, if $F$ is chosen such that $F_{\bot}$ is the midpoint of $\overline{DE}$, then
\begin{equation}
\label{eq_lower_bound_3}
1+\mu\bfgamma^T\bfdelta\geqslant1+\mu\bar{\bfgamma}^T\bfdelta(1-4\sqrt{2}\sqrt{\bar{\kappa}}).
\end{equation}
\end{lemma}
\begin{proof}
By Lemma \ref{lem2_3}, \eqref{eqn_psi_iT_values_bounds}, and Remark \ref{rem2_1}, we have
%shows that $\|\mathbf{ n}(F)-\mathbf{ n}(l)\|<Ch$ and $\|F-F_{\bot}\|<Ch^2$. So we have $\| \nabla \psi_{i,T}(F)-\nabla \psi_{i,T}(F_{\bot})\|<C$, for $\forall i\in\mathcal{I}$. Then, we have
\begin{equation*}
\begin{split}
\label{}
& |\nabla\psi^P_{i,T}(F)\cdot\mathbf{ n}(F) -  \nabla\psi^P_{i,T}(F_{\bot})\cdot \bar{\bfn}| \\
\leqslant& |(\nabla\psi^P_{i,T}(F) -  \nabla\psi^P_{i,T}(F_{\bot}))\cdot\mathbf{ n}(F)|+|\nabla\psi^P_{i,T}(F_{\bot})\cdot(\mathbf{ n}(F)- \bar{\bfn})| \\
 \leqslant& \|(\nabla\psi^P_{i,T}(F) -  \nabla\psi^P_{i,T}(F_{\bot}))\|\cdot \|\mathbf{ n}(F)\|+ \|\nabla\psi^P_{i,T}(F_{\bot})\|\cdot \|(\mathbf{ n}(F)- \bar{\bfn})\| \\
\leqslant&  4(1+(1-2\epsilon^2)^{-3/2})\kappa, ~i \in \mathcal{I}
\end{split}
\end{equation*}
which implies $\norm{\bfgamma - \bar{\bfgamma}} \leq 2(1+2(1-2\epsilon^2)^{-3/2})\kappa$. For all the types of finite elements considered in this article, by $\|L\|_{\infty,T} \leq \sqrt{2}h$, we have $\norm{\bfdelta} \leq \sqrt{2}h$ and therefore,
\begin{equation*}
\begin{split}
\label{}
 \mu \bfgamma^T\bfdelta = \mu \bar{\bfgamma}^T\bfdelta + \mu(\bfgamma^T - \bar{\bfgamma}^T)\bfdelta \geq \mu \bar{\bfgamma}^T\bfdelta- 4\sqrt{2}(1+(1-2\epsilon^2)^{-3/2})\kappa \mu h,
\end{split}
\end{equation*}
which yields \eqref{eqn4_22} by the assumption \eqref{h_small_assump}. Furthermore, for linear finite elements, if $F$ is chosen such that $\mathbf{ n}(F)=\bar{\mathbf{ n}}$, then $\bar{\bfgamma}=\bfgamma$; thus Lemma \ref{lemma4_3} and \eqref{bound_for_k} imply $1+\mu \,\bar{\bfgamma}^T\bfdelta\geqslant1$. For the bilinear finite elements, if $F$ is chosen such that $F_{\bot}$ is  the mid point of $\overline{DE}$, then by \eqref{eq_Ldelgam_2} we have
\begin{eqnarray*}
 \mu \bfgamma^T\bfdelta = \mu \bar{\bfgamma}^T\bfdelta + \mu(\bfgamma^T - \bar{\bfgamma}^T)\bfdelta
&\geq& \mu \bar{\bfgamma}^T\bfdelta - \mu 4(1+(1-2\epsilon^2)^{-3/2})\kappa (2h\bar{\bfgamma}^T\bfdelta) \\
&\geq& \mu\bar{\bfgamma}^T\bfdelta ( 1 - 8(1+(1-2\epsilon^2)^{-3/2})\kappa h ),
% \geqslant  \mu\bar{\bfgamma}^T\bfdelta(1-4\sqrt{2}\sqrt{\bar{\kappa}})
\end{eqnarray*}
which leads to \eqref{eq_lower_bound_3} by the assumption \eqref{h_small_assump}.
%Hence, $1+ \mu \bfgamma^T \bfdelta \geq 1 + \mu \bar{\bfgamma}^T\bfdelta-C\kappa \mu h$ which leads to \eqref{eqn4_22} by Lemma \ref{lemma4_3} and $\mu \geq 0$.
\end{proof}


\begin{thm}[$\mathbf{ Unisolvence}$]
\label{thm4_1}
Let $\mathcal{T}_h$ be a mesh satisfying \eqref{h_small_assump} with $\bar{\kappa}$ specified therein for linear and Crouzeix-Raviart $\psi_{i, T}, i\in\mathcal{I}$ and
\begin{equation}
\sqrt{\bar{\kappa}} < \frac{1}{4\sqrt{2}}, \text{~~for bilinear~~} \psi_{i, T}, i\in\mathcal{I}, \label{rem_unisol_eq2_1}
\end{equation}
and for some $\lambda \in (0, 1)$
\begin{equation}
\sqrt{\bar{\kappa}} \leq \frac{\rho(1- \lambda)}{4 - (3 + \lambda)\rho}, \text{~~for rotated-$Q_1$~~} \psi_{i, T}, i\in\mathcal{I}. \label{rem_unisol_eq1_1}
\end{equation}
In addition, we assume $F$ and $F_{\bot}$ are chosen such that estimates given by Lemma \ref{lemma4_4} hold. Then, given any vector $v=(v_1,v_2,v_3,v_4)\in \mathbb{R}^4$ for the bilinear and rotated $Q_1$ case (or $v=(v_1,v_2,v_3)\in \mathbb{R}^3$ for the linear and Crouzeix-Raviart case), there exists one and only one IFE function $\phi^P_T$ in the form of \eqref{eq4_10} satisfying \eqref{eq4_8}-\eqref{eq4_12}.
\end{thm}
\begin{proof}
For linear and Crouzeix-Raviart $\psi_{i, T}, i\in\mathcal{I}$, by \eqref{eq_lower_bound_2} in Lemma \ref{lemma4_4}, we have $1+\mu \bfgamma^T\bfdelta \neq 0$.
For bilinear $\psi_{i, T}, i\in\mathcal{I}$, by \eqref{rem_unisol_eq2_1}, we have
$
1 - 4\sqrt{2}\sqrt{\bar{\kappa}} > 0
$
which leads to $1+\mu \bfgamma^T\bfdelta \neq 0$ because of \eqref{eq_lower_bound_3} in Lemma \ref{lemma4_4}. Similarly,
for rotated-$Q_1$ $\psi_{i, T}, i\in\mathcal{I}$, by \eqref{rem_unisol_eq1_1}, we have
\begin{align}
1 - 4\mu \sqrt{\bar{\kappa}} \geq 1-\frac{4\sqrt{\bar{\kappa}}}{1-\bar{\kappa}}\left(\frac{1}{\rho}-1\right) \geq \lambda \label{rem_unisol_eq3}
\end{align}
which, by \eqref{eqn4_22} in Lemma \ref{lemma4_4}, leads to $1+\mu \bfgamma^T\bfdelta \neq 0$ again.
Hence, by the well known results about the \textit{Sherman-Morrison} formula, the matrix in the linear system \eqref{eqn4_15} is nonsingular which together with \eqref{eqn4_18} lead to the existence and uniqueness for coefficients $c_i, i \in \mathcal{I}^-$ and $c_0$ of $\phi_T^P(X)$.
\end{proof}

\begin{rem}
\label{rem_unisol}
Theorem \ref{thm4_1} provides guidelines on the choice for the angle allowance parameter $\bar{\kappa}$ needed in \eqref{h_small_assump} for
bilinear and rotated-$Q_1$ $\psi_{i, T}, i\in\mathcal{I}$. In the bilinear case, condition \eqref{rem_unisol_eq2_1} suggests an upper bound for $\bar{\kappa}$
%\begin{eqnarray*}
%\sqrt{\bar{\kappa}} < \frac{1}{4\sqrt{2}}
%\end{eqnarray*}
which is nevertheless independent of $\rho$. In the rotated-$Q_1$ case, condition \eqref{rem_unisol_eq1_1} leads to the following upper bound for $\bar{\kappa}$:
$\sqrt{\bar{\kappa}}<\frac{\rho}{4- 3\rho}$
%\begin{eqnarray*}
%\sqrt{\bar{\kappa}}<\frac{\rho}{4- 3\rho}
%\end{eqnarray*}
which depends on $\rho$, and this restriction on $\bar{\kappa}$ becomes severer when $\rho$ approaches $0$.
\end{rem}

\begin{rem}
\label{rem4_1}
Under the conditions given in Theorem \ref{thm4_1}, we can apply the \textit{Sherman-Morrison} formula to express the solution to \eqref{eqn4_15} explicitly as
\begin{equation}
\label{eqn4_23}
\bfc = \bfb - \mu \frac{(\bfgamma^T \bfb)\bfdelta}{1 + \mu \bfgamma^T\bfdelta}
\end{equation}
which facilitates both analysis and computation for these IFE spaces.
\end{rem}

%unisolvency


% define basis functions
On each interface element $T$, Theorem \ref{thm4_1} guarantees the existence and uniqueness of the IFE shape functions $\phi^P_{i,T}$, $i\in \mathcal{I}$
satisfying \eqref{eq4_10}-\eqref{eq4_9} such that
\begin{equation}
\label{eq4_21}
\phi^P_{i,T}(M_j)=\delta_{ij},~~i, j \in \mathcal{I},
\end{equation}
where $\delta_{ij}$ is the Kronecker delta function. Therefore, we can define the local IFE space on each interface element $T\in\mathcal{T}^i_h$ as
\begin{equation}
\label{eq4_22}
S^{P}_h(T)=\textrm{Span}\{ \phi^P_{i,T}:\;i\in \mathcal{I} \}.
\end{equation}

%Bound of basis functions
%We emphasize that the construction of basis functions and the argument for unisolvency is similar for other types of interface elements. Moreover we have the following properties of basis functions.


%\begin{equation}
%\label{eq4_6}
%\phi^{-}_T\in \Pi_{T}, \;\;\; and \;\;\;  \phi^{+}_T\in \Pi_{T},
%\end{equation}
%\begin{equation}
%\label{eq4_7}
%d(\phi^{-}_T)=d(\phi^{+}_T),
%\end{equation}
%\begin{equation}
%\label{eq4_8}
%\phi^{-}_T|_{l}=\phi^{+}_T|_{l},
%\end{equation}
%\begin{equation}
%\label{eq4_9}
%\beta^-\nabla\phi^{-}_T(F)\cdot\mathbf{ n}(F)=\beta^+\nabla\phi^{+}_T(F)\cdot\mathbf{ n}(F).
%\end{equation}

\commentout{
As in \cite{2013ZhangTHESIS}, depending on whether the interface cuts two adjacent or opposite edges, there are two types of interface elements.
According the location of $M_i, i \in \mathcal{I}$, a Type I element can be in $3$ cases and a Type II element can be in $2$ cases, see Figure \ref{fig_noncon_P_TypeI} and Figure \ref{fig_noncon_P_TypeII} for
illustrations of typical configurations.



 Let $F$ be an arbitrary point on $\Gamma\cap T$. We emphasize that, once chosen, $F$ should be fixed for each interface element.
%graph
\begin{figure}[H]
\centering
\includegraphics[width=6in]{figure5_1-eps-converted-to.pdf}
\caption{ \textbf{Type I} interface element and 3 cases }
\label{fig_noncon_P_TypeI}
\end{figure}

%graph
\begin{figure}[H]
\centering
\includegraphics[width=4in]{figure5_2-eps-converted-to.pdf}
\caption{ \textbf{Type II} interface element and 2 cases }
\label{fig_noncon_P_TypeII}
\end{figure}
}

%To be more specific, we assume $T$ is of \textbf{Type I} to exemplify the construction of IFE basis functions.

%\left\{
%\begin{aligned}
% \phi^{P\;-}_T(X) & = \phi^{P\;+}_T(X)+c_0L(X) & \text{if} \;\; X\in T^- \\
% \phi^{P\;+}_T(X)&  = \sum_{i\in\mathcal{I}^-}c_i\psi^P_{i,T}(X)+\sum_{i\in\mathcal{I}^+}v_i\psi^P_{i,T}(X)& \text{if} \;\; X\in T^+
%\end{aligned}
%\right..








\subsection{Properties of the IFE Shape Functions}

In this section, we present some fundamental properties of IFE shape functions.
%Unless otherwise clarified, we generally assume that the mesh size
%$h$ is small enough such that all the estimates established in the previous sections hold.

\begin{thm}[Bounds of IFE shape functions]
\label{thm4_2}
Under the conditions given in Theorem \ref{thm4_1}, we have the following estimates:
\begin{itemize}
\item
For rotated-$Q_1$ and Crouzeix-Raviart $\phi^P_{i,T}$, $i\in \mathcal{I},~\forall \,T \in \mathcal{T}^i_h$,
\begin{equation}
\label{eq4_24}
|\phi^P_{i,T}|_{k,\infty,T}\leqslant  \frac{C}{\rho}h^{-k}, k= 0, 1, %~\forall \,T \in \mathcal{T}^i_h,
\end{equation}
where $C$ depends also on $\lambda$ for rotated-$Q_1$ case.
\item
For linear $\phi^P_{i,T}$, $i\in \mathcal{I}, ~\forall \,T \in \mathcal{T}^i_h$,
%\begin{subequations}
%\label{linear_basis_bounds}
%\begin{equation}
%\label{linear_basis_bounds_1}
%|\phi^P_{i,T}|_{0,\infty,T}\leqslant  \frac{C}{\sqrt{\rho}}, ~~|\phi^P_{i,T}|_{0,\infty,T^-}\leqslant  \frac{C}{\rho}.
%\end{equation}
%\begin{equation}
%\label{linear_basis_bounds_2}
%|\phi^P_{i,T}|_{1,\infty,T^+}\leqslant  \frac{C}{\sqrt{\rho}}h^{-1}, ~~~ |\phi^P_{i,T}|_{1,\infty,T^-}\leqslant  \frac{C}{\rho}h^{-1}.
%\end{equation}
%\end{subequations}
\begin{equation}
\label{linear_basis_bounds}
|\phi^P_{i,T}|_{k,\infty,T^+}\leqslant  \frac{C}{\sqrt{\rho}}h^{-k}, ~~|\phi^P_{i,T}|_{k,\infty,T^-}\leqslant  \frac{C}{\rho}h^{-k}, ~~ k=0,1.
\end{equation}
\item
For bilinear $\phi^P_{i,T}$, $i\in \mathcal{I}, ~\forall \,T \in \mathcal{T}^i_h$,
%\begin{subequations}
%\begin{equation}
%\label{bilinear_basis_bounds_1}
%|\phi^P_{i,T}|_{0,\infty,T}\leqslant  C,
%\end{equation}
%\begin{equation}
%\label{bilinear_basis_bounds_2}
%|\phi^P_{i,T}|_{1,\infty,T^+}\leqslant  Ch^{-1}, ~~~ |\phi^P_{i,T}|_{1,\infty,T^-}\leqslant  \frac{C}{\rho}h^{-1}
%\end{equation}
%\end{subequations}
\begin{equation}
\label{bilinear_basis_bounds_2}
|\phi^P_{i,T}|_{k,\infty,T^+}\leqslant  Ch^{-k}, ~~~ |\phi^P_{i,T}|_{k,\infty,T^-}\leqslant  \frac{C}{\rho}h^{-k}, ~~ k=0,1.
\end{equation}
\end{itemize}
\end{thm}
\begin{proof}
%Since $|\psi^P_{i,T}|_{k,\infty,T}\leqslant Ch^{-k}$, $i\in \mathcal{I}$, we only need to prove the coefficients in \eqref{eq4_13} are uniformly bounded.
For convenience, we let $\mathbf{ e}=(e_{j})_{j\in\mathcal{I}}$ be the unit vector constructing the basis functions $\phi^P_{i,T}$, which could be $(1,0\cdots,0)$, $(0,1\cdots,0)$, $\cdots$, $(0,0\cdots,1)$, and $\mathbf{ e}^-=(e_{j})_{j\in\mathcal{I}^-}$. Let $w=\sum_{j\in\mathcal{I}^+} \nabla\psi^P_{j,T}(F)\cdot\mathbf{ n}(F)e_j$. Then \eqref{eqn4_17} implies $\bfb=\mathbf{ e}^- - \mu w \bfdelta$ and plugging it into the \textit{Sherman}-\textit{Morrison} formula \eqref{eqn4_23} leads to
\begin{equation}
\label{bounded_1}
\begin{split}
\bfc %=& \frac{\bfb+k(\bfgamma^T\bfdelta\bfb - (\bfgamma^T\bfb)\bfdelta )}{1+k\bfgamma^T\bfdelta} \\
 %= &\frac{ \mathbf{e}^--k\omega\bfdelta+k(\bfgamma^T\bfdelta\mathbf{ e}^- -\bfgamma^T\mathbf{e}^-\bfdelta ) }{ 1+k\bfgamma^T\bfdelta }
 = \mathbf{e}^- - \mu \frac{(w + \bfgamma^T\mathbf{ e}^-)\bfdelta }{1 + \mu \bfgamma^T\bfdelta },
\end{split}
\end{equation}
and plugging \eqref{bounded_1} into \eqref{eqn4_18} yields
\begin{equation}
\label{bounded_2}
c_0=\frac{\mu (w + \bfgamma^T\mathbf{ e}^-)}{1 + \mu \bfgamma^T\bfdelta}.
\end{equation}
Since
$
\|\nabla\psi^P_{i,T}\|_{\infty,T}\leqslant Ch^{-1}, ~i\in \mathcal{I}, ~\|L\|_{\infty,T}<Ch, ~|\nabla L\|_{\infty,T}<C \label{eq_thm4_2_1}
$
for some constants $C$ independent of the location of the interface and $\beta^{\pm}$, we have $\norm{\bfgamma} \leq Ch^{-1}, ~\norm{\bfdelta} \leq Ch$, $\norm{\bfb} \leq C$ and $|w|\leqslant Ch^{-1}$. \\

When $\phi^P_{i,T}$, $i\in \mathcal{I}$ are rotated-$Q_1$ polynomials, we can apply \eqref{bound_for_k} and \eqref{rem_unisol_eq1_1} to \eqref{bounded_1} to obtain:
\begin{equation}
\begin{split}
\label{eqn3_13}
\|\bfc\| \leq \|\mathbf{ e}^-\|+ \frac{1}{1- \bar{\kappa}}\left( \frac{1}{\rho}-1 \right) \frac{ (|w|+\| \bfgamma \| \| \mathbf{ e}^-\|)\|\bfdelta\| }{1-\frac{4\sqrt{\bar{\kappa}}}{1-\bar{\kappa}}\left( \frac{1}{\rho}-1 \right)} \leq \frac{C}{\rho}.
\end{split}
\end{equation}
% since if $\rho\leqslant1$, then $\|c\|<C+C(\frac{1}{\rho}-1)=C\frac{1}{\rho}$ and if $\rho>1$, then $\|c\|<C+C(\rho-1)=C\rho$.
Applying similar arguments to \eqref{bounded_2} we have $\|c_0\|<C/\rho$. Constant $C$ in these inequalities depends on $\lambda$. Finally, \eqref{eq4_24} follows from applying these bounds for $\bfc$ and $\bfc_0$ and the bounds for standard finite element basis functions $\psi_{i,T}$ to the formula of $\phi^P_{i,T}$ given in \eqref{eq4_13}. When $\phi^P_{i,T}$, $i\in \mathcal{I}$ are
Crouzeix-Raviart polynomials, we apply \eqref{bound_for_k} and \eqref{eq_lower_bound_2} to \eqref{bounded_1} and \eqref{bounded_2} to obtain:
\begin{eqnarray}
\|\bfc\| &\leq&  \|\mathbf{ e}^-\|+ \frac{1}{1- \bar{\kappa}}\left( \frac{1}{\rho}-1 \right) (|w|+\| \bfgamma \| \| \mathbf{ e}^-\|)\|\bfdelta\|  \leq \frac{C}{\rho}, \nonumber \\
| c_0 | &\leq& \frac{1}{1- \bar{\kappa}}\left( \frac{1}{\rho}-1 \right) (|w|+\| \bfgamma \| \| \mathbf{ e}^-\|)\|\bfdelta\| \leq  \frac{C}{\rho}. \label{bounded_7}
\end{eqnarray}
Then \eqref{eq4_24} in this case follows from the same arguments used above for the rotated-$Q_1$ case.

When $\phi^P_{i,T}$, $i\in \mathcal{I}$ are linear polynomials, by using \eqref{eq_Ldelgam_1} and Lemma \ref{lemma4_4} in \eqref{bounded_1}, we have
\begin{equation}
\label{bounded_3}
\| \bfc \| \leq C+C\frac{\mu\sqrt{\bar{\bfgamma}^T\bfdelta}}{1+\mu\bar{\bfgamma}^T\bfdelta} \leq C+C\max{ (1, \sqrt{\mu}) } \leq \frac{C}{\sqrt{\rho}}.
\end{equation}
%Also, since $L(X)$ is a linear function, by \eqref{eq_Ldelgam_1} and \eqref{bounded_2} we have
%\begin{equation}
%\label{bounded_4}
%\abs{c_0L(X)} \leq \max_{i \in \mathcal{I}^-}\abs{c_0L(M_i)} \leq C\frac{\mu h^{-1}|L(M_i)|}{1+\mu\bar{\bfgamma}^T\bfdelta} \leq
%C\frac{\mu \sqrt{\bar{\bfgamma}^T\bfdelta}}{1+\mu\bar{\bfgamma}^T\bfdelta} \leq \frac{C}{\sqrt{\rho}},~~\forall X \in T^-.
%\end{equation}
Also, estimate in \eqref{bounded_7} is valid. Then, estimates in \eqref{linear_basis_bounds} follow from applying \eqref{bounded_3} and \eqref{bounded_7} to \eqref{eq4_13}.
%Then, applying \eqref{eqn_psi_iT_values_bounds} and \eqref{bounded_3} to \eqref{eq4_13} yields the inequality on $T^+$ in \eqref{linear_basis_bounds_2}.
%
%And from \eqref{eq4_13}, on $T^-$, we only need to estimate $c_0L(X)$ which reduces to the estimation of $c_0L(M_i)$, $i\in\mathcal{I}^-$, since $L(X)$ is a linear function and vanishes on $l$. And by \eqref{bounded_2} and similar arguments above, we obtain
%
%Besides, \eqref{bounded_3} also shows \eqref{linear_basis_bounds_2} on $T^+$. For estimates under $|\cdot|_{1,\infty,T^-}$, using \eqref{bounded_2} and by similar arguments as \eqref{eqn3_13}, we could show that $\|c_0\nabla L\|<C/\rho$. \\

Finally, when $\phi^P_{i,T}$, $i\in \mathcal{I}$ are bilinear polynomials, we can apply \eqref{eq_Ldelgam_2} and Lemma \ref{lemma4_4} in \eqref{bounded_1} to have
\begin{equation}
\label{bounded_5}
\| \bfc \|<C+C\frac{\mu\bar{\bfgamma}^T\bfdelta}{1 + (1-4\sqrt{2}\sqrt{\bar{\kappa}}) \mu\bar{\bfgamma}^T\bfdelta}<C,
\end{equation}
because $1-4\sqrt{2}\sqrt{\bar{\kappa}}$ is a positive constant by the condition \eqref{rem_unisol_eq2_1}.  Also, the estimate for $c_0$ is similar to \eqref{bounded_7}. Then, \eqref{bilinear_basis_bounds_2} follows from applying these bounds for $\bfc$ and $c_0$ to \eqref{eq4_13}.
%using \eqref{eq_Ldelgam_2} in \eqref{bounded_2} we have
%\begin{equation}
%\label{bounded_6}
%\abs{c_0L(X)} \leq \max_{i \in \mathcal{I}}\abs{c_0L(M_i)} \leq C\frac{\mu h^{-1}|L(M_i)|}{1+\mu\bar{\bfgamma}^T\bfdelta} \leq
%C\frac{\mu \bar{\bfgamma}^T\bfdelta}{1+\mu\bar{\bfgamma}^T\bfdelta} \leq C,~~\forall X \in T.
%\end{equation}
\end{proof}


%partition of utility
\begin{lemma} [Partition of Unity]
\label{lem:POU}
On every interface element $T \in \mathcal{T}_h^i$, we have
\begin{equation}
\label{eq4_25}
\sum_{i\in \mathcal{I}}\phi^P_{i,T}(X) = 1,~~\forall X \in T,
\end{equation}
\begin{equation}
\label{eq4_26}
\nabla \left( \sum_{i\in\mathcal{I}}\phi^P_{i,T}(X) \right)=\sum_{i\in\mathcal{I}}\nabla\phi^P_{i,T}(X)=0, ~~\forall X \in T.
\end{equation}
\end{lemma}
\begin{proof}
Let $p(X) = \sum_{i\in \mathcal{I}}\phi^P_{i,T}(X)$ and $q(X) = 1$. Obviously both $p(X)$ and $q(X)$ are in the format of \eqref{eq4_10} and they satisfy \eqref{eq4_8} and
\eqref{eq4_9}. Furthermore, it is easy to verify that $p(M_i) = 1 = q(M_i), i \in \mathcal{I}$. Hence,
\begin{eqnarray*}
\sum_{i\in \mathcal{I}}\phi^P_{i,T}(X) = p(X) = q(X) = 1
\end{eqnarray*}
according to the unisolvence stated in Theorem \ref{thm4_1}. Property \eqref{eq4_26} follows from \eqref{eq4_25} directly.
\end{proof}

Now, on every $T \in \mathcal{T}_h^i$, choosing arbitrary points $\overline{X}_i \in l, i \in \mathcal{I}$, we can
construct two vector functions as follows:
\begin{subequations}
\label{eq4_27}
\begin{equation}\label{eq4_27a}
\bfLambda_1(X)=\sum_{i\in\mathcal{I}}(M_i-X)\phi^{P\;+}_{i,T}(X)+\sum_{i\in \mathcal{I}^-}(\overline{M}^+(F) - I)^T(M_i-\overline{X}_i)\phi^{P\;+}_{i,T}(X), \;\;\; \text{if}\; X\in T^+,
\end{equation}
\begin{equation}\label{eq4_27b}
\bfLambda_2(X)=\sum_{i\in\mathcal{I}}(M_i-X)\phi^{P\;-}_{i,T}(X)+\sum_{i\in \mathcal{I}^+}(\overline{M}^-(F) - I)^T(M_i-\overline{X}_i)\phi^{P\;-}_{i,T}(X), \;\;\; \text{if}\; X\in T^-.
\end{equation}
\end{subequations}
It follows from the Lemma \ref{lem2_5} that the two functions in \eqref{eq4_27} are independent of the location of $\overline{X}_i\in l, i \in \mathcal{I}$. Furthermore, from the partition of unity stated in Lemma \ref{lem:POU}, we have
\begin{subequations}
\label{eq4_27_modified}
\begin{align}
\label{}
   \bfLambda_1(X)&=\sum_{i\in \mathcal{I}}M_i\phi^{P\;+}_{i,T}(X) - X + \sum_{i \in \mathcal{I}^-}(\overline{M}^+(F) -I)^T(M_i-\overline{X}_i)\phi^{P\;+}_{i,T}(X),    \\
   \bfLambda_2(X)&=\sum_{i\in \mathcal{I}}M_i\phi^{P\;-}_{i,T}(X) - X + \sum_{i \in \mathcal{I}^+}(\overline{M}^-(F) -I)^T(M_i-\overline{X}_i)\phi^{P\;-}_{i,T}(X),
\end{align}
\end{subequations}
which imply that each component of $\bfLambda_1(X)$ and $\bfLambda_2(X)$ is a polynomial in $\Pi_T$ because $\phi^{P\;s}_{i,T}(X)\in \Pi_T, s = \pm$, for $i\in\mathcal{I}$.
We consider two auxiliary vector functions
\begin{equation}
\begin{split}
\label{eq4_29}
&\bfLambda^+(X)=\bfLambda_1(X), \;\; \text{and} \;\;\;\; \bfLambda^-(X)=(\overline{M}^+(F))^T\bfLambda_2(X).
\end{split}
\end{equation}
Let $d(\bfLambda^s), s = \pm$ be the vector of the coefficients of the second degree term in each component of $\bfLambda^s$.
\begin{lemma}
\label{lem4_2}
$\bfLambda^+$ and $\bfLambda^-$ are such that $d(\bfLambda^+)=d(\bfLambda^-)$.
\end{lemma}
\begin{proof}
Let $d(\phi^{P\; +}_{i,T})=d(\phi^{P\; -}_{i,T})=d_i$, $i\in \mathcal{I}$. By the partition of unity, $\sum_{i\in\mathcal{I}}d_i=0$.
By using \eqref{eq2_10} given in Lemma \ref{lem2_4} and using Lemma \ref{lem2_5} to interchange $\overline{X}_i, i \in \mathcal{I}$ with a fixed $\overline{X} \in l$ and, we have
%By Lemma \ref{lem2_5}, we can take $\overline{X}_i$ be the same point $\overline{X}$ on $l$.
\begin{equation*}
\begin{split}
\label{}
d(\bfLambda^-)&=\sum_{i\in\mathcal{I}}(\overline{M}^+(F))^TM_id_i+\sum_{i\in \mathcal{I}^+}(I-\overline{M}^+(F))^T(M_i-\overline{X})d_i\\
=&\sum_{i\in\mathcal{I}^-}(\overline{M}^+(F))^TM_id_i+\sum_{i\in\mathcal{I}^+}M_id_i - (I-\overline{M}^+(F))^T\overline{X}\sum_{i\in \mathcal{I}^+}d_i\\
=&\sum_{i\in\mathcal{I}}M_id_i+\sum_{i\in\mathcal{I}^-}(\overline{M}^+(F) - I)^TM_i d_i
+ (I-\overline{M}^+(F))^T\overline{X}\sum_{i\in \mathcal{I}^-}d_i \\
=& \sum_{i\in\mathcal{I}}M_id_i+\sum_{i\in\mathcal{I}^-}(\overline{M}^+(F) - I)^T (M_i - \overline{X}_i) d_i
\end{split}
\end{equation*}
which is exactly $d(\bfLambda^+)$.
\end{proof}

\begin{lemma}
\label{lem4_3}
$\bfLambda^+$ and $\bfLambda^-$ satisfy the condition \eqref{eq4_8}, i.e., $\bfLambda^+(X)|_l=\bfLambda^-(X)|_l$.
\end{lemma}
\begin{proof}
Since $\bfLambda^s(X), s = \pm$ are independent of points $\overline{X}_i \in l, i \in \mathcal{I}$, we can replace these points by an arbitrary $\overline{X} \in l$ without
changing values of $\bfLambda^s(X), s = \pm$. Then, since $\phi^P_{i,T}$ satisfies \eqref{eq4_8}, applying \eqref{eq2_10}, we have
\begin{equation*}
\begin{split}
%\label{ }
\bfLambda^+(\overline{X})=&\sum_{i\in\mathcal{I}}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})+\sum_{i\in\mathcal{I}^-}(\overline{M}^+(F)-I)^T(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})\\
=&\sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X}) + (\overline{M}^+(F))^T\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})\\
&= (\overline{M}^+(F))^T \left((\overline{M}^-(F))^T \sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})+\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X}) \right)\\
&=(\overline{M}^+(F))^T \left(\sum_{i\in\mathcal{I}}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})
+ (\overline{M}^-(F) - I)^T \sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})\right)\\
\end{split}
\end{equation*}
%Again, from \eqref{eq2_10}, one can verify that
%\begin{equation*}
%\begin{split}
%%\label{ }
%\bfLambda^+(\overline{X})&= (\overline{M}^+)^T \left((\overline{M}^-)^T \sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})+\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X}) \right)\\
%&=(\overline{M}^+)^T \left((\overline{M}^--I)^T \sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X})+\sum_{i\in\mathcal{I}}(M_i-\overline{X})\phi^{P\;-}_{i,T}(\overline{X}) \right)\\
%\end{split}
%\end{equation*}
which is exactly $\bfLambda^-(\overline{X})$.
\end{proof}

\begin{lemma}
\label{lem4_4}
$\bfLambda^+$ and $\bfLambda^-$ satisfy the condition \eqref{eq4_9}, i.e., $\beta^+\nabla\bfLambda^+(F)\cdot\mathbf{ n}(F)=\beta^-\nabla\bfLambda^-(F)\cdot\mathbf{ n}(F)$, where the gradient operator is understood as the gradient on each component.
\end{lemma}
\begin{proof}
Again, Lemma \ref{lem2_5} allows us to exchange $\overline{X}_i, i \in \mathcal{I}$ for an arbitrary point $\overline{X} \in l$ in the discussion below. By \eqref{eq4_9}, \eqref{eq4_26}, \eqref{eq2_10}, and \eqref{eq2_12} we have
\begin{equation*}
\begin{split}
\label{}
& \beta^+ \nabla \bfLambda^+(F)\cdot\mathbf{ n}(F) = \sum_{i\in\mathcal{I}}M_i\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot\mathbf{ n}(F) \\
& \hspace{0.5in} + (\overline{M}^+(F) - I)^T\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot\mathbf{ n}(F)
- \beta^+\mathbf{ n}(F) \\
=&\sum_{i\in\mathcal{I}^+}M_i\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot\mathbf{ n}(F) + \sum_{i\in\mathcal{I}^-}\overline{X}\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot\mathbf{ n}(F)  \\
& \hspace{0.5in}+ (\overline{M}^+(F))^T\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot\mathbf{ n}(F)
- \beta^+\mathbf{ n}(F)\\
=&\sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot\mathbf{ n}(F)
+(\overline{M}^+(F))^T\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot\mathbf{ n}(F) - \beta^+\mathbf{ n}(F)\\
=& (\overline{M}^+(F))^T \left( (\overline{M}^-(F))^T\sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot \bfn(F) \right. \\
&\hspace{1in} \left. +\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot \bfn(F) - \beta^+(\overline{M}^-(F))^T\cdot \bfn(F) \right)\\
=& (\overline{M}^+(F))^T \left( (\overline{M}^-(F))^T\sum_{i\in\mathcal{I}^+}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot \bfn(F) \right. \\
&\hspace{1in} \left. +\sum_{i\in\mathcal{I}^-}(M_i-\overline{X})\beta^-\nabla\phi^{P\;-}_{i,T}(F)\cdot \bfn(F) - \beta^-\cdot \bfn(F) \right)
= \beta^-\nabla\bfLambda^-(F)\cdot\mathbf{ n}(F).
\end{split}
\end{equation*}

\end{proof}

%% fundamental identity
\begin{thm}
\label{thm4_3}
On every interface element $T \in \mathcal{T}_h^i$ we have
\begin{subequations}
\label{eq_thm4_3}
\begin{align}
&\sum_{i\in\mathcal{I}}(M_i-X)\phi^{P\;-}_{i,T}(X)+\sum_{i\in \mathcal{I}^+}(\overline{M}^-(F) - I)^T(M_i-\overline{X}_i)\phi^{P\;-}_{i,T}(X) = \bfLambda_1(X) = 0, ~\forall X\in T^-,
\label{eq_thm4_3_1} \\
&\sum_{i\in\mathcal{I}}(M_i-X)\phi^{P\;+}_{i,T}(X)+\sum_{i\in \mathcal{I}^-}(\overline{M}^+(F) - I)^T(M_i-\overline{X}_i)\phi^{P\;+}_{i,T}(X) = \bfLambda_2(X) = 0, ~\forall X\in T^+,
\label{eq_thm4_3_2}
\end{align}
\end{subequations}
and
\begin{subequations}
\label{eq_thm4_4}
\begin{align}
&\sum_{i\in\mathcal{I}}(M_i-X)\partial_d\phi^{P\;-}_{i,T}(X)+\sum_{i\in \mathcal{I}^+}\left[ (\overline{M}^-(F) - I)^T(M_i-\overline{X}_i)\partial_s\phi^{P\;-}_{i,T}(X) \right] - \mathbf{ e}_d = 0, \forall X \in T^-, \label{eq4_30} \\
&\sum_{i\in\mathcal{I}}(M_i-X)\partial_d\phi^{P\;+}_{i,T}(X)+\sum_{i\in \mathcal{I}^-}\left[ (\overline{M}^+(F) - I)^T(M_i-\overline{X}_i)\partial_s\phi^{P\;+}_{i,T}(X) \right] - \mathbf{ e}_d = 0, \forall X \in T^+, \label{eq4_31}
\end{align}
\end{subequations}
where $d = 1, 2$, $\partial_1 = \partial_x, \partial_2 = \partial_y$ are partial differential operators, and $\mathbf{ e}_d, d = 1, 2$ is the standard $d$-th unit vector in
$\mathbb{R}^2$.

%$\bfLambda_1(X)\equiv0$ and $\bfLambda_2(X)\equiv0$ for all $X \in T$.
\end{thm}
\begin{proof}
We define a piecewise vector polynomial on $T$ as
\begin{equation*}
%\label{ }
\bfLambda(X)=
\left\{\begin{array}{cc}
\bfLambda^+(X) & \text{if} \; X\in T^+, \\
\bfLambda^-(X) & \text{if} \; X\in T^-.
\end{array}\right.
\end{equation*}
First, the restriction of each component of $\bfLambda$ to $T^s, s = \pm$ is a polynomial in $\Pi_T$. By Lemmas \ref{lem4_2}-\ref{lem4_4}, the components of
$\bfLambda$ also satisfy \eqref{eq4_8} and \eqref{eq4_9}. In addition, we can easily see that $\bfLambda(M_i)= {\bf 0}, i\in \mathcal{I}$. Therefore,
by the unisolvence stated in Theorem \ref{thm4_1}, we have $\bfLambda^+(X)=\bfLambda_1(X)\equiv0$ and $\bfLambda^-(X)=(\overline{M}^+(F))^T\bfLambda_2(X)\equiv0$. Since $(\overline{M}^+(F))^T$ is nonsingular, we have $\bfLambda_2(X)\equiv0$. Therefore \eqref{eq_thm4_3} is proved. The proof for \eqref{eq_thm4_4} can be accomplished by differentiating \eqref{eq_thm4_3} and applying \eqref{eq4_27_modified}.
\end{proof}

%\begin{thm}
%\label{thm4_4}
%For every interface element $T \in \mathcal{T}_h^i$, the following identities hold:
%\begin{align}
%&\sum_{i\in\mathcal{I}}(M_i-X)\partial_s\phi^{P,-}_{i,T}(X)+\sum_{i\in \mathcal{I}^+}\left[ (\overline{M}^-(F) - I)^T(M_i-\overline{X}_i)\partial_s\phi^{P,-}_{i,T}(X) \right] - \mathbf{ e}_s = 0, \forall X \in T^-, \label{eq4_30} \\
%&\sum_{i\in\mathcal{I}}(M_i-X)\partial_s\phi^{P,+}_{i,T}(X)+\sum_{i\in \mathcal{I}^-}\left[ (\overline{M}^+(F) - I)^T(M_i-\overline{X}_i)\partial_s\phi^{P,+}_{i,T}(X) \right] - \mathbf{ e}_s = 0, \forall X \in T^+, \label{eq4_31}
%\end{align}
%%\begin{equation}
%%\label{eq4_30}
%%\sum_{i\in\mathcal{I}}(M_i-X)\partial_x\phi^{P,+}_{i,T}(X)+\sum_{i\in \mathcal{I}^-}\left[ (\overline{M}^+-I)^T(M_i-\overline{X}_i)\partial_x\phi^{P,+}_{i,T}(X) \right] - \mathbf{ e}_1 \equiv0
%%\end{equation}
%%\begin{equation}
%%\label{eq4_31}
%%\sum_{i\in\mathcal{I}}(M_i-X)\partial_y\phi^{P,+}_{i,T}(X)+\sum_{i\in \mathcal{I}^-}\left[ (\overline{M}^+-I)^T(M_i-\overline{X}_i)\partial_y\phi^{P,+}_{i,T}(X) \right] - \mathbf{ e}_2 \equiv0
%%\end{equation}
%where $s = x$ or $s = y$, $\mathbf{ e}_x=(1,0)^T$, and $\mathbf{ e}_y=(0,1)^T$. %And for $\phi^{P\;-}_{i,T}$, there are similar results.
%\end{thm}
%\begin{proof}
%The proof can be accomplished by a direct application of Theorem \ref{thm4_3} and \eqref{eq4_27_modified}.
%\end{proof}

%%%%%%%%%%%%%%%%%%%%%%%
%subsection
\subsection{Optimal Approximation Capabilities of IFE Spaces}

As usual, the local IFE spaces on elements can be employed to define the IFE function space globally on $\Omega$. As an example, we consider
\begin{equation}
\begin{split}
\label{eq4_23}
S_h(\Omega)=&\left\{ v\in L^2(\Omega):v|_{T}\in S^{P}_h(T); \right. \\
& \hspace{0.2in} \left. v|_{T_1}(M)=v|_{T_2}(M)~\forall M \in \mathcal{N}_h, \forall\,T_1, T_2 \in \mathcal{T}_h \text{~such that~} M \in T_1\cap T_2 \right\},
%S_h(\Omega)=&\left\{ v\in L^2(\Omega):v|_{T}\in S^{P}_h(T) \right. \\
%& \left. \text{if} \; T_1\cap T_2=M\in \mathcal{N}_h,\; \text{then}\; v|_{T_1}(M)=v|_{T_2}(M)  \right\}
\end{split}
\end{equation}
where $\mathcal{N}_h$ is the set of nodes in the mesh $\mathcal{T}_h$, and this implies every IFE
function in this space is continuous at every node in the mesh. However, as observed in \cite{2015LinLinZhang}, every $v\in S_h(\Omega)$ is usually discontinuous across the interface edges. When the interface is a generic curve, the discontinuity of $v\in S_h(\Omega)$ also occurs along the interface curve because these IFE shape functions are defined according to the actual interface. These features ensure the $H^1$ regularity of $v\in S_h(\Omega)$ in the subdomain
of $\Omega$ minus the union of interface elements.
%The continuity at the mesh nodes enables us to use a partial penalty Galerkin scheme for producing an
%IFE solution whose degree of freedom is not more than the number of nodes in the mesh.

\par
We proceed to show that these IFE spaces formed above by linear, bilinear, CR and rotated-$Q_1$ polynomials have the optimal approximation property from the point of view how
well the interpolation of a function $u$ in these IFE spaces can approximate $u$. First we define the interpolation operator on an element $T \in \mathcal{T}_h$ as the mapping $I_{h,T}: C^0(T) \rightarrow S_h^P(T)$ such that
\begin{equation}
\label{eq4_33}
I^P_{h,T}u(X)= \begin{cases}
\sum_{i\in\mathcal{I}}u(M_i)\psi^P_{i,T}(X), & \text{if~} T \in \mathcal{T}_h^n, \\
\sum_{i\in\mathcal{I}}u(M_i)\phi^P_{i,T}(X), & \text{if~} T \in \mathcal{T}_h^i.
\end{cases}
\end{equation}
Furthermore, the global IFE interpolation $I^P_{h}: C^0(\Omega)\rightarrow S_h(\Omega)$ can be defined piecewisely as
\begin{equation}
\label{eq4_34}
(I^P_hu)|_{T}=I^P_{h,T}u, \;\;\; \forall T\in \mathcal{T}_h.
\end{equation}

On every non-interface element $T\in\mathcal{T}^n_h$, the standard scaling argument \cite{1994BrennerScott,1978Ciarlet,1992RanacherTurek}
yields the following error estimate for the local interpolation $I_{h,T}^Pu$ on $T$:
\begin{equation}
\label{eq4_32}
\| I^P_{h,T}u-u\|_{0,T}+h| I^P_{h,T}u-u|_{1,T}\leqslant Ch^2|u|_{2,T}, ~\forall u \in H^2(T).
\end{equation}
However, how to use the scaling argument to derive an error bound for the interpolation on an interface element is unclear because the local IFE space $S^{P}_h(T)$ is interface dependent and it is not even a subspace of $H^1(T)$ in general.  Instead, we will use the multi-point Taylor expansion method \cite{2009HeTHESIS,2008HeLinLin,2004LiLinLinRogers,1982Xu,2013ZhangTHESIS}
to derive estimates for the IFE interpolation error.
%We will present detailed arguments for the error estimation on $T^-$ for each interface element $T\in \mathcal{T}_h^i$,
%and all the arguments apply to $T^+$ and lead to estimates on the whole $T$.

%expansion
\begin{thm}
\label{thm4_5}
Let $T \in \mathcal{T}_h^i$, assume $u\in PC_{int}^2(T)$. Then for any $\overline{X}_i\in l, i \in \mathcal{I}$, we have
%Assume $X\in T^-$ and , define $Y_{i}(t,X)=tM_i+(1-t)X$, $\forall i\in \mathcal{I}$. Note that for $i\in \mathcal{I}^+$, there exists $\tilde{t}_i$ such that $\widetilde{Y}_j=Y_j(\tilde{t}_j,X) \in \Gamma\cap T$. Let $\overline{X}_i\in l$, then
\begin{subequations}\label{eq4_35}
\begin{equation}
\label{eq4_35a}
I^P_{h,T}u(X)-u(X)=\sum_{i\in\mathcal{I}^{s'}}(E_i^s+F_i^s)\phi^P_{i,T}(X)+\sum_{i\in\mathcal{I}}R_i^s\phi^P_{i,T}(X),~\forall X \in T^s, ~s = \pm,
\end{equation}
\begin{equation}
\label{eq4_35b}
\partial_d( I^P_{h,T}u(X)-u(X))=\sum_{i\in\mathcal{I}^{s'}}(E_i^s + F_i^s)\partial_d\phi^P_{i,T}(X)+\sum_{i\in\mathcal{I}}R_i^s \partial_d\phi^P_{i,T}(X), ~\forall X \in T^s, ~s = \pm,
\end{equation}
\end{subequations}
%\begin{equation}
%\label{eq4_35}
%I^P_{h,T}u(X)-u(X)=\sum_{i\in\mathcal{I}^+}(E_i+F_i)\phi^P_{i,T}(X)+\sum_{i\in\mathcal{I}}R_i\phi^P_{i,T}(X)
%\end{equation}
%and
%\begin{equation}
%\label{eq4_36}
%\partial_s( I^P_{h,T}u(X)-u(X))=\sum_{i\in\mathcal{I}^+}(E_i+F_i)\partial_s\phi^P_{i,T}(X)+\sum_{i\in\mathcal{I}}R_i\partial_s\phi^P_{i,T}(X)
%\end{equation}
where $d=1$ or $2$, $R_i^s$ are given by \eqref{eq3_2} and \eqref{eq3_4}, and
\begin{equation}
\label{eq4_37}
\begin{split}
&E_i^s = \left( (M^s(\widetilde{Y}_i)-\overline{M}^s(F))\nabla u^s(X) \right)\cdot (M_i-\widetilde{Y}_i), \;\;\; s = \pm, ~i\in\mathcal{I}^{s'}, \\
&F_i^s = -\left( (\overline{M}^s(F) - I)\nabla u^s(X) \right) \cdot (\widetilde{Y}_i-\overline{X}_i), \;\;\; s = \pm, ~i\in\mathcal{I}^{s'}.
\end{split}
\end{equation}
\end{thm}
\begin{proof}
For $X \in T^s, s = \pm$, substituting the expansion \eqref{eq3_1} and \eqref{eq3_3} into the IFE interpolation \eqref{eq4_33}, we have
\begin{equation}
\begin{split}
\label{eq4_38}
I^P_{h,T}u(X) &= u^s(X)\sum_{i\in\mathcal{I}}\phi^P_{i,T}(X)+ \nabla u^s(X)\cdot \sum_{i\in\mathcal{I}}(M_i-X)\phi^P_{i,T}(X)\\
+ &\sum_{i\in \mathcal{I}^{s'}}\left( \left( M^s(\widetilde{Y}_i)-I \right)\nabla u^-(X) \right)\cdot (M_i-\widetilde{Y}_i)\phi^P_{i,T}(X) + \sum_{i\in\mathcal{I}}R_i^s\phi^P_{i,T}(X),
~X \in T^s, ~s = \pm.
\end{split}
\end{equation}
From Theorem \ref{thm4_3}, we have
\begin{equation}
\begin{split}
\label{eq4_39}
\sum_{i\in\mathcal{I}}(M_i-X)\phi^P_{i,T}(X)=-\sum_{i\in \mathcal{I}^{s'}}(\overline{M}^s(F) - I)^T(M_i-\overline{X}_i)\phi^{P}_{i,T}(X), ~\forall X \in T^s, ~s = \pm.
\end{split}
\end{equation}
Then, applying \eqref{eq4_39} and the partition of unity to \eqref{eq4_38} leads to
\begin{equation}
\begin{split}
\label{eq4_40}
I^P_{h,T}u(X) &= u^s(X)- \sum_{i\in\mathcal{I}^{s'}}\left( (\overline{M}^s(F) - I)\nabla u^s(X) \right) \cdot(M_i-\overline{X}_i)\phi^P_{i,T}(X)\\
+& \sum_{i\in \mathcal{I}^{s'}}\left( \left( M^s(\widetilde{Y}_i)-I \right)\nabla u^s(X) \right)\cdot (M_i-\widetilde{Y}_i)\phi^P_{i,T}(X)+\sum_{i\in\mathcal{I}}R_i^s\phi^P_{i,T}(X),
~X \in T^s, ~s = \pm
\end{split}
\end{equation}
from which \eqref{eq4_35a} follows by using $M_i-\overline{X}_i=(M_i- \widetilde{Y}_i)+(\widetilde{Y}_i -\overline{X}_i)$. For \eqref{eq4_35b}, applying the expansions \eqref{eq3_1} and \eqref{eq3_3}
in $\partial_d I^P_{h,T}u(X)=\sum_{i\in\mathcal{I}}u(M_i) \partial_d\phi^P_{i,T}(X), d = 1, 2$, we have
\begin{equation}
\begin{split}
\label{eq4_42}
& \partial_dI^P_{h,T}u(X) = u^s(X)\sum_{i\in\mathcal{I}}\partial_d\phi^P_{i,T}(X)+ \nabla u^s(X)\cdot \sum_{i\in\mathcal{I}}(M_i-X)\partial_d\phi^P_{i,T}(X)\\
&~~+ \sum_{i\in \mathcal{I}^{s'}}\left( \left( M^s(\widetilde{Y}_i)-I \right)\nabla u^s(X)\right) \cdot (M_i-\widetilde{Y}_i)\partial_d\phi^P_{i,T}(X)+\sum_{i\in\mathcal{I}}R_i^s\partial_d\phi^P_{i,T}(X), ~\forall X \in T^s, ~s = \pm.
\end{split}
\end{equation}
Then, applying \eqref{eq4_26} and Theorem \ref{thm4_3} to \eqref{eq4_42} we have
%Note that $\sum_{i\in\mathcal{I}}\partial_x\phi^P_{i,T}\equiv0$ from \eqref{eq4_26}. And by the similar argument above, we apply \eqref{eq4_30} to get
\begin{equation*}
\begin{split}
%\label{ }
&\partial_dI^P_{h,T}u(X) = \nabla u^s(X)\cdot \mathbf{ e}_d - \sum_{i\in\mathcal{I}^{s'}}\left( \left(\overline{M}^s(F) - I\right)\nabla u^s(X) \right)\cdot(M_i-\overline{X}_i)\partial_d\phi^P_{i,T}(X)\\
&~~+\sum_{i\in \mathcal{I}^{s'}}\left( \left( M^s(\widetilde{Y}_i)-I \right)\nabla u^s(X)\right)\cdot (M_i-\widetilde{Y}_i)\partial_d\phi^P_{i,T}(X)+\sum_{i\in\mathcal{I}}R_i^s\partial_x\phi^P_{i,T}(X), ~\forall X \in T^s, ~s = \pm \\
\end{split}
\end{equation*}
which leads to \eqref{eq4_35b} because $\nabla u^s(X)\cdot \mathbf{ e}_d=\partial_d u^s(X)$ and $M_i-\overline{X}_i=(M_i-\tilde{Y}_i)+(\tilde{Y}_i-\overline{X}_i)$.
\end{proof}

%estimation
By an argument similar to that used in \cite{2013ZhangTHESIS}, we can estimate $E_i^s$ and $F_i^s$ in \eqref{eq4_37} by geometric properties established in Section 3.
\begin{lemma}
\label{lem4_5}
There exist constants $C>0$ independent of the interface location and $\beta^{\pm}$ such that the following estimates hold for
every $T\in\mathcal{T}^i_h$ and $u\in PC^2_{int}(T)$:
\begin{subequations}
\begin{equation}
\label{eq_lem4_5_1}
\| E_i^-\|_{0,T^-}\leqslant C\kappa h^2 |u|_{1,T^-}, \;\;\;\; i\in \mathcal{I}^{+}, ~~~~\| E_i^+\|_{0,T^+}\leqslant C\frac{\kappa}{\rho}h^2 |u|_{1,T^+}, \;\;\;\; i\in \mathcal{I}^{-}.
\end{equation}
\begin{equation}
\label{eq_lem4_5_2}
\| F_i^-\|_{0,T^-}\leqslant Ch^2 |u|_{1,T^-}, \;\;\;\; i\in \mathcal{I}^{+}, ~~~~\| F_i^+\|_{0,T^+}\leqslant \frac{C}{\rho}h^2 |u|_{1,T^+}, \;\;\;\; i\in \mathcal{I}^{-}.
\end{equation}
\end{subequations}
%\begin{subequations}\label{eq_lem4_5_1}
%\begin{align}
%&\| E_i^s\|_{0,T^s}\leqslant Ch^2 |u|_{1,T^s}, \;\;\;\; i\in \mathcal{I}^{s'}, \label{eq_lem4_5_1a} \\
%&\| F_i^s\|_{0,T^s}\leqslant Ch^2 |u|_{1,T^s}, \text{~~when~~} \overline{X}_i=\widetilde{Y}_{i\bot}, \;\;\;\; i\in \mathcal{I}^{s'}. \label{eq_lem4_5_1b}
%\end{align}
%\end{subequations}
\end{lemma}
\begin{proof}
By $\|M_i-\widetilde{Y}_i\|\leqslant h$ and Lemma \ref{lem2_6}, we have
\begin{equation}
%\label{ }
\|E_i^s\|_{0,T^s}\leqslant\|(M^s(\widetilde{Y}_i)-\overline{M}^s(F))\| \; \|\nabla u^s\|_{0,T^s} \; \|(M_i-\widetilde{Y}_j)\|
\end{equation}
which implies the \eqref{eq_lem4_5_1} by \eqref{eq2_13} in Lemma \ref{lem2_6}. By Lemma \ref{lem2_5}, $F_i^s$ is independent of the choice of $\overline{X}_i \in l$. Hence,
by taking $\overline{X}_i=\widetilde{Y}_{i\bot}$ in $F_i^s$ and applying Lemma \ref{lem2_3}, we have
\begin{equation}
%\label{ }
\|F^s_i\|_{0,T^s}\leqslant\|\overline{M}^s(F) - I\| \; \|\nabla u^s\|_{0,T^s} \; \|\widetilde{Y}_i-\widetilde{Y}_{i\bot}\|
\end{equation}
which establishes \eqref{eq_lem4_5_2} by \eqref{bounds_for_Mbar}.
\end{proof}

Now we are ready to prove the main result in this section.
\begin{thm}
\label{thm4_6}
Assume all the conditions required by Theorem \ref{thm4_2} hold and $u\in PH^2_{int}(T)$. Then on every $T\in\mathcal{T}^i_h$ the following hold.
\begin{itemize}
\item
For rotated-$Q_1$ and Crouzeix-Raviart finite elements,
\begin{subequations}
\label{eq4_47}
\begin{equation}
\label{eq4_47_a}
\| I^P_{h,T}u-u\|_{0,T^-} + h| I^P_{h,T}u-u |_{1,T^-}\leqslant C\frac{1+\kappa}{\rho}h^2 (|u|_{1,T}+|u|_{2,T}),
\end{equation}
\begin{equation}
\label{eq4_47_a}
\| I^P_{h,T}u-u\|_{0,T^+} + h| I^P_{h,T}u-u |_{1,T^+}\leqslant C\frac{1+\kappa}{\rho^2}h^2 (|u|_{1,T}+|u|_{2,T}),
\end{equation}
\end{subequations}
where $C$ depends on $\lambda$ chosen for rotated-$Q_1$ case in \eqref{rem_unisol_eq1_1}.
\item
For linear finite elements,
\begin{subequations}
\label{linear_local_est}
\begin{equation}
\label{linear_local_est_a}
\| I^P_{h,T}u-u\|_{0,T^-} +  h| I^P_{h,T}u-u |_{1,T^-}\leqslant C\frac{1+\kappa}{\rho}h^2 (|u|_{1,T}+|u|_{2,T}),
\end{equation}
\begin{equation}
\label{linear_local_est_b}
\| I^P_{h,T}u-u\|_{0,T^+} + h| I^P_{h,T}u-u |_{1,T^+}\leqslant C\frac{1+\kappa}{\rho^{3/2}}h^2 (|u|_{1,T}+|u|_{2,T}).
\end{equation}
\end{subequations}

\item
For bilinear finite elements,
\begin{equation}
\label{bilinear_local_est}
\| I^P_{h,T}u-u\|_{0,T} +  h| I^P_{h,T}u-u |_{1,T}\leqslant C \frac{1+\kappa}{\rho}h^2 (|u|_{1,T}+|u|_{2,T}).
\end{equation}
\end{itemize}
\end{thm}

\begin{proof}
On each $T \in \mathcal{T}_h^i$, by Theorem \ref{thm4_5}, for every $u\in PC^2_{int}(T)$ we have
\begin{equation}
\label{}
| I^P_{h,T}u-u |_{k,T^s}\leqslant \sum_{i \in\mathcal{I}^{s'}}\left(\|E^s_i\|_{0,T^s}+\|F^s_i\|_{0,T^s} \right)|\phi^P_{i,T}|_{k,\infty,T^s} +\sum_{i\in \mathcal{I}}\|R^s_i\|_{0,T^s}|\phi^P_{i,T}|_{k,\infty,T^s},
\end{equation}
where $k=0,1$.
%\begin{equation}
%%\label{ }
%\| \partial_d(I^P_{h,T}u-u)\|_{0,T^s}\leqslant \frac{C}{h}\left(\sum_{i\in\mathcal{I}^{s'}}\left(\|E^s_i\|_{0,T^s}+\|F^s_i\|_{0,T^s}\right)+\sum_{i\in \mathcal{I}}\|R^s_i\|_{0,T^s}\right),
%~d = 1, 2.
%\end{equation}
Then, applying Lemmas \ref{lem3_3}, \ref{lem3_5}, \ref{lem4_5} and Theorem \ref{thm4_2} for corresponding IFE functions to the inequality above,
%\begin{eqnarray*}
%\| I^P_{h,T}u-u\|_{0,T^s} + h| I^P_{h,T}u-u |_{1,T^s}\leqslant Ch^2 (|u|_{1,T}+|u|_{2,T}), ~s = \pm
%\end{eqnarray*}
we obtain \eqref{eq4_47}-\eqref{bilinear_local_est} for $u\in PC_{int}^2(T)$. Finally,
the density hypothesis \textbf{(H4)} shows that \eqref{eq4_47}-\eqref{bilinear_local_est} also hold for any $u\in PH^2_{int}(T)$.
\end{proof}

The local estimate in Theorem \ref{thm4_6} leads to the following global estimate for the IFE interpolation directly.
\begin{thm}
\label{thm5_7}
For any $u\in PH^2_{int}(\Omega)$, the following estimate of interpolation error holds
\begin{equation}
\label{eq4_49}
\| I^P_{h}u-u\|_{0,\Omega^s} + h| I^P_{h}u-u |_{1,\Omega^s}\leqslant \widetilde{C}^sh^2 (|u|_{1,\Omega}+|u|_{2,\Omega}), ~~ s=\pm.
\end{equation}
The constants $\widetilde{C}$ depending on $\kappa$ and $\rho$ are specified as the following:
\begin{itemize}
\item for the rotated-$Q_1$ and Crouzeix-Raviart IFE space,
\begin{equation}
\label{interp_estimates_Q1_CR}
\widetilde{C}^-=C\frac{1+\kappa}{\rho},~~~\widetilde{C}^+=C\frac{1+\kappa}{\rho^2},
\end{equation}
where $C$ depends on $\lambda$ for rotated-$Q_1$ case;
\item for the linear IFE space,
\begin{equation}
\label{interp_estimates_linear}
\widetilde{C}^-=C\frac{1+\kappa}{\rho},~~~\widetilde{C}^+=C\frac{1+\kappa}{\rho^{3/2}};
\end{equation}
\item for the bilinear IFE space,
\begin{equation}
\label{interp_estimates_linear}
\widetilde{C}^-=\widetilde{C}^+=C\frac{1+\kappa}{\rho}.
\end{equation}
\end{itemize}
\begin{proof}
Estimate \eqref{eq4_49} follows directly from combining the estimates \eqref{eq4_47} to \eqref{bilinear_local_est} and \eqref{eq4_32}.
\end{proof}
\end{thm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%section



\section{Numerical Examples}

In this section we use numerical examples to demonstrate the approximation capability of the IFE spaces by IFE interpolation and IFE solutions. In generating numerical results, all computations involving integrations on interface elements, such as the assemblage
of local matrices and vectors with IFE shape functions or assessing the errors with integral norms, are handled by the numerical quadratures based on
the transfinite mapping between the reference straight edge triangle/square and the physical curved edge triangles/quadrilaterals. More
details about quadrature techniques on curved-edge domains can be found in \cite{2009Kopriva,2011SevillaSonia}.

All the numerical results to be presented are generated in the domain $\Omega = (-1,1)\times (-1,1)$ in which the interface curve $\Gamma$ is a circle with radius $r_0=\pi/6.28$ which divides $\Omega$ into two subdomains $\Omega^-$ and $\Omega^+$ with
$$
\Omega^-=\{ (x,y): x^2+y^2 < r_0^2 \}.
$$
The function to be approximated is
\begin{equation}
\label{example_u}
u(x,y) =\left\{
\begin{aligned}
&\frac{1}{\beta^-}r^{\alpha},  &\; (x,y)\in \Omega^-,\\
&\frac{1}{\beta^+}r^{\alpha}+\left( \frac{1}{\beta^-} -\frac{1}{\beta^+} \right)r_0^{\alpha}, & \; (x,y) \in \Omega^+,
\end{aligned}
\right.
\end{equation}
where $r=\sqrt{x^2+y^2}$ and $\alpha=5$. It is easy to verify that $u$ satisfies the interface jump condition \eqref{eq1_3} and \eqref{eq1_4}. We note that this is the same interface problem for the numerical examples given in \cite{2015LinLinZhang}. Numerical examples presented here are generated with the bilinear IFE space developed in Section \ref{sec:IFE Spaces and Their Properties}, and we note that numerical results with other IFE spaces developed in Section \ref{sec:IFE Spaces and Their Properties} are similar which are therefore not presented in order to avoid redundancy.

Note that the curvature of the interface in this interface problem is $\kappa\approx2$.
Condition \eqref{rem_unisol_eq2_1} allows us to use $\bar{\kappa} \approx 0.031$. Then using $\epsilon \approx 0.4$ in
\eqref{h_small_assump} leads to a suggested bound for the mesh size $h \approx 0.0273$. Therefore, our numerical experiments presented in this
section are all on meshes whose sizes are not larger than $1/40 = 0.025$ which can sufficiently satisfy the conditions in the error estimation in the previous section. \\
\commentout{
clear
kappa_b = 0.031
lambda = 1 - 4*sqrt(2)*sqrt(kappa_b)
kappa = 2
ep = 0.4
[sqrt(kappa_b)/(sqrt(2)*(1 + (1-2*ep^2)^(-3/2)*kappa)), ep/kappa]
1./[20, 40]
}

%\setcounter{table}{1}

\noindent
{\bf The convergence of IFE interpolation}: Table \ref{table:bilinearIterpolationError_1_10000} and Table \ref{table:bilinearIterpolationError_10000_1} present interpolation error
$u - I_hu$ in both the $L^2$ and the semi-$H^1$ norms over a sequence of meshes whose mesh size is $h$. In these tables, the rate is the estimated
values of $r$ such that $\|u - I_hu\|_{0,\Omega} \approx Ch^r$ or $|u - I_hu|_{1,\Omega} \approx Ch^r$ with numerical results generated on two consecutive meshes. The estimated
values for $r$ clearly demonstrate the optimal convergence of $I_hu$. We note that this example involves a coefficient $\beta$ with a jump quite large. Our numerical experiments show that these IFE spaces converge optimally also when $\beta$ has a moderate jump such as $\beta^- = 1, \beta^+ = 10$.

%Here we only conisder two case that one is the moderate data $\beta^-=1$, $\beta^+=10$ and another one is the challenging data $\beta^-=1$, $\beta^+=10000$ and $\alpha=5$.

%\begin{table}[H]
%\begin{center}
%\begin{tabular}{|c |c c|c c|}
%\hline
%$N$ & $\|\cdot\|_{0,\Omega}$ & rate & $|\cdot|_{1,\Omega}$ & rate  \\
%\hline
%10  &  1.4712E-2 &            & 1.8003E-1 &  \\
%\hline
%20 & 3.7703E-3 & 1.9642 &  9.1797E-2 & 0.9718 \\
%\hline
%40 & 9.4937E-4 & 1.9896 &  4.6187E-2 & 0.9910\\
%\hline
%80 & 2.3829E-4 & 1.9943 &  2.3170E-2 & 0.9952\\
%\hline
%160 & 5.9702E-5 & 1.9968 &  1.1606E-2 &0.9974\\
%\hline
%320 & 1.4944E-5 & 1.9982  &  5.8079E-3 & 0.9988\\
%\hline
%640 & 3.7382E-6 & 1.9992  &  2.9053E-3 & 0.9993\\
%\hline
%1280 & 9.3482E-7 & 1.9996  &  1.4530E-3 & 0.9997\\
%\hline
%\end{tabular}
%\caption{Interpolation errors and rates for IFE function with $\beta^-=1$ and $\beta^+=10$}
%\end{center}
%\end{table}


\begin{table}[H]
\begin{center}
\begin{tabular}{|c |c c|c c|}
\hline
$h$    & $\|u - I_hu\|_{0,\Omega}$ & rate   & $|u - I_hu|_{1,\Omega}$ & rate   \\ \hline
%10  & 2.965E-3                  &        & 4.5040E-2               &        \\ \hline
%1/20  & 1.0260E-3                 &        & 2.7724E-2               &        \\ \hline
1/40   & 2.7681E-4                 &        & 1.4482E-2               & \\ \hline
1/80   & 7.2447E-4                 & 1.9339 & 7.4468E-3               & 0.9596 \\ \hline
1/160  & 1.8580E-5                 & 1.9632 & 3.7827E-3               & 0.9772 \\ \hline
1/320  & 4.7122E-6                 & 1.9793 & 1.9061E-3               & 0.9888 \\ \hline
1/640  & 1.1858E-6                 & 1.9906 & 9.5723E-4               & 0.9937 \\ \hline
1/1280 & 2.9744E-7                 & 1.9952 & 4.7965E-4               & 0.9969 \\ \hline
\end{tabular}
\end{center}
\caption{Interpolation errors and rates for the bilinear IFE function, $\beta^-=1$ and $\beta^+=10000$}
\label{table:bilinearIterpolationError_1_10000}
\end{table}

\begin{table}[H]
\begin{center}
\begin{tabular}{|c |c c|c c|}
\hline
$h$    & $\|u - I_hu\|_{0,\Omega}$ & rate   & $|u - I_hu|_{1,\Omega}$ & rate   \\ \hline
%10  & 2.965E-3                  &        & 4.5040E-2               &        \\ \hline
%1/20  & 3.6177E-2                 &        & 8.7489E-1               &        \\ \hline
1/40   & 9.0663E-3                 &        & 4.3850E-1               &  \\ \hline
1/80   & 2.2680E-3                 & 1.9991 & 2.1939E-1               & 0.9991 \\ \hline
1/160  & 5.6711E-4                 & 1.9997 & 1.0971E-1               & 0.9998 \\ \hline
1/320  & 1.4179E-4                 & 1.9999 & 5.4859E-2               & 0.9999 \\ \hline
1/640  & 3.5447E-5                 & 2.0000 & 2.7430E-2               & 1.0000 \\ \hline
1/1280 & 8.8618E-6                 & 2.0000 & 1.3715E-2               & 1.0000 \\ \hline
\end{tabular}
\end{center}
\caption{Interpolation errors and rates for the bilinear IFE function, $\beta^-=10000$ and $\beta^+=1$}
\label{table:bilinearIterpolationError_10000_1}
\end{table}

\noindent
{\bf The convergence of the IFE solution}: Let $u_h$ be the IFE solution generated by the bilinear IFE space applied in the partially penalized method in
\cite{2015LinLinZhang} for the interface problem \eqref{eq1_1}-\eqref{eq1_4} where $f$ and $g$ are chosen such that $u$ given by \eqref{example_u} is its exact solution. The errors in the bilinear IFE solution generated by the symmetric partially penalized IFE (SPP IFE) method on a sequence of meshes are listed in
Tables \ref{table:bilinearIFESolutionError_SymPen1_10000} and \ref{table:bilinearIFESolutionError_SymPen10000_1}. The values of numerically estimated rate $r$ in these tables clearly indicate the optimal convergence of the bilinear IFE solution gauged in either the $L^2$ norm or $H^1$ norm. We also have carried out extensive numerical experiments by applying the IFE spaces developed in Section
\ref{sec:IFE Spaces and Their Properties} to the partially penalized IFE methods in \cite{2015LinLinZhang} with all the popular penalties, and we have
observed similar optimal convergence in the related IFE solution for this interface problem.

%\begin{table}[H]
%\begin{center}
%\begin{tabular}{| c |c c|c c|c c|c c|}
%\hline
%      &  \multicolumn{2}{|c|}{Classic IFE} & \multicolumn{2}{|c|}{SPP IFE} & \multicolumn{2}{|c|}{IPP IFE} & \multicolumn{2}{|c|}{NPP IFE} \\
%\hline
%$N$ & $\|\cdot\|_{0,\infty,\Omega}$ & rate & $\|\cdot\|_{0,\infty,\Omega}$ & rate & $\|\cdot\|_{0,\infty,\Omega}$ & rate & $\|\cdot\|_{0,\infty,\Omega}$ & rate  \\
%\hline
%10  &  1.7083E-2  &            & 1.6921E-2   &            & 1.6847E-2 &           & 1.6648E-2 &    \\
%\hline
%20 &  4.4551E-3 & 1.9116 &  4.4376E-3 & 1.9309 & 4.4416E-3 & 1.9233 & 4.4372E-3 & 1.9076  \\
%\hline
%40 & 1.0966E-3 & 2.0225 &  1.1088E-3  & 2.0008 & 1.1072E-3 & 2.0041 & 1.0961E-3 & 2.0173 \\
%\hline
%80 & 2.6894E-4 & 2.0266 &  2.7675E-4  & 2.0023 & 2.7596E-4 & 2.0045 & 2.7221E-4 & 2.0095 \\
%\hline
%160 & 6.6147E-5 & 2.0245 &  7.0911E-5 & 1.9645 &  6.8817E-5 & 2.0036 & 6.9427E-5 & 1.9711 \\
%\hline
%320 & 1.6481E-5 & 2.0049  &  1.7879E-5 & 1.9877 & 1.7387E-5 & 1.9847 & 1.7581E-5 & 1.9815 \\
%\hline
%640 & 4.3519E-6 & 1.9211  &  4.3539E-6 & 2.0379 & 4.3356E-6 & 2.0037 & 4.2737E-6 & 2.0404 \\
%\hline
%1280 & 1.3858E-6 & 1.6509  & 1.0865E-6  & 2.0026 & 1.1038E-6 & 1.9737 & 1.0882E-6 & 1.9736 \\
%\hline
%\end{tabular}
%\caption{Comparison of $\|u_h-u\|_{L^2(\Omega)}$ for different IFE methods with $\beta^-=1$, $\beta^+=10$. }
%\end{center}
%\end{table}

%\begin{table}[H]
%\begin{center}
%\begin{tabular}{| c |c c|c c|c c|c c|}
%\hline
%      &  \multicolumn{2}{|c|}{Classic IFE} & \multicolumn{2}{|c|}{SPP IFE} & \multicolumn{2}{|c|}{IPP IFE} & \multicolumn{2}{|c|}{NPP IFE} \\
%\hline
%$N$ & $|\cdot|_{1,\infty,\Omega}$ & rate & $|\cdot|_{1,\infty,\Omega}$ & rate & $|\cdot|_{1,\infty,\Omega}$ & rate & $|\cdot|_{1,\infty,\Omega}$ & rate  \\
%\hline
%10  &  1.8028E-1  &            & 1.7987E-1   &            & 1.7995E-1 &           & 1.7959E-1 &  \\
%\hline
%20 &  9.1823E-2 & 0.9716 &  9.1817E-2 & 0.9702 & 9.1818E-2 & 0.9707 & 9.1777E-2 & 0.9685 \\
%\hline
%40 & 4.6387E-2 & 0.9851 &  4.6196E-2  & 0.9910 & 4.6196E-2 & 0.9910 & 4.6183E-2 & 0.9908 \\
%\hline
%80 & 2.3469E-2 & 0.9830 &  2.3173E-2  & 0.9953 & 2.3173E-2 & 0.9953 & 2.3169E-2 & 0.9952 \\
%\hline
%160 & 1.1885E-2 & 0.9816 &  1.1609E-2 & 0.9972 &  1.1607E-2 & 0.9975 & 1.1608E-2 & 0.9971 \\
%\hline
%320 & 6.0182E-3 & 0.9818  &  5.8100E-3 & 0.9987 & 5.8080E-3 & 0.9989 & 5.8097E-3 & 0.9985 \\
%\hline
%640 & 3.1481E-3 & 0.9349  &  2.9053E-3 & 0.9998 & 2.9053E-3 & 0.9993 & 2.9053E-3 & 0.9998 \\
%\hline
%1280 & 1.6677E-3 & 0.9166  &  1.4530E-3 & 0.9997 & 1.4538E-3 & 0.9989 & 1.4542E-3 & 0.9984 \\
%\hline
%\end{tabular}
%\caption{ Comparison of $|u_h-u|_{H^1(\Omega)}$ for different IFE methods with $\beta^-=1$, $\beta^+=10$. }
%\end{center}
%\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{table}[H]
\begin{center}
\begin{tabular}{| c |c c|c c|}\hline
       & \multicolumn{2}{|c|}{SPP IFE}   & \multicolumn{2}{|c|}{SPP IFE}     \\ \hline
$h$    & $\|u-u_h\|_{0,\Omega}$ & rate   & $|u-u_h|_{1,\Omega}$    & rate    \\ \hline
%1/20  & 1.4807E-3              &        & 2.9032E-2               &         \\ \hline
1/40   & 3.7917E-4              &        & 1.5276E-2               &   \\ \hline
1/80   & 1.0409E-4              & 1.8650 & 7.9599E-3               & 0.9405  \\ \hline
1/160  & 2.5628E-5              & 2.0220 & 3.9096E-3               & 1.0257  \\ \hline
1/320  & 6.6828E-6              & 1.9392 & 1.9501E-3               & 1.0035  \\ \hline
1/640  & 1.7806E-6              & 1.9081 & 9.7745E-4               & 0.9964  \\ \hline
1/1280 & 4.0278E-7              & 2.1443 & 4.8374E-4               & 1.0148  \\ \hline
\end{tabular}
\end{center}
\caption{Errors in the bilinear IFE solution generated by the symmetric partially penalized IFE method, $\beta^-=1$, $\beta^+=10000$.}
\label{table:bilinearIFESolutionError_SymPen1_10000}
\end{table}


\begin{table}[H]
\begin{center}
\begin{tabular}{| c |c c|c c|}\hline
       & \multicolumn{2}{|c|}{SPP IFE}   & \multicolumn{2}{|c|}{SPP IFE}     \\ \hline
$h$    & $\|u-u_h\|_{0,\Omega}$ & rate   & $|u-u_h|_{1,\Omega}$    & rate    \\ \hline
%1/20  & 4.2019E-2              &        & 8.7973E-1               &         \\ \hline
1/40   & 1.0734E-2              &        & 4.4052E-1               &   \\ \hline
1/80   & 2.5715E-3              & 2.0616 & 2.1966E-1               & 1.0040  \\ \hline
1/160  & 6.2918E-4              & 2.0310 & 1.0974E-1               & 1.0012  \\ \hline
1/320  & 1.5709E-4              & 2.0019 & 5.4864E-2               & 1.0001  \\ \hline
1/640  & 4.0137E-5              & 1.9686 & 2.7431E-2               & 1.0000  \\ \hline
1/1280 & 9.8101E-6              & 2.0326 & 1.3715E-2               & 1.0000  \\ \hline
\end{tabular}
\end{center}
\caption{Errors in the bilinear IFE solution generated by the symmetric partially penalized IFE method, $\beta^-=10000$, $\beta^+=1$.}
\label{table:bilinearIFESolutionError_SymPen10000_1}
\end{table}
%\begin{table}[H]
%\begin{center}
%\begin{tabular}{| c |c c|c c|c c|}\hline
%       & \multicolumn{2}{|c|}{SPP IFE} & \multicolumn{2}{|c|}{IPP IFE} & \multicolumn{2}{|c|}{NPP IFE} \\ \hline
%$h$    & $|u-u_h|_{1,\Omega}$ & rate   & $|u-u_h|_{1,\Omega}$ & rate   & $|u-u_h|_{1,\Omega}$ & rate  \\ \hline
%1/10   & 2.9032E-2            &        & 2.9110E-2            &        & 2.9029E-2            &        \\ \hline
%1/20   & 1.5276E-2            & 0.9263 & 1.5276E-2            & 0.9303 & 1.5227E-2            & 0.9308 \\ \hline
%1/40   & 7.9599E-3            & 0.9405 & 7.9612E-3            & 0.9402 & 7.8361E-3            & 0.9585 \\ \hline
%1/80   & 3.9096E-3            & 1.0257 & 3.9106E-3            & 1.0256 & 3.8655E-3            & 1.0195 \\ \hline
%1/160  & 1.9501E-3            & 1.0035 & 1.9506E-3            & 1.0034 & 1.9266E-3            & 1.0046 \\ \hline
%1/320  & 9.7745E-4            & 0.9964 & 9.7774E-4            & 0.9964 & 9.6178E-4            & 1.0022 \\ \hline
%1/640  & 4.8374E-4            & 1.0148 & 4.8388E-4            & 1.0148 & 4.8055E-4            & 1.0010 \\ \hline
%\end{tabular}
%\caption{ Comparison of $|u_h-u|_{H^1(\Omega)}$ for different IFE methods with $\beta^-=1$, $\beta^+=10000$. }
%\end{center}
%\label{table:bilinearIFESolutionH1Error_1_10000}
%\end{table}


\begin{thebibliography}{10}

\bibitem{2015AdjeridChaabaneLin}
Slimane Adjerid, Nabil Chaabane, and Tao Lin.
\newblock An immersed discontinuous finite element method for stokes interface
  problems.
\newblock {\em Comput. Methods Appl. Mech. Engrg.}, 293:170--190, 2015.
%\newblock in press.

\bibitem{1970Babuska}
Ivo Babu{\v{s}}ka.
\newblock The finite element method for elliptic equations with discontinuous
  coefficients.
\newblock {\em Computing (Arch. Elektron. Rechnen)}, 5:207--213, 1970.

\bibitem{2003BabuskaBanerjeeOsborn}
Ivo Babu{\v{s}}ka, Uday Banerjee, and John~E. Osborn.
\newblock Survey of meshless and generalized finite element methods: a unified
  approach.
\newblock {\em Acta Numer.}, 12:1--125, 2003.

\bibitem{1998BabuskaZhang}
Ivo Babu{\v{s}}ka and Zhimin Zhang.
\newblock The partition of unity method for the elastically supported beam.
\newblock {\em Comput. Methods Appl. Mech. Engrg.}, 152(1-2):1--18, 1998.
\newblock Symposium on Advances in Computational Mechanics, Vol. 5 (Austin, TX,
  1997).

\bibitem{1987BarrettElliott}
John~W. Barrett and Charles~M. Elliott.
\newblock Fitted and unfitted finite-element methods for elliptic equations
  with smooth interfaces.
\newblock {\em IMA J. Numer. Anal.}, 7(3):283--300, 1987.

\bibitem{2009BastianEngwer}
Peter Bastian and Christian Engwer.
\newblock An unfitted finite element method using discontinuous {G}alerkin.
\newblock {\em Internat. J. Numer. Methods Engrg.}, 79(12):1557--1576, 2009.

\bibitem{2010Bordas}
Stephane Pierre~Alain Bordas, T.~Rabczuk, J.J. Rodenas, Pierre Kerfriden,
  M.~Moumnassi, and S.~Belouettar.
\newblock Recent advances towards reducing the meshing and re-meshing burden in
  computational sciences.
\newblock {\em Computational technology reviews}, 2:51--82, 2010.

\bibitem{2001Braess}
Dietrich Braess.
\newblock {\em Finite elements}.
\newblock Cambridge University Press, Cambridge, second edition, 2001.
\newblock Theory, fast solvers, and applications in solid mechanics, Translated
  from the 1992 German edition by Larry L. Schumaker.

\bibitem{1996BrambleKing}
James~H. Bramble and J.~Thomas King.
\newblock A finite element method for interface problems in domains with smooth
  boundaries and interfaces.
\newblock {\em Adv. Comput. Math.}, 6(2):109--138, 1996.

\bibitem{1994BrennerScott}
Susanne~C. Brenner and L.~Ridgway Scott.
\newblock {\em The mathematical theory of finite element methods}, volume~15 of
  {\em Texts in Applied Mathematics}.
\newblock Springer-Verlag, New York, 1994.

\bibitem{2015BurmanClaus}
Erik Burman, Susanne Claus, Peter Hansbo, Mats~G. Larson, Andr{\'e} Massing,
  Produktutveckling JTH, Tekniska H{\"o}gskolan, H{\"o}gskolan
  i~J{\"o}nk{\"o}ping, and JTH. Forskningsmilj{\"o}
  Produktutveckling~Simulering och optimering.
\newblock Cutfem: Discretizing geometry and partial differential equations.
\newblock {\em International Journal for Numerical Methods in Engineering},
  104(7):472--501, 2015.

\bibitem{2016BurmanGuzmanSanchezSarkisHcontrast}
Erik Burman, Johnny Guzm\'{a}n, Manuel~A. S\'{a}nchez, and Marcus Sarkis.
\newblock Robust flux error estimation of nitsche's method for high contrast
  interface problems.
\newblock {\em arXiv:1602.00603v1}, 2016.

\bibitem{1998ChenZou}
Zhiming Chen and Jun Zou.
\newblock Finite element methods and their convergence for elliptic and
  parabolic interface problems.
\newblock {\em Numer. Math.}, 79(2):175--202, 1998.

\bibitem{2010ChuGrahamHou}
C.-C. Chu, I.~G. Graham, and T.-Y. Hou.
\newblock A new multiscale finite element method for high-contrast elliptic
  interface problems.
\newblock {\em Math. Comp.}, 79(272):1915--1955, 2010.

\bibitem{1978Ciarlet}
Philippe~G. Ciarlet.
\newblock {\em The finite element method for elliptic problems}.
\newblock North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978.
\newblock Studies in Mathematics and its Applications, Vol. 4.

\bibitem{1966CloughTocher}
Ray.~W. Clough and James~L. Tocher.
\newblock Finite element stiffness matrices for analysis of plate bending.
\newblock In {\em Matrix Methods in Structual Mechanics}, pages 515--545, 1966.

\bibitem{1973CrouzeixRaviart}
M.~Crouzeix and P.-A. Raviart.
\newblock Conforming and nonconforming finite element methods for solving the
  stationary {S}tokes equations. {I}.
\newblock {\em Rev. Fran\c caise Automat. Informat. Recherche Op\'erationnelle
  S\'er. Rouge}, 7(R-3):33--75, 1973.

\bibitem{2001DolbowMoesBelytschko}
John Dolbow, Nicolas Mo{\"e}s, and Ted Belytschko.
\newblock An extended finite element method for modeling crack growth with
  frictional contact.
\newblock {\em Comput. Methods Appl. Mech. Engrg.}, 190(51-52):6825--6846,
  2001.

\bibitem{2009EfendievHou}
Yalchin Efendiev and Thomas~Y. Hou.
\newblock {\em Multiscale finite element methods}, volume~4 of {\em Surveys and
  Tutorials in the Applied Mathematical Sciences}.
\newblock Springer, New York, 2009.
\newblock Theory and applications.

\bibitem{2007GongLiLi}
Yan Gong, Bo~Li, and Zhilin Li.
\newblock Immersed-interface finite-element methods for elliptic interface
  problems with nonhomogeneous jump conditions.
\newblock {\em SIAM J. Numer. Anal.}, 46(1):472--495, 2007/08.

\bibitem{2009GuyomarchLeeJeon}
Gr{\'e}gory Guyomarc'h, Chang-Ock Lee, and Kiwan Jeon.
\newblock A discontinuous {G}alerkin method for elliptic interface problems
  with application to electroporation.
\newblock {\em Comm. Numer. Methods Engrg.}, 25(10):991--1008, 2009.

\bibitem{2015GuzmanSanchezSarkisP1}
Johnny Guzm\'{a}n, Manuel~A. S\'{a}nchez, and Marcus Sarkis.
\newblock A finite element method for high-contrast interface problems with
  error estimates independent of contrast.
\newblock {\em arXiv:1507.03873v2}, 2015.

\bibitem{2015GuzmanSanchezSarkisMathComp}
Johnny Guzm\'{a}n, Manuel~A. S\'{a}nchez, and Marcus Sarkis.
\newblock On the accuracy of finite element approximations to a class of
  interface problems.
\newblock {\em Mathematics of Computation}, 2015.
\newblock (to appear).

\bibitem{2002HansboHansbo}
Anita Hansbo and Peter Hansbo.
\newblock An unfitted finite element method, based on {N}itsche's method, for
  elliptic interface problems.
\newblock {\em Comput. Methods Appl. Mech. Engrg.}, 191(47-48):5537--5552,
  2002.

\bibitem{2009HeTHESIS}
Xiaoming He.
\newblock {\em Bilinear immersed finite elements for interface problems}.
\newblock PhD thesis, Virginia Polytechnic Institute and State University,
  2009.

\bibitem{2008HeLinLin}
Xiaoming He, Tao Lin, and Yanping Lin.
\newblock Approximation capability of a bilinear immersed finite element space.
\newblock {\em Numer. Methods Partial Differential Equations},
  24(5):1265--1300, 2008.

\bibitem{2009HeLinLin}
Xiaoming He, Tao Lin, and Yanping Lin.
\newblock A bilinear immersed finite volume element method for the diffusion
  equation with discontinuous coefficient.
\newblock {\em Commun. Comput. Phys.}, 6(1):185--202, 2009.

\bibitem{2013HeLinLinZhang}
Xiaoming He, Tao Lin, Yanping Lin, and Xu~Zhang.
\newblock Immersed finite element methods for parabolic equations with moving
  interface.
\newblock {\em Numer. Methods Partial Differential Equations}, 29(2):619--646,
  2013.

\bibitem{2012HouLiWangWang}
Songming Hou, Zhilin Li, Liqun Wang, and Wei Wang.
\newblock A numerical method for solving elasticity equations with interfaces.
\newblock {\em Commun. Comput. Phys.}, 12(2):595--612, 2012.

\bibitem{2005KafafyWangLin}
R.~Kafafy, J.~Wang, and T.~Lin.
\newblock A hybrid-grid immersed-finite-element particle-in-cell simulation
  model of ion optics plasma dynamics.
\newblock {\em Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms},
  12(Suppl. Vol. 12b):1--16, 2005.

\bibitem{2007KafafyWang}
Raed Kafafy and Joseph Wang.
\newblock Whole ion optics gridlet simulations using a hybrid-grid
  immersed-finite-element particle-in-cell code.
\newblock {\em J. Propulsion Power}, 23(1):59--68, 2007.

\bibitem{2009Kopriva}
David~A. Kopriva.
\newblock {\em Implementing spectral methods for partial differential
  equations: algorithms for scientists and engineers}.
\newblock Springer, 2009.

\bibitem{2010KwakWeeChang}
Do~Y. Kwak, Kye~T. Wee, and Kwang~S. Chang.
\newblock An analysis of a broken {$P_1$}-nonconforming finite element method
  for interface problems.
\newblock {\em SIAM J. Numer. Anal.}, 48(6):2117--2134, 2010.

\bibitem{2004LiLinLinRogers}
Z.~Li, T.~Lin, Y.~Lin, and R.~C. Rogers.
\newblock An immersed finite element space and its approximation capability.
\newblock {\em Numer. Methods Partial Differential Equations}, 20(3):338--367,
  2004.

\bibitem{1998Li}
Zhilin Li.
\newblock The immersed interface method using a finite element formulation.
\newblock {\em Appl. Numer. Math.}, 27(3):253--267, 1998.

\bibitem{2003LiLinWu}
Zhilin Li, Tao Lin, and Xiaohui Wu.
\newblock New {C}artesian grid methods for interface problems using the finite
  element formulation.
\newblock {\em Numer. Math.}, 96(1):61--98, 2003.

\bibitem{2005LiYang}
Zhilin Li and Xingzhou Yang.
\newblock An immersed finite element method for elasticity equations with
  interfaces.
\newblock In {\em Recent advances in adaptive computation}, volume 383 of {\em
  Contemp. Math.}, pages 285--298. Amer. Math. Soc., Providence, RI, 2005.

\bibitem{2012LianBordas}
H.~Lian, S.~P.~A. Bordas, R.~Sevilla, and R.~N. Simpson.
\newblock Recent developments in cad/analysis integration.
\newblock {\em arXiv:1210.8216v1}, 2012.

\bibitem{2011LinLinSunWang}
T.~Lin, Y.~Lin, W.-W. Sun, and Z.~Wang.
\newblock Immersed finite element methods for 4th order differential equations.
\newblock {\em J. Comput. Appl. Math.}, 235(13):3953--3964, 2011.

\bibitem{2001LinLinRogersRyan}
Tao Lin, Yanping Lin, Robert Rogers, and M.~Lynne Ryan.
\newblock A rectangular immersed finite element space for interface problems.
\newblock In {\em Scientific computing and applications ({K}ananaskis, {AB},
  2000)}, volume~7 of {\em Adv. Comput. Theory Pract.}, pages 107--114. Nova
  Sci. Publ., Huntington, NY, 2001.

\bibitem{2013LinLinZhang1}
Tao Lin, Yanping Lin, and Xu~Zhang.
\newblock A method of lines based on immersed finite elements for parabolic
  moving interface problems.
\newblock {\em Adv. Appl. Math. Mech.}, 5(4):548--568, 2013.

\bibitem{2015LinLinZhang}
Tao Lin, Yanping Lin, and Xu~Zhang.
\newblock Partially penalized immersed finite element methods for elliptic
  interface problems.
\newblock {\em SIAM J. Numer. Anal.}, 53(2):1121--1144, 2015.

\bibitem{2013LinSheenZhang}
Tao Lin, Dongwoo Sheen, and Xu~Zhang.
\newblock A locking-free immersed finite element method for planar elasticity
  interface problems.
\newblock {\em J. Comput. Phys.}, 247:228--247, 2013.

\bibitem{2002LinWang}
Tao Lin and Joseph Wang.
\newblock An immersed finite element electric field solver for ion optics
  modeling.
\newblock In {\em Proceedings of AIAA Joint Propulsion Conference},
  Indianapolis, IN, Jul 2002. AIAA.

\bibitem{1996BabuskaMelenk}
J.~M. Melenk and I.~Babu{\v{s}}ka.
\newblock The partition of unity finite element method: basic theory and
  applications.
\newblock {\em Comput. Methods Appl. Mech. Engrg.}, 139(1-4):289--314, 1996.

\bibitem{2016KihyoMoon_thesis}
Kihyo Moon.
\newblock {\em Immersed Discontinuous Galerkin Methods for Acoustic Wave
  Propagation in Inhomogeneous Media}.
\newblock PhD thesis, Virginia Tech, 2016.

\bibitem{2014MoumnassiFigueredo}
M.~Moumnassi, S.~P.~A. Bordas, R.~Figueredo, and P.~Sansen.
\newblock Analysis using higher-order xfem: implicit representation of
  geometrical features from a given parametric representation.
\newblock {\em Mechanics \& Industry}, 15(5):443--448, 2014.

\bibitem{2011MohammedSalim}
Mohammed Moumnassi, Salim Belouettar, {\'E}ric B{\'e}chet, St{\'e}phane P.~A.
  Bordas, Didier Quoirin, and Michel Potier-Ferry.
\newblock Finite element analysis on implicitly defined domains: An accurate
  representation based on arbitrary parametric surfaces.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  200(5):774--796, 2011.

\bibitem{2013NadalAlbelda}
E.~Nadal, J.~J. R{\'o}denas, J.~Albelda, M.~Tur, J.~E. Taranc{\'o}n, and F.~J.
  Fuenmayor.
\newblock Efficient finite element methodology based on cartesian grids:
  Application to structural shape optimization.
\newblock {\em Abstract and Applied Analysis}, pages 1--19, 2013.

\bibitem{2015NguyenAnitescu}
Vinh~P. Nguyen, Cosmin Anitescu, St{\'e}phane P.~A. Bordas, and Timon Rabczuk.
\newblock Isogeometric analysis: An overview and computer implementation
  aspects.
\newblock {\em Mathematics and Computers in Simulation}, 117:89--116, 2015.

\bibitem{1992RanacherTurek}
R.~Rannacher and S.~Turek.
\newblock Simple nonconforming quadrilateral {S}tokes element.
\newblock {\em Numer. Methods Partial Differential Equations}, 8(2):97--111,
  1992.

\bibitem{2011SevillaSonia}
Ruben Sevilla and Sonia Fern{\'a}ndez-M{\'e}ndez.
\newblock Numerical integration over 2d nurbs-shaped domains with applications
  to nurbs-enhanced fem.
\newblock {\em Finite Elements in Analysis \& Design}, 47(10):1209--1220, 2011.

\bibitem{2008SevillaSoniaHuerta}
Ruben Sevilla, Sonia Fern{\'a}ndez-M{\'e}ndez, and Antonio Huerta.
\newblock Nurbs-enhanced finite element method (nefem).
\newblock {\em International Journal for Numerical Methods in Engineering},
  76(1):56--83, 2008.

\bibitem{2010VallaghePapadopoulo}
Sylvain Vallagh{\'e} and Th{\'e}odore Papadopoulo.
\newblock A trilinear immersed finite element method for solving the
  electroencephalography forward problem.
\newblock {\em SIAM J. Sci. Comput.}, 32(4):2379--2394, 2010.

\bibitem{2016WangXiaoXu}
Fei Wang, Yuanming Xiao, and Jinchao Xu.
\newblock High-order extended finite element methods for solving interface
  problems.
\newblock {\em arXiv:1604.06171}, 2016.

\bibitem{1982Xu}
Jinchao Xu.
\newblock Estimate of the convergence rate of the finite element solutions to
  elliptic equation of second order with discontinuous coefficients.
\newblock {\em Natural Science Journal of Xiangtan University}, 1:1--5, 1982.

\bibitem{2003YangLiLi}
Xingzhou Yang, Bo~Li, and Zhilin Li.
\newblock The immersed interface method for elasticity problems with
  interfaces.
\newblock {\em Dyn. Contin. Discrete Impuls. Syst. Ser. A Math. Anal.},
  10(5):783--808, 2003.
\newblock Progress in partial differential equations (Pullman, WA, 2002).

\bibitem{2013ZhangTHESIS}
Xu~Zhang.
\newblock {\em Nonconforming {I}mmersed {F}inite {E}lement {M}ethods for
  {I}nterface {P}roblems}.
\newblock PhD thesis, Virginia Polytechnic Institute and State University,
  2013.

\end{thebibliography}

%\bibliographystyle{plain}
%\clearpage
%\bibliography{RuchiBib}

\end{document}


1612.00919.extracted_bib
{2015AdjeridChaabaneLin Slimane Adjerid, Nabil Chaabane, and Tao Lin. \newblock An immersed discontinuous finite element method for stokes interface problems. \newblock {\em Comput. Methods Appl. Mech. Engrg., 293:170--190, 2015. %\newblock in press.
1612.00919.extracted_bib
{1970Babuska Ivo Babu{\v{ska. \newblock The finite element method for elliptic equations with discontinuous coefficients. \newblock {\em Computing (Arch. Elektron. Rechnen), 5:207--213, 1970.
1612.00919.extracted_bib
{2003BabuskaBanerjeeOsborn Ivo Babu{\v{ska, Uday Banerjee, and John~E. Osborn. \newblock Survey of meshless and generalized finite element methods: a unified approach. \newblock {\em Acta Numer., 12:1--125, 2003.
1612.00919.extracted_bib
{1998BabuskaZhang Ivo Babu{\v{ska and Zhimin Zhang. \newblock The partition of unity method for the elastically supported beam. \newblock {\em Comput. Methods Appl. Mech. Engrg., 152(1-2):1--18, 1998. \newblock Symposium on Advances in Computational Mechanics, Vol. 5 (Austin, TX, 1997).
1612.00919.extracted_bib
{1987BarrettElliott John~W. Barrett and Charles~M. Elliott. \newblock Fitted and unfitted finite-element methods for elliptic equations with smooth interfaces. \newblock {\em IMA J. Numer. Anal., 7(3):283--300, 1987.
1612.00919.extracted_bib
{2009BastianEngwer Peter Bastian and Christian Engwer. \newblock An unfitted finite element method using discontinuous {Galerkin. \newblock {\em Internat. J. Numer. Methods Engrg., 79(12):1557--1576, 2009.
1612.00919.extracted_bib
{2010Bordas Stephane Pierre~Alain Bordas, T.~Rabczuk, J.J. Rodenas, Pierre Kerfriden, M.~Moumnassi, and S.~Belouettar. \newblock Recent advances towards reducing the meshing and re-meshing burden in computational sciences. \newblock {\em Computational technology reviews, 2:51--82, 2010.
1612.00919.extracted_bib
{2001Braess Dietrich Braess. \newblock {\em Finite elements. \newblock Cambridge University Press, Cambridge, second edition, 2001. \newblock Theory, fast solvers, and applications in solid mechanics, Translated from the 1992 German edition by Larry L. Schumaker.
1612.00919.extracted_bib
{1996BrambleKing James~H. Bramble and J.~Thomas King. \newblock A finite element method for interface problems in domains with smooth boundaries and interfaces. \newblock {\em Adv. Comput. Math., 6(2):109--138, 1996.
1612.00919.extracted_bib
{1994BrennerScott Susanne~C. Brenner and L.~Ridgway Scott. \newblock {\em The mathematical theory of finite element methods, volume~15 of {\em Texts in Applied Mathematics. \newblock Springer-Verlag, New York, 1994.
1612.00919.extracted_bib
{2015BurmanClaus Erik Burman, Susanne Claus, Peter Hansbo, Mats~G. Larson, Andr{\'e Massing, Produktutveckling JTH, Tekniska H{\"ogskolan, H{\"ogskolan i~J{\"onk{\"oping, and JTH. Forskningsmilj{\"o Produktutveckling~Simulering och optimering. \newblock Cutfem: Discretizing geometry and partial differential equations. \newblock {\em International Journal for Numerical Methods in Engineering, 104(7):472--501, 2015.
1612.00919.extracted_bib
{2016BurmanGuzmanSanchezSarkisHcontrast Erik Burman, Johnny Guzm\'{an, Manuel~A. S\'{anchez, and Marcus Sarkis. \newblock Robust flux error estimation of nitsche's method for high contrast interface problems. \newblock {\em arXiv:1602.00603v1, 2016.
1612.00919.extracted_bib
{1998ChenZou Zhiming Chen and Jun Zou. \newblock Finite element methods and their convergence for elliptic and parabolic interface problems. \newblock {\em Numer. Math., 79(2):175--202, 1998.
1612.00919.extracted_bib
{2010ChuGrahamHou C.-C. Chu, I.~G. Graham, and T.-Y. Hou. \newblock A new multiscale finite element method for high-contrast elliptic interface problems. \newblock {\em Math. Comp., 79(272):1915--1955, 2010.
1612.00919.extracted_bib
{1978Ciarlet Philippe~G. Ciarlet. \newblock {\em The finite element method for elliptic problems. \newblock North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978. \newblock Studies in Mathematics and its Applications, Vol. 4.
1612.00919.extracted_bib
{1966CloughTocher Ray.~W. Clough and James~L. Tocher. \newblock Finite element stiffness matrices for analysis of plate bending. \newblock In {\em Matrix Methods in Structual Mechanics, pages 515--545, 1966.
1612.00919.extracted_bib
{1973CrouzeixRaviart M.~Crouzeix and P.-A. Raviart. \newblock Conforming and nonconforming finite element methods for solving the stationary {Stokes equations. {I. \newblock {\em Rev. Fran\c caise Automat. Informat. Recherche Op\'erationnelle S\'er. Rouge, 7(R-3):33--75, 1973.
1612.00919.extracted_bib
{2001DolbowMoesBelytschko John Dolbow, Nicolas Mo{\"es, and Ted Belytschko. \newblock An extended finite element method for modeling crack growth with frictional contact. \newblock {\em Comput. Methods Appl. Mech. Engrg., 190(51-52):6825--6846, 2001.
1612.00919.extracted_bib
{2009EfendievHou Yalchin Efendiev and Thomas~Y. Hou. \newblock {\em Multiscale finite element methods, volume~4 of {\em Surveys and Tutorials in the Applied Mathematical Sciences. \newblock Springer, New York, 2009. \newblock Theory and applications.
1612.00919.extracted_bib
{2007GongLiLi Yan Gong, Bo~Li, and Zhilin Li. \newblock Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions. \newblock {\em SIAM J. Numer. Anal., 46(1):472--495, 2007/08.
1612.00919.extracted_bib
{2009GuyomarchLeeJeon Gr{\'egory Guyomarc'h, Chang-Ock Lee, and Kiwan Jeon. \newblock A discontinuous {Galerkin method for elliptic interface problems with application to electroporation. \newblock {\em Comm. Numer. Methods Engrg., 25(10):991--1008, 2009.
1612.00919.extracted_bib
{2015GuzmanSanchezSarkisP1 Johnny Guzm\'{an, Manuel~A. S\'{anchez, and Marcus Sarkis. \newblock A finite element method for high-contrast interface problems with error estimates independent of contrast. \newblock {\em arXiv:1507.03873v2, 2015.
1612.00919.extracted_bib
{2015GuzmanSanchezSarkisMathComp Johnny Guzm\'{an, Manuel~A. S\'{anchez, and Marcus Sarkis. \newblock On the accuracy of finite element approximations to a class of interface problems. \newblock {\em Mathematics of Computation, 2015. \newblock (to appear).
1612.00919.extracted_bib
{2002HansboHansbo Anita Hansbo and Peter Hansbo. \newblock An unfitted finite element method, based on {Nitsche's method, for elliptic interface problems. \newblock {\em Comput. Methods Appl. Mech. Engrg., 191(47-48):5537--5552, 2002.
1612.00919.extracted_bib
{2009HeTHESIS Xiaoming He. \newblock {\em Bilinear immersed finite elements for interface problems. \newblock PhD thesis, Virginia Polytechnic Institute and State University, 2009.
1612.00919.extracted_bib
{2008HeLinLin Xiaoming He, Tao Lin, and Yanping Lin. \newblock Approximation capability of a bilinear immersed finite element space. \newblock {\em Numer. Methods Partial Differential Equations, 24(5):1265--1300, 2008.
1612.00919.extracted_bib
{2009HeLinLin Xiaoming He, Tao Lin, and Yanping Lin. \newblock A bilinear immersed finite volume element method for the diffusion equation with discontinuous coefficient. \newblock {\em Commun. Comput. Phys., 6(1):185--202, 2009.
1612.00919.extracted_bib
{2013HeLinLinZhang Xiaoming He, Tao Lin, Yanping Lin, and Xu~Zhang. \newblock Immersed finite element methods for parabolic equations with moving interface. \newblock {\em Numer. Methods Partial Differential Equations, 29(2):619--646, 2013.
1612.00919.extracted_bib
{2012HouLiWangWang Songming Hou, Zhilin Li, Liqun Wang, and Wei Wang. \newblock A numerical method for solving elasticity equations with interfaces. \newblock {\em Commun. Comput. Phys., 12(2):595--612, 2012.
1612.00919.extracted_bib
{2005KafafyWangLin R.~Kafafy, J.~Wang, and T.~Lin. \newblock A hybrid-grid immersed-finite-element particle-in-cell simulation model of ion optics plasma dynamics. \newblock {\em Dyn. Contin. Discrete Impuls. Syst. Ser. B Appl. Algorithms, 12(Suppl. Vol. 12b):1--16, 2005.
1612.00919.extracted_bib
{2007KafafyWang Raed Kafafy and Joseph Wang. \newblock Whole ion optics gridlet simulations using a hybrid-grid immersed-finite-element particle-in-cell code. \newblock {\em J. Propulsion Power, 23(1):59--68, 2007.
1612.00919.extracted_bib
{2009Kopriva David~A. Kopriva. \newblock {\em Implementing spectral methods for partial differential equations: algorithms for scientists and engineers. \newblock Springer, 2009.
1612.00919.extracted_bib
{2010KwakWeeChang Do~Y. Kwak, Kye~T. Wee, and Kwang~S. Chang. \newblock An analysis of a broken {$P_1$-nonconforming finite element method for interface problems. \newblock {\em SIAM J. Numer. Anal., 48(6):2117--2134, 2010.
1612.00919.extracted_bib
{2004LiLinLinRogers Z.~Li, T.~Lin, Y.~Lin, and R.~C. Rogers. \newblock An immersed finite element space and its approximation capability. \newblock {\em Numer. Methods Partial Differential Equations, 20(3):338--367, 2004.
1612.00919.extracted_bib
{1998Li Zhilin Li. \newblock The immersed interface method using a finite element formulation. \newblock {\em Appl. Numer. Math., 27(3):253--267, 1998.
1612.00919.extracted_bib
{2003LiLinWu Zhilin Li, Tao Lin, and Xiaohui Wu. \newblock New {Cartesian grid methods for interface problems using the finite element formulation. \newblock {\em Numer. Math., 96(1):61--98, 2003.
1612.00919.extracted_bib
{2005LiYang Zhilin Li and Xingzhou Yang. \newblock An immersed finite element method for elasticity equations with interfaces. \newblock In {\em Recent advances in adaptive computation, volume 383 of {\em Contemp. Math., pages 285--298. Amer. Math. Soc., Providence, RI, 2005.
1612.00919.extracted_bib
{2012LianBordas H.~Lian, S.~P.~A. Bordas, R.~Sevilla, and R.~N. Simpson. \newblock Recent developments in cad/analysis integration. \newblock {\em arXiv:1210.8216v1, 2012.
1612.00919.extracted_bib
{2011LinLinSunWang T.~Lin, Y.~Lin, W.-W. Sun, and Z.~Wang. \newblock Immersed finite element methods for 4th order differential equations. \newblock {\em J. Comput. Appl. Math., 235(13):3953--3964, 2011.
1612.00919.extracted_bib
{2001LinLinRogersRyan Tao Lin, Yanping Lin, Robert Rogers, and M.~Lynne Ryan. \newblock A rectangular immersed finite element space for interface problems. \newblock In {\em Scientific computing and applications ({Kananaskis, {AB, 2000), volume~7 of {\em Adv. Comput. Theory Pract., pages 107--114. Nova Sci. Publ., Huntington, NY, 2001.
1612.00919.extracted_bib
{2013LinLinZhang1 Tao Lin, Yanping Lin, and Xu~Zhang. \newblock A method of lines based on immersed finite elements for parabolic moving interface problems. \newblock {\em Adv. Appl. Math. Mech., 5(4):548--568, 2013.
1612.00919.extracted_bib
{2015LinLinZhang Tao Lin, Yanping Lin, and Xu~Zhang. \newblock Partially penalized immersed finite element methods for elliptic interface problems. \newblock {\em SIAM J. Numer. Anal., 53(2):1121--1144, 2015.
1612.00919.extracted_bib
{2013LinSheenZhang Tao Lin, Dongwoo Sheen, and Xu~Zhang. \newblock A locking-free immersed finite element method for planar elasticity interface problems. \newblock {\em J. Comput. Phys., 247:228--247, 2013.
1612.00919.extracted_bib
{2002LinWang Tao Lin and Joseph Wang. \newblock An immersed finite element electric field solver for ion optics modeling. \newblock In {\em Proceedings of AIAA Joint Propulsion Conference, Indianapolis, IN, Jul 2002. AIAA.
1612.00919.extracted_bib
{1996BabuskaMelenk J.~M. Melenk and I.~Babu{\v{ska. \newblock The partition of unity finite element method: basic theory and applications. \newblock {\em Comput. Methods Appl. Mech. Engrg., 139(1-4):289--314, 1996.
1612.00919.extracted_bib
{2016KihyoMoon_thesis Kihyo Moon. \newblock {\em Immersed Discontinuous Galerkin Methods for Acoustic Wave Propagation in Inhomogeneous Media. \newblock PhD thesis, Virginia Tech, 2016.
1612.00919.extracted_bib
{2014MoumnassiFigueredo M.~Moumnassi, S.~P.~A. Bordas, R.~Figueredo, and P.~Sansen. \newblock Analysis using higher-order xfem: implicit representation of geometrical features from a given parametric representation. \newblock {\em Mechanics \& Industry, 15(5):443--448, 2014.
1612.00919.extracted_bib
{2011MohammedSalim Mohammed Moumnassi, Salim Belouettar, {\'Eric B{\'echet, St{\'ephane P.~A. Bordas, Didier Quoirin, and Michel Potier-Ferry. \newblock Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces. \newblock {\em Computer Methods in Applied Mechanics and Engineering, 200(5):774--796, 2011.
1612.00919.extracted_bib
{2013NadalAlbelda E.~Nadal, J.~J. R{\'odenas, J.~Albelda, M.~Tur, J.~E. Taranc{\'on, and F.~J. Fuenmayor. \newblock Efficient finite element methodology based on cartesian grids: Application to structural shape optimization. \newblock {\em Abstract and Applied Analysis, pages 1--19, 2013.
1612.00919.extracted_bib
{2015NguyenAnitescu Vinh~P. Nguyen, Cosmin Anitescu, St{\'ephane P.~A. Bordas, and Timon Rabczuk. \newblock Isogeometric analysis: An overview and computer implementation aspects. \newblock {\em Mathematics and Computers in Simulation, 117:89--116, 2015.
1612.00919.extracted_bib
{1992RanacherTurek R.~Rannacher and S.~Turek. \newblock Simple nonconforming quadrilateral {Stokes element. \newblock {\em Numer. Methods Partial Differential Equations, 8(2):97--111, 1992.
1612.00919.extracted_bib
{2011SevillaSonia Ruben Sevilla and Sonia Fern{\'andez-M{\'endez. \newblock Numerical integration over 2d nurbs-shaped domains with applications to nurbs-enhanced fem. \newblock {\em Finite Elements in Analysis \& Design, 47(10):1209--1220, 2011.
1612.00919.extracted_bib
{2008SevillaSoniaHuerta Ruben Sevilla, Sonia Fern{\'andez-M{\'endez, and Antonio Huerta. \newblock Nurbs-enhanced finite element method (nefem). \newblock {\em International Journal for Numerical Methods in Engineering, 76(1):56--83, 2008.
1612.00919.extracted_bib
{2010VallaghePapadopoulo Sylvain Vallagh{\'e and Th{\'eodore Papadopoulo. \newblock A trilinear immersed finite element method for solving the electroencephalography forward problem. \newblock {\em SIAM J. Sci. Comput., 32(4):2379--2394, 2010.
1612.00919.extracted_bib
{2016WangXiaoXu Fei Wang, Yuanming Xiao, and Jinchao Xu. \newblock High-order extended finite element methods for solving interface problems. \newblock {\em arXiv:1604.06171, 2016.
1612.00919.extracted_bib
{1982Xu Jinchao Xu. \newblock Estimate of the convergence rate of the finite element solutions to elliptic equation of second order with discontinuous coefficients. \newblock {\em Natural Science Journal of Xiangtan University, 1:1--5, 1982.
1612.00919.extracted_bib
{2003YangLiLi Xingzhou Yang, Bo~Li, and Zhilin Li. \newblock The immersed interface method for elasticity problems with interfaces. \newblock {\em Dyn. Contin. Discrete Impuls. Syst. Ser. A Math. Anal., 10(5):783--808, 2003. \newblock Progress in partial differential equations (Pullman, WA, 2002).
1612.00919.extracted_bib
{2013ZhangTHESIS Xu~Zhang. \newblock {\em Nonconforming {Immersed {Finite {Element {Methods for {Interface {Problems. \newblock PhD thesis, Virginia Polytechnic Institute and State University, 2013.