text:
\section{Introduction}
\label{intro}
Let $\Omega \subset \mathbb{R}^d$, $d \in \{2,3\}$, be a bounded polyhedral domain. Furthermore, we assume that $\Omega$ is convex if necessary. We consider the Poisson problem as follows. Find $u: \Omega \to \mathbb{R}$ such that
\begin{align}
\displaystyle
- \varDelta u = f \quad \text{in $\Omega$}, \quad u = 0 \quad \text{on $\partial \Omega$}, \label{intro1}
\end{align}
where $f \in L^2(\Omega)$ is a given function. This paper gives error estimates for the first-order Crouzeix--Raviart (CR) finite-element approximation on anisotropic meshes in three dimensions. Anisotropic meshes have different mesh sizes in different directions. The shape regularity assumption on triangulations $\mathbb{T}_h$ is no longer valid on these meshes; see for example \cite{Ape99}. Furthermore, we do not impose the maximum-angle condition proposed in \cite{BabAzi76} during mesh partitioning. In many instances, the discussion also relates to two dimensions. We therefore discuss the problem here as uniformly valid in an arbitrary number of dimensions.

CR finite error estimates for the non-homogeneous Dirichlet Poisson problem are known. Let $CR_{h0}^1$ be the CR finite-element space, to be defined in Section \ref{FEspaces}. Let $u \in H_0^1(\Omega)$ and $u_h^{CR} \in CR_{h0}^1$ be the exact and CR finite-element solutions, respectively. In \cite[Corollary 2.2]{Gud10}, adopting medius analysis, the estimate
\begin{align}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
\leq c_0 \left( \inf_{v_h \in CR_{h0}^1} | u - v_h |_{H^1(\mathbb{T}_h)} +  Osc_1(f) \right), \label{intro2}
\end{align}
is given, where $| \cdot |_{H^1(\mathbb{T}_h)}$ denotes the broken (piecewise) $H^1$-semi norm defined in Section \ref{regularmesh}, and $c_0$ a positive constant independent of $h$. Here, the oscillation $Osc_1(f)$ is expressed as
\begin{align*}
\displaystyle
Osc_1(f) := \left( \sum_{T \in \mathbb{T}_h} h_T^2
\left [ \inf_{\bar{f} \in \mathcal{P}^0(T)} \| f - \bar{f} \|^2_{L^2(T)} \right] \right)^{1/2},
\end{align*}
where $\mathcal{P}^0(T)$ denotes the space of constant functions on $T$. Suppose that $u \in H^2(\Omega)$ and oscillation $Osc_1(f)$ vanishes. Let $I_h u \in CR_{h0}^1$ be the nodal interpolation of $u$ at the midpoints of the faces. Then, from the standard interpolation error estimate (see for example \cite[Corollary 1.109]{ErnGue04}), we have
\begin{align*}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)} \leq c_0 | u - I_h u |_{H^1(\mathbb{T}_h)}
\leq c_1 h |u|_{H^2(\Omega)},
\end{align*}
where $c_1$ represents a positive constant independent of $h$ and $u$ but depending on the parameter of the simplicial mesh; see for example \cite[Definition 1.107]{ErnGue04}. This parameter is bounded if the simplicial mesh sequence is shape regular. However, the situation is different without the shape-regular condition. The aim of the present paper is to deduce an analogous error estimate on anisotropic finite-element meshes. Note that very flat elements might be included in the mesh sequence. In many papers reporting on such investigations, the maximum-angle condition instead of the shape-regular condition is imposed. However, the maximum-angle condition is not necessarily needed to obtain error estimates. Recently, in the two-dimensional instance, the CR finite-element analysis of the non-homogeneous Dirichlet-Poisson problem has been investigated under a more relaxed mesh condition, \cite{KobTsu18b}. The present paper extends previous research to a three-dimensional setting.

However, it may not be easy to use the estimate \eqref{intro2} on anisotropic finite-element meshes. To overcome this difficulty, we use the interpolation error estimates obtained in \cite{IshKobTsu}. In that paper, the CR and Raviart--Thomas (RT) interpolation errors are bounded in terms of $h$ and the new parameter $H$, see Corollary \ref{modify=coro2}, \ref{modify=coro3}.

The CR finite-element space is not in $H_0^1(\Omega)$. Hence, an error between the exact solution and the CR finite-element approximation solution with a $H^1$-broken seminorm is divided into two parts (\cite{Bre15,ErnGue04}). One is an approximation error that measures how well the exact solution is approximated by the CR finite-element functions, the other is a nonconformity error term. For the former, the CR interpolation error estimates (Corollary \ref{modify=coro2}) are used. In the latter, the standard scaling argument is often used to obtain the error estimates. However, in this way, we are unable to derive the correct order on anisotropic meshes. To overcome this difficulty, we shall use the lowest-order RT interpolation error estimates on anisotropic meshes (Corollary \ref{modify=coro3}). By this technique, we consequently have the error estimates in the $H^1$-broken seminorm (Theorem \ref{CR=thr4}) and the $L^2$ norm (Theorem \ref{CR=thr5}) on anisotropic meshes.

Furthermore, we present an error estimate for the first-order RT finite-element approximation of the Poisson problem \eqref{intro1} based on the dual mixed formulation (Theorem \ref{mix=thr7}). In the proof, we again use Corollary \ref{modify=coro3}. We again emphasise that we do not impose either the shape-regular or the maximum-angle condition during mesh partitioning.

We next present the equivalence of the enriched piecewise linear CR finite-element method introduced by \cite{HuMa15} and the first-order RT finite-element method. In two dimensions, the work \cite{ArnBre85} represents pioneering research. Marini \cite{Mar85} further found an expression relating RT and CR finite-element methods:
\begin{align}
\displaystyle
\bar{\sigma}_h^{RT}|_T &= \nabla \bar{u}_h^{CR} - \frac{f_T^0}{2} (x - x_T) \quad \text{on $T$}, \label{intro3}
\end{align}
where $T$ denotes a mesh element, $x_i$ ($i=1,2,3$) the vertices of triangle $T$, $x_T$ the barycentre of $T$ such that $x_T := \frac{1}{3}(x_1 + x_2 + x_3)$, and $\bar{\sigma}_h^{RT}$ and $ \bar{u}_h^{CR}$ respectively denote the RT and CR finite-element solutions with a given external piecewise-constant function $f_T^0$. It was recently proved \cite{HuMa15} that the enriched piecewise-linear CR finite-element method is identical to the first-order RT finite-element method for both the Poisson and Stokes problems in any number of dimensions. In the present paper, we extend Marini's results to three dimensions (Lemma \ref{RTCR=lem12}).

The remainder of the present paper is organised as follows. Section 2 introduces the weak form of the continuous problem \eqref{intro1}, the finite-element meshes, and finite-element spaces. Furthermore, we propose a parameter $H$. Section 3 introduces discrete settings of the CR finite-element method for \eqref{intro1} and proposes error estimates. Section 4 proves error estimates for the first-order RT finite-element method based on the dual mixed formulation of the Poisson problem. Section 5 gives the equivalence of the RT and CR finite-element problems. Finally, Section 6 presents numerical results obtained using the Lagrange P1 element and the first-order CR element.

\begin{comment}
Let $\Omega \subset \mathbb{R}^d$, $d \in \{2,3\}$, be a bounded polyhedral domain. Further, we assume that $\Omega$ is convex if necessary. We consider the Poisson problem as follows. Find $u: \Omega \to \mathbb{R}$ such that
\begin{align}
\displaystyle
- \varDelta u = f \quad \text{in $\Omega$}, \quad u = 0 \quad \text{on $\partial \Omega$}, \label{intro1}
\end{align}
where $f \in L^2(\Omega)$ is a given function. Let $\ell$ be a nonnegative integer. $H^{\ell}(\Omega)$ is a Hilbert space with the scalar product $(\varphi,\psi)_{H^{\ell}(\Omega)} := \sum_{|\beta| \leq \ell} (\partial^{\beta} \varphi , \partial^{\beta} \psi)_{L^2(\Omega)}$ and the norm $\|  \varphi\|_{H^{\ell}(\Omega)} := \sqrt{(\varphi,\varphi)_{H^{\ell}(\Omega)}}$. We set $L^2(\Omega) := H^0(\Omega)$ with $(\cdot,\cdot) := (\cdot,\cdot)_{H^0(\Omega)}$ and $\| \cdot \| := \| \cdot \|_{H^0(\Omega)}$.

This paper gives error estimates for the first-order Crouzeix-Raviart (CR) finite element approximation on anisotropic meshes in the three-dimensional case. Anisotropic meshes have different mesh sizes in different directions. The shape regularity assumption on triangulations $\mathbb{T}_h$ is no longer valid on these meshes; for example, see \cite{Ape99}. Further, we do not impose the maximum-angle condition proposed in \cite{BabAzi76} for the mesh partition. In many cases, the discussion also relates to two dimensions. We will therefore discuss the problem here as uniformly valid in any number of dimensions as possible.

A CR finite error estimate for the non-homogeneous Dirichlet Poisson problem are known as follows. Let $CR_{h0}^1$ be the CR finite element space defined in Section \ref{FEspaces}. Let $u \in H_0^1(\Omega)$ and $u_h^{CR} \in CR_{h0}^1$ be the exact and CR finite element solutions, respectively. In \cite[Corollary 2.2]{Gud10}, adopting medius analysis, the estimate
\begin{align}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
\leq c \left( \inf_{v_h \in CR_{h0}^1} | u - v_h |_{H^1(\mathbb{T}_h)} +  Osc_1(f) \right), \label{intro2}
\end{align}
is given, where $| \cdot |_{H^1(\mathbb{T}_h)}$ is the broken (piecewise) $H^1$-semi norm defined in Section \ref{regularmesh}, and $c$ is a positive constant independent of $h$. Here, the oscillation $Osc_1(f)$ is described by
\begin{align*}
\displaystyle
Osc_1(f) := \left( \sum_{T \in \mathbb{T}_h} h_T^2 \left [ \inf_{\bar{f} \in \mathcal{P}^0(T)} \| f - \bar{f} \|^2_{L^2(T)}  \right] \right)^{1/2},
\end{align*}
where $\mathcal{P}^0(T)$ is the piecewise constant space in $T$. Suppose that $u \in H^2(\Omega)$ and the oscillation $Osc_1(f)$ vanishes. Let $I_h u \in CR_{h0}^1$ be the nodal interpolation of $u$ at the midpoints of the faces. Then, from the standard interpolation error estimate (e.g., see \cite[Corollary 1.109]{ErnGue04}), we have
\begin{align*}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)} \leq | u - I_h u |_{H^1(\mathbb{T}_h)}
\leq c h |u|_{H^2(\Omega)},
\end{align*}
where $c$ is a positive constant independent of $h$ and $u$ but dependent on the parameter of the simplicial mesh; for example, see \cite[Definition 1.107]{ErnGue04}. This parameter is bounded if the simplicial mesh sequence is shape regular. However, the situation is different without the shape regular condition. The aim of the present paper is to deduce an analogous error estimate on anisotropic finite element meshes. Note that very flat elements might be included in the mesh sequence. In many papers, the maximum-angle condition instead of the shape regular condition is imposed to investigate. However, the maximum-angle condition is not necessary needed to get error estimates. Recently, in the two-dimensional case, the CR finite element analysis of the non-homogeneous Dirichlet Poisson problem has been investigated under the more relaxed mesh condition, \cite{KobTsu18b}. The present paper extends previous research to the three-dimensional case.

 It however might not be easy to use the estimate \eqref{intro1} on anisotropic finite element meshes. The CR finite element space is not in $H_0^1(\Omega)$. So, an error between the exact solution and the CR finite element approximation solution in $H^1$-broken seminorm is divided into two parts. One is an approximation error which measures how well the exact solution can be approximated by the CR finite element functions. The other is a nonconformity error term. For the former estimate, the CR interpolation error estimates are often used. In the latter estimate, the standard scaling argument is often used to obtain the error estimates. However, in this way, we cannot get the correct order on anisotropic meshes. To overcome this difficulty, we shall use the lowest order Raviart--Thomas (RT) interpolation error estimates on anisotropic meshes. By this technique, we consequently have the error estimates in the $H^1$-broken seminorm (Theorem \ref{CR=thr4}) and the $L^2$ norm (Theorem \ref{CR=thr5}) on anisotropic meshes.

 We next present the equivalence of the enriched piecewise linear CR finite element method introduced by \cite{HuMa15} and the first-order RT finite element method. In the two-dimensional case, the work \cite{ArnBre85} is pioneering research. Marini \cite{Mar85} further found the relation between the RT and CR finite element methods:
\begin{align}
\displaystyle
\bar{\sigma}_h^{RT}|_T &= \nabla \bar{u}_h^{CR} - \frac{f_T^0}{2} (x - x_T) \quad \text{on $T$}, \label{intro3}
\end{align}
where $T$ is a mesh element, $x_i$ ($i=1,2,3$) denotes vertices of the triangle $T$, $x_T$ is the barycentre of $T$ such that $x_T := \frac{1}{3}(x_1 + x_2 + x_3)$, and $\bar{\sigma}_h^{RT}$ and $ \bar{u}_h^{CR}$ respectively denote the RT and CR finite element solutions with a piecewise constant given external function $f_T^0$. It was recently proved \cite{HuMa15} that the enriched piecewise linear CR finite element method is identical to the first-order RT finite element method for the Poisson problem and Stokes problem in any number of dimensions. In the present paper, we extend Marini's results to the three-dimensional case. To this end, we present an error estimate for the first order RT finite element approximation of the Poisson problem \eqref{intro1} based on the dual mixed formulation. We again emphasise that we do not impose either the shape-regular condition or the maximum-angle condition for the mesh partition.

The remainder of the present paper is organised as follows. Section 2 introduces the weak form of the continuous problem \eqref{intro1}, finite element meshes and finite element spaces. Furthermore, we propose a parameter $H$. Section 3 introduces discrete settings of the CR finite element method for \eqref{intro1} and proposes error estimates. Section 4 proves error estimates for the first-order RT finite element method based on a dual mixed formulation of the Poisson problem. Section 5 gives the equivalence of the RT and CR finite element problems. Finally, Section 6 presents numerical results obtained using the Lagrange P1 element and the first-order CR element.
\end{comment}

\section{Preliminaries}
\subsection{Weak formulation}
The variational formulation for the Poisson problem \eqref{intro1} is then as follows. Find $u \in H_0^1(\Omega)$ such that
\begin{align}
\displaystyle
a_0(u,\varphi) = (f , \varphi) \quad \forall \varphi \in H_0^1(\Omega), \label{pre1}
\end{align}
where $a_0: H^1(\Omega) \times H^1(\Omega) \to \mathbb{R}$ denotes a bilinear form defined by
\begin{align*}
\displaystyle
a_0(u,\varphi) := (\nabla u , \nabla \varphi).
\end{align*}
Here, we define  $H_0^1(\Omega)$ as the closure of $C_0^{\infty}(\Omega)$ in the semi-norm  $| \cdot |_{H^1(\Omega)}$. By the Lax--Milgram lemma, there exists a unique solution $u \in H_0^1(\Omega)$ for any $f \in L^2(\Omega)$ and it holds that
\begin{align*}
\displaystyle
| u |_{H^1(\Omega)} \leq C_P(\Omega) \| f \|,
\end{align*}
where $C_P(\Omega)$ is the Poincar$\rm{\acute{e}}$ constant depending on $\Omega$.
Furthermore, if $\Omega$ is convex, then $u \in H^2(\Omega)$ and
\begin{align}
\displaystyle
| u |_{H^2(\Omega)} \leq \| \varDelta u \|. \label{pre2}
\end{align}
The proof can be found in, for example, \cite[Theorem 3.1.1.2, Theorem 3.2.1.2]{Gri11}.

\subsection{Meshes, Mesh faces, Averages and Jumps} \label{regularmesh}
Let $\mathbb{T}_h = \{ T \}$ be a simplicial mesh of $\overline{\Omega}$, made up of closed $d$-simplices, such as
\begin{align*}
\displaystyle
\overline{\Omega} = \bigcup_{T \in \mathbb{T}_h} T,
\end{align*}
with $h := \max_{T \in \mathbb{T}_h} h_{T}$, where $ h_{T} := \mathop{\mathrm{diam}}(T)$. We assume that each face of any $d$-simplex $T_1$ in $\mathbb{T}_h$ is either a subset of the boundary $\partial \Omega$ or a face of another $d$-simplex $T_2$ in  $\mathbb{T}_h$. That is, $\mathbb{T}_h$ is a simplicial mesh of $\overline{\Omega}$ without hanging nodes.

\begin{defi} \label{pre=defi1}
For any $T \in \mathbb{T}_h$, we define the parameter $H_{T}$ as
\begin{align*}
\displaystyle
H_{T} := \frac{h_{T}^2}{|T|} \min_{1 \leq i \leq 3} |L_i|  \quad \text{if $d=2$},
\end{align*}
where $L_i$ $(i=1,2,3)$ denotes edges of the triangle $T$. Further, we define the parameter $H_{T}$ as
\begin{align*}
\displaystyle
H_{T} := \frac{h_{T}^2}{|T|} \min_{1 \leq i , j \leq 6, i \neq j} |L_i| |L_j| \quad \text{if $d=3$},
\end{align*}
where $L_i$ $(i=1,\ldots,6)$ denotes edges of the tetrahedra $T$. Here, $|T|$ denotes the measure of $T$. Furthermore, we set
\begin{align*}
\displaystyle
H := H(h) := \max_{T \in \mathbb{T}_h} H_{T}.
\end{align*}
\end{defi}
We impose the following assumption.
 \begin{ass}  \label{pre=ass1}
 We assume that $\{ \mathbb{T}_h \}_{h \> 0}$ is a sequence of  triangulations of $\Omega$ such that
\begin{align*}
\displaystyle
\lim_{h \to 0} H(h) = 0.
\end{align*}
 \end{ass}

 \begin{rem}
 The parameter $H_T$ was introduced, and the interpolation errors are bounded (locally) in terms of $h_T$ and $H_T$ on anisotropic meshes without any geometric conditions in \cite{IshKobTsu}. In two-dimensional case, the parameter $H_T$ is equivalent to the circumradius of $T$.  Hence, the maximum-angle condition or the semiregular condition holds if and only if there exists a constant $\sigma_0$ such that $H_T / h_T \leq \sigma_0$. In three-dimensional case, it is conjectured that the maximum-angle condition holds if and only if the quantity $H_T / h_T$ is bounded.
 \end{rem}

We adopt the concepts of mesh faces, averages and jumps in the analysis of RT and CR finite element method. Let $\mathcal{F}_h^i$ be the set of interior faces and $\mathcal{F}_h^{\partial}$ the set of the faces on the boundary $\partial \Omega$. Let $\mathcal{F}_h := \mathcal{F}_h^i \cup \mathcal{F}_h^{\partial}$. For any $F \in \mathcal{F}_h$, we define the unit normal $n_F$ to $F$ as follows: (\roman{sone}) If  $F \in \mathcal{F}_h^i$ with $F = T_1 \cap T_2$, $T_1,T_2 \in \mathbb{T}_h$, let $n_1$ and $n_2$ be the outward unit normals of $T_1$ and $T_2$, respectively. Then, $n_F$ is either of $\{ n_1 , n_2\}$; (\roman{stwo}) If $F \in \mathcal{F}_h^{\partial}$, $n_F$ is the unit outward normal $n$ to $\partial \Omega$.

Let $k$ be a positive integer. We then define the broken (piecewise) Sobolev space as
\begin{align*}
\displaystyle
H^k(\mathbb{T}_h) &:= \left\{ \varphi \in L^2(\Omega); \ \varphi|_T \in H^k(T) \ \forall T \in \mathbb{T}_h  \right\}
\end{align*}
with the norm
\begin{align*}
\displaystyle
| \varphi |_{H^1(\mathbb{T}_h)} &:= \left( \sum_{T \in \mathbb{T}_h}\| \nabla \varphi \|^2_{L^2(T)^d} \right)^{1/2} \quad \varphi \in H^1(\mathbb{T}_h).
\end{align*}

Let $\varphi \in H^k(\mathbb{T}_h)$. Suppose that $F \in \mathcal{F}_h^i$ with $F = T_1 \cap T_2$, $T_1,T_2 \in \mathbb{T}_h$. Set $\varphi_1 := \varphi{|_{T_1}}$ and $\varphi_2 := \varphi{|_{T_2}}$. The jump and the average of $\varphi$ across $F$ is then defined as
\begin{align*}
\displaystyle
 [[ \varphi ]]_F := (\varphi_1 n_1 + \varphi_2 n_2) \cdot n_F, \quad \{ \! \{ \varphi \} \! \}_F := \frac{1}{2} (\varphi_1 + \varphi_2).
\end{align*}
For a boundary face $F \in \mathcal{F}_h^{\partial}$ with $F = \partial T \cap \partial \Omega$, $[[\varphi]]_F := \varphi|_T$ and $\{ \! \{ \varphi \} \! \}_F := \varphi |_T$. When $v$ is an $\mathbb{R}^d$-valued function, we use the notation
\begin{align*}
\displaystyle
[[ v \cdot n ]]_F := (v_1 - v_2) \cdot n_F, \quad \{ \! \{ v\} \! \}_F := \frac{1}{2} (v_1  + v_2 )
\end{align*}
for the jump of the normal component of $v$. For a boundary face $F \in \mathcal{F}_h^{\partial}$ with $F = \partial T \cap \partial \Omega$, $[[v \cdot n]]_F := v|_T \cdot n $ and $\{\! \{ v  \}\! \}_F := v|_T $. Whenever no confusion can arise, we simply write $[[\varphi]]$, $\{\! \{ \varphi \}\! \}$, $[[v \cdot n]]$ and $\{\! \{ v \} \! \}$, respectively.

 Suppose that $F \in \mathcal{F}_h^i$ with $F = T_1 \cap T_2$, $T_1,T_2 \in \mathbb{T}_h$. For $v \in H^1(\mathbb{T}_h)^d$ and $\varphi \in H^1(\mathbb{T}_h)$, it holds that
\begin{align*}
\displaystyle
[[ (v \varphi) \cdot n ]]_F = \{ \! \{ v \} \! \}_F \cdot n_F [[ \varphi ]]_F + [[ v \cdot n]]_F  \{ \! \{ \varphi \} \! \}_F.
\end{align*}

We here define a broken gradient operator as follows.
\begin{defi} \label{pre=defi2}
For $\varphi \in H^1(\mathbb{T}_h)$, the broken gradient $\nabla_h:H^1(\mathbb{T}_h) \to L^2(\Omega)^d$ is defined by
\begin{align*}
\displaystyle
(\nabla_h \varphi)|_T &:= \nabla (\varphi|_T) \quad \forall T \in \mathbb{T}_h.
\end{align*}
\end{defi}
\noindent
Note that $H^1(\Omega )\subset H^1(\mathbb{T}_h)$ and the broken gradient coincides with the distributional gradient in $H^1(\Omega)$.

\subsection{Finite Element Spaces and Interpolations Error Estimates} \label{FEspaces}
This section introduce the piecewise-constant, CR and RT finite element spaces.

Let $T \in \mathbb{T}_h$. For any $k \in \mathbb{N}_0$, let $\mathcal{P}^k(T)$ be the space of polynomials with degree at most $k$ in $T$.

\begin{thr}[Poincar\'e inequality] \label{modify=thr1}
Let $D \subset \mathbb{R}^d$  be a convex domain with diameter $\mathop{\mathrm{diam}}(D)$. It then holds that, for $\varphi \in H^1(D)$ with $\int_D \varphi dx = 0$,
\begin{align}
\displaystyle
\| \varphi \|_{L^2(D)} \leq \frac{\mathop{\mathrm{diam}} (D)}{\pi} |\varphi|_{H^1(D)}.\label{modify3}
\end{align}
\end{thr}

\begin{pf*}
The proof is found in \cite[Theorem 3.2]{Mar03}, also see \cite{PayWei60}.
\qed
\end{pf*}

\subsubsection{Piecewise-constant finite element space}
We define the standard piecewise constant space as
\begin{align*}
\displaystyle
M_h^0 := \left\{ q_h \in L^2(\Omega); \ q_h|_T \in \mathcal{P}^0(T) \ \forall T \in \mathbb{T}_h \right\}.
\end{align*}
The local interpolation $\Pi_{T}^0$ from $L^2(T)$ into the space $\mathcal{P}^0(T)$ is defined by
\begin{align*}
\displaystyle
\int_{T} (\Pi_{T}^0 q - q) dx = 0 \quad \forall q \in L^2(T).
\end{align*}
Note that $\Pi_{T}^0 q$ is the constant function equal to $\frac{1}{|T|} \int_T q dx$. We also define the global interpolation $\Pi_h^0$ to the space $M_{h}^0$ by
\begin{align*}
\displaystyle
(\Pi_h^0 q)|_T = \Pi_T^0 (q|_T) \quad \forall T \in \mathbb{T}_h, \quad \forall q \in L^2(\Omega).
\end{align*}

The Poincar\'e inequality \eqref{modify3} directly yields the following error estimate of the local $L^2$-projection $\Pi^0_T$.

\begin{thr}  \label{modify=thr2}
We have the error estimate of the local $L^2$-projection such that
\begin{align}
\displaystyle
\| \Pi^0_T q - q \|_{L^2(T)} &\leq \frac{h_T}{\pi} |q|_{H^1(T)} \quad \forall T \in \mathbb{T}_h, \quad \forall q \in H^1(T). \label{modify6}
\end{align}
\end{thr}

\begin{pf*}
For any $q \in H^1(T)$, we set $w := \Pi^0_T q - q$. It then holds that
\begin{align*}
\displaystyle
\int_T w dx = \int_T ( \Pi^0_T q - q ) dx = \frac{1}{|T|} \int_T q dx |T| - \int_T q dx = 0.
\end{align*}
Therefore, using the Poincar\'e inequality \eqref{modify3}, we conclude \eqref{modify6}.
\qed
\end{pf*}

The global error estimate of the $L^2$-projection is obtained as follows.
\begin{coro}  \label{modify=coro1}
Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. It then holds that
\begin{align}
\displaystyle
\| \Pi^0_h q - q \|_{L^p(\Omega)} &\leq \frac{h}{\pi} |q|_{H^{1}(\Omega)} \quad \forall q \in H^{1}(\Omega). \label{modify7}
\end{align}
\end{coro}

\subsubsection{CR finite element space}
We define the following CR finite element space as
\begin{align*}
\displaystyle
CR_{h0}^1 &:= \left \{ \varphi_h \in L^2(\Omega); \ \varphi_h|_T \in \mathcal{P}^1(T) \ \forall T \in \mathbb{T}_h, \ \int_F [[ \varphi_h ]]_F ds = 0 \ \forall F \in \mathcal{F}_h \right \}.
\end{align*}
Using the barycentric coordinates $ \lambda_i: \mathbb{R}^d \to \mathbb{R}$, $i=1,\ldots,d+1$, we define the local basis functions as
\begin{align*}
\displaystyle
\theta_i(x) := d \left( \frac{1}{d} - \lambda_i(x) \right), \quad 1 \leq i \leq d+1.
\end{align*}
For $i=1,\ldots,d+1$, let $F_i$ be the face of $T$ and $x_{F_i}$ the barycentre of the face $F_i$. We then define the local CR interpolation operator as
\begin{align*}
\displaystyle
I_{T}^{CR}: H^{1}(T) \ni \varphi  \mapsto I_{T}^{CR} \varphi := \sum_{i=1}^{d+1} \left( \frac{1}{|F_i|} \int_{F_i} \varphi ds \right)  \theta_i \in \mathcal{P}^1.
\end{align*}
Furthermore, it holds that
\begin{align*}
\displaystyle
 \frac{1}{|F_i|} \int_{F_i} \left( I_{T}^{CR} \varphi -  \varphi \right) ds = 0, \quad i=1,\ldots,d+1, \quad \forall \varphi \in H^1(T).
\end{align*}
We define the global CR interpolation $I_h^{CR}:H_0^{1}(\Omega) \to CR_{h0}^1$ by
\begin{align*}
\displaystyle
(I_h^{CR} \varphi )|_T = I_T^{CR} (\varphi|_T) \quad \forall T \in \mathbb{T}_h, \quad \forall \varphi \in H_0^{1}(\Omega).
\end{align*}

We give the local CR interpolation error estimate.

\begin{thr}  \label{modify=thr3}
We have the following estimates such that for $m \in \{ 0 , 1 \}$,
\begin{align}
\displaystyle
 |\varphi - {I}^{CR}_{T} \varphi|_{H^{m}(T)}
&\leq C_I^{CR,m} \left( \frac{H_{T}}{h_{T}} \right)^m h_{T}^{2-m} | \varphi |_{H^{2}(T)} \ \forall T \in \mathbb{T}_h, \ \forall \varphi \in H^{2}(T). \label{modify10}
\end{align}
Here,  $C_I^{CR,m}$ is a positive constant independent of $h_T$ and $H_T$.
\end{thr}

\begin{pf*}
The proof is found in \cite[Theorem 2]{IshKobTsu}.
\qed
\end{pf*}

The global CR interpolation error estimates are obtained as follows.

\begin{coro}  \label{modify=coro2}
Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Then, there exists constants $C_G^{CR,0},  C_G^{CR,1}\> 0$, independent of $H$ and $h$, such that
\begin{align}
\displaystyle
 \|\varphi - {I}^{CR}_{h} \varphi \|_{L^{2}(\Omega)}
&\leq C_G^{CR,0} h^2 | \varphi |_{H^{2}(\Omega)} \quad  \forall \varphi \in H_0^1(\Omega) \cap H^{2}(\Omega),\label{modify11} \\
 |\varphi - {I}^{CR}_{h} \varphi|_{H^{1}(\mathbb{T}_h)}
&\leq C_G^{CR,1} H | \varphi |_{H^{2}(\Omega)} \quad  \forall \varphi \in H_0^1(\Omega) \cap H^{2}(\Omega). \label{modify12}
\end{align}
\end{coro}

The inequality \eqref{modify10} with $m=1$ can be improved by replacing $H_{T}$ with $h_{T}$. To this end, we use the Poincar\'e inequality \eqref{modify3}.

\begin{thr}  \label{modify=thr4}
It then holds that
\begin{align}
\displaystyle
| I_{T}^{CR} \varphi - \varphi |_{H^1(T)}  &\leq \frac{h_{T}}{\pi} |\varphi|_{H^2(T)} \quad \forall T \in \mathbb{T}_h, \quad \forall \varphi \in H^2(T). \label{modify13}
\end{align}
\end{thr}

\begin{pf*}
Let $F_i$, $i = 1,\ldots,d+1$ be the faces of the element $T$.
We set $\psi := I_{T}^{CR} \varphi - \varphi \in H^2(T)$. From Green's formula and the property of the CR interpolation, we have
\begin{align*}
\displaystyle
\int_{T} \frac{\partial \psi}{\partial x_j} dx = \int_{\partial T} \psi n_{T}^{(j)} ds = \sum_{i=1}^{d+1} n_{T}^{(j)} \int_{F_i} \psi ds = 0,
\end{align*}
where $ n_{T}^{(j)}$ denotes the $j$th component of the outer unit normal vector $n_{T}$. From  the Poincar\'e inequality \eqref{modify3}, we have
\begin{align*}
\displaystyle
| I_{T}^{CR} \varphi - \varphi |_{H^1(T)}^2
&= \sum_{j=1}^d \left \|  \frac{\partial}{\partial x_j} ( I_{T}^{CR} \varphi - \varphi)  \right\|^2_{L^2(T)} \\
&\leq \left( \frac{h_{T}}{\pi} \right)^2 \sum_{j=1}^d \left | \frac{\partial}{\partial x_j} ( I_{T}^{CR} \varphi - \varphi) \right |^2_{H^1(T)} \\
&=  \left( \frac{h_{T}}{\pi} \right)^2 \sum_{j,k=1}^d  \left \|  \frac{\partial^2}{\partial x_j \partial x_k} ( I_{T}^{CR} \varphi - \varphi) \right \|^2_{L^2(T)} \\
&= \left( \frac{h_{T}}{\pi} \right)^2 |\varphi|^2_{H^2(T)},
\end{align*}
which conclude \eqref{modify13}.
\qed
\end{pf*}

\begin{rem} \label{modify=rem1}
For $i=1,\ldots,d+1$, let $x_{F_i}$ the barycentre of face $F_i$. If we choose the domain of the local CR interpolation operator as $W^{\ell,p}(T) \subset \mathcal{C}^0(T)$ with $1 \leq p \< \infty$ and $d \< \ell p$, it is possible to define
\begin{align*}
\displaystyle
I_{T}^{CR,S}:W^{\ell,p}(T) \ni \varphi  \mapsto I_{T}^{CR,S} \varphi := \sum_{i=1}^{d+1} \varphi (x_{F_i}) \theta_i \in \mathcal{P}^1.
\end{align*}
However, the estimate \cite[Theorem 2]{IshKobTsu}
\begin{align*}
\displaystyle
| I_{T}^{CR,S} \varphi - \varphi |_{H^{1}(T)}  &\leq C_I^{CR,1} H_{T} |\varphi|_{H^{2}(T)}  \quad \forall \varphi \in H^{2}(T)
\end{align*}
can not be improved by replacing $H_{T}$ with $h_{T}$.

As a counter example, let us consider $T$ with vertices $x_1 := (0,0,0)^T$, $x_2 := (h,0,0)^T$, $x_3 := (\frac{h}{2}, h^{\gamma},0)^T$ and  $x_4 := (\frac{h}{2}, 0, \frac{h}{2})^T$, where $h := \frac{1}{N}$, $N \in \mathbb{N}$ and $\gamma \in \mathbb{R}$, $1 \< \gamma \leq 2$. Let $\varphi$ be a function such that
\begin{align*}
\displaystyle
\varphi(x,y,z) := x^2 + y^2 + z ^2.
\end{align*}
If an exact solution $\varphi$ is known, the error $e_h := \varphi - \varphi_h$ and $e_{h/2} := \varphi - \varphi_{h/2}$ are computed numerically for two mesh sizes $h$ and $h/2$, where $\varphi_h := I_{T}^{CR,S} \varphi$. The convergence indicator $r$ is defined by
\begin{align*}
\displaystyle
r = \frac{1}{\log(2)} \log \left( \frac{\| e_h \|_X}{\| e_{h/2} \|_X} \right).
\end{align*}
The parameter $H_{T}$ is then $H_{T} = \mathcal{O}( h^{2 - \gamma} )$. We compute the convergence order with respect to the $H_0^1$ norm defined by
\begin{align*}
\displaystyle
&Err_h^{CR,S}(H^1) := \frac{| \varphi - I_{T}^{CR,S} \varphi |_{H^1(T)}}{| \varphi |_{H^2(T)}},
\end{align*}
for the case: $\gamma = 1.5$ (Table \ref{modify=table1}).

\begin{table}[htb]
\caption{Error of the local CR interpolation operator ($\gamma = 1.5$)}
\centering
\begin{tabular}{l | l | l | l | l } \hline
$N$ &  $h$ & $H_{T}$ & $Err_h^{CR,S}(H^1)$ & $r$  \\ \hline \hline
128 & 7.8125e-03 & 3.8081e-01  & 2.8183e-03  &    \\
256 & 3.9062e-03  & 2.6723e-01  &  1.7641e-03 & 0.68  \\
512 & 1.9531e-03  & 1.8823e-01  & 1.1587e-03   &  0.61  \\
1024 & 9.7656e-04  & 1.3284e-01  &  7.8625e-04  & 0.60  \\
2048 & 4.8828e-04 & 9.3842e-02  &  5.4390e-04 &  0.53 \\
4096 & 2.4414e-04  &  6.6324e-02 &  3.8026e-04  & 0.52 \\
\hline
\end{tabular}
\label{modify=table1}
\end{table}
\end{rem}

\subsubsection{RT finite element space}
The lowest order RT finite element space is defined by
\begin{align*}
\displaystyle
RT^0(T) := \{ v ; \ v(x) = p + x q, \ p \in \mathcal{P}^0(T)^d, \ q \in \mathcal{P}^0(T), \ x \in \mathbb{R}^d \}.
\end{align*}
The functionals are defined by, for any $v \in RT^0(T)$,
 \begin{align*}
\displaystyle
{\chi}_{i} (v) := \frac{1}{|F_i|} \int_{F_i} v \cdot n_{i} ds, \quad F_i \subset \partial T, \quad 1 \leq i \leq d+1,
\end{align*}
where $n_i$ denotes the outer unit normal vector of $T$ along $F_i$. We set $\sum := \{ {\chi}_{i} \}_{i=1}^{d+1}$. Note that $\dim RT^0(T) = d+1$. The triple $\{ T , RT^0(T) , \Sigma \}$ is then a finite element. We define the RT finite element space by
\begin{align*}
\displaystyle
RT^0_{h} &:= \{ v_h \in L^2(\Omega)^d;  \ v_h |_T \in RT^0(T), \ \forall T \in \mathbb{T}_h, \ [[ v_h \cdot n ]]_F = 0, \ \forall F \in \mathcal{F}_h^i \}.
\end{align*}
Note that $RT^0_{h} \subset H(\mathop{\mathrm{div}};\Omega)  := \left\{ v \in L^2(\Omega)^d ; \ \mathop{\mathrm{div}} v \in L^2(\Omega) \right\}$.

We next define the local RT interpolation as
\begin{align*}
\displaystyle
I_T^{RT}: H^1(T)^d \to RT^0(T),
\end{align*}
using
\begin{align*}
\displaystyle
 \int_{F_i} (v - I_T^{RT} v) \cdot n_i ds = 0, \quad F_i \subset \partial T, \ i \in \{ 1 , \ldots, d+1\} \quad \forall v \in H^1(T)^d.
\end{align*}
Further, we define the global RT interpolation $I_h^{RT} : H^1(\Omega)^d \to RT^0_{h}$ by
\begin{align*}
\displaystyle
(I_h^{RT} v )|_T = I_T^{RT} (v|_T) \quad \forall T \in \mathbb{T}_h, \quad \forall v \in H^1(\Omega)^d.
\end{align*}

The local RT interpolation error estimate is as follows.

\begin{thr}  \label{modify=thr5}
We have the following estimates such that
\begin{align}
\displaystyle
\| I_{T}^{RT} v - v \|_{L^2(T)^d}
&\leq  C_I^{RT} H_{T} |v|_{H^{1}(T)^d} \quad \forall T \in \mathbb{T}_h, \quad  \forall v \in H^{1}(T)^d. \label{modify14}
\end{align}
Here,  $C_I^{RT}$ is a positive constant independent of $H_T$.
\end{thr}

\begin{pf*}
The proof is found in \cite[Theorem 3]{IshKobTsu}.
\qed
\end{pf*}

The global RT interpolation error estimates are obtained as follows.

\begin{coro}  \label{modify=coro3}
Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Then, there exists a constant $C_G^{RT} \> 0$, independent of $H$, such that
\begin{align}
\displaystyle
\| I_{h}^{RT} v - v \|_{L^2(\Omega)^d}
&\leq  C_G^{RT} H |v|_{H^{1}(\Omega)^d} \quad  \forall v \in H^{1}(\Omega)^d. \label{modify15}
\end{align}
\end{coro}

\begin{comment}

\subsubsection{RT finite element space}
Let $T \in \mathbb{T}_h$. For any $k \in \mathbb{N}_0$, let $\mathcal{P}^k(T)$ be the space of polynomials with degree at most $k$ in $T$.

The lowest order RT finite element space is defined by
\begin{align*}
\displaystyle
RT^0(T) := \{ v ; \ v(x) = p + x q, \ p \in \mathcal{P}^0(T)^d, \ q \in \mathcal{P}^0(T), \ x \in \mathbb{R}^d \}.
\end{align*}
The functionals are defined by, for any $v \in RT^0(T)$,
 \begin{align}
\displaystyle
{\chi}_{i} (v) := \frac{1}{|F_i|} \int_{F_i} v \cdot n_{i} ds, \quad F_i \subset \partial T, \quad 1 \leq i \leq d+1, \label{pre3}
\end{align}
where $n_i$ denotes the outer unit normal vector of $T$ along $F_i$. We set $\sum := \{ {\chi}_{i} \}_{i=1}^{d+1}$. Note that $\dim RT^0(T) = d+1$. The triple $\{ T , RT^0(T) , \Sigma \}$ is then a finite element. We define the global RT finite element space by
\begin{align*}
\displaystyle
RT^0_{h} &:= \{ v_h \in L^2(\Omega)^d;  \ v_h |_T \in RT^0(T), \ \forall T \in \mathbb{T}_h, \ [[ v_h \cdot n ]]_F = 0, \ \forall F \in \mathcal{F}_h^i \}.
\end{align*}
Note that $RT^0_{h} \subset H(\mathop{\mathrm{div}};\Omega)  := \left\{ v \in L^2(\Omega)^d ; \ \mathop{\mathrm{div}} v \in L^2(\Omega) \right\}$.

We next define the local RT interpolation as
\begin{align}
\displaystyle
I_T^{RT}: H^1(T)^d \to RT^0(T), \label{pre4}
\end{align}
using
\begin{align}
\displaystyle
 \int_{F_i} (v - I_T^{RT} v) \cdot n_i ds = 0, \quad F_i \subset \partial T, \ i \in \{ 1 , \ldots, d+1\} \quad \forall v \in H^1(T)^d. \label{pre5}
\end{align}
Further, we define the global RT interpolation $I_h^{RT} : H(\mathop{\mathrm{div}};\Omega) \cap H^1(\mathbb{T}_h)^d \to RT^0_{h}$ by
\begin{align}
\displaystyle
(I_h^{RT} v )|_T = I_T^{RT} (v|_T) \quad \forall T \in \mathbb{T}_h, \quad \forall v \in H^1(\Omega)^d. \label{pre6}
\end{align}

We give the local RT interpolation error estimate.
\begin{thr}  \label{pre=thr1}
We have the following estimate such that
\begin{align}
\displaystyle
\| I_T^{RT} v - v \|_{L^2(T)^d} &\leq  C_I^{RT} H_T |v|_{H^1(T)^d} \quad  \forall T \in \mathbb{T}_h, \quad \forall v \in H^1(T)^d, \label{pre7}
\end{align}
where $C_I^{RT}$ is a positive constant independent of $H_T$.
\end{thr}

\begin{pf*}
The proof can be found in \cite[Theorem 3]{IshKobTsu}.
\qed
\end{pf*}

The global RT interpolation error estimate is obtained as follows.
\begin{coro}  \label{pre=coro1}
Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Then, there exists $C_{g}^{RT} \> 0$, independent of $H$, such that
\begin{align}
\displaystyle
\| I_h^{RT} v - v \|_{L^2(\Omega)^d} &\leq  C_{g}^{RT} H |v|_{H^1(\Omega)^d} \quad \forall v \in H^1(\Omega)^d. \label{pre8}
\end{align}
\end{coro}

\subsubsection{CR finite element space}
In introducing a nonconforming method, we define the following CR finite element space as
\begin{align*}
\displaystyle
CR_{h0}^1 &:= \left \{ \varphi_h \in L^2(\Omega); \ \varphi_h|_T \in \mathcal{P}^1(T) \ \forall T \in \mathbb{T}_h, \ \int_F [[ \varphi_h ]]_F ds = 0 \ \forall F \in \mathcal{F}_h \right \}.
\end{align*}
Using the barycentric coordinates $ \lambda_i: \mathbb{R}^d \to \mathbb{R}$, $i=1,\ldots,d+1$, we define the local basis functions as
\begin{align}
\displaystyle
\theta_i(x) := d \left( \frac{1}{d} - \lambda_i(x) \right), \quad 1 \leq i \leq d+1. \label{pre9}
\end{align}

For $i=1,\ldots,d+1$, let $F_i$ be the face of $T$ and $x_{F_i}$ the barycentre of face $F_i$. We then define the local CR interpolation operator as
\begin{align}
\displaystyle
I_T^{CR}: W^{1,1}(T) \ni \varphi  \mapsto I_T^{CR} \varphi := \sum_{i=1}^{d+1} \left( \frac{1}{|F_i|} \int_{F_i} \varphi ds \right)  \theta_i \in \mathcal{P}^1. \label{pre10}
\end{align}
Because the trace of a function in $W^{1,1}(T)$ is in $L^1(\partial T)$, $ \frac{1}{|F_i|} \int_{F_i} \varphi ds$ is meaningful.
Further, it holds that
\begin{align}
\displaystyle
 \frac{1}{|F_i|} \int_{F_i} \left( I_T^{CR} \varphi -  \varphi \right) ds = 0, \quad i=1,\ldots,d+1. \label{pre11}
\end{align}
We define the global CR interpolation $I_h^{CR}:W^{1,1}(\Omega) \to CR_{h0}^1$ by
\begin{align}
\displaystyle
(I_h^{CR} \varphi )|_T = I_T^{CR} (\varphi|_T) \quad \forall T \in \mathbb{T}_h, \quad \forall \varphi \in W^{1,1}(\Omega). \label{pre12}
\end{align}

 We give the local CR interpolation error estimate.
\begin{thr}  \label{pre=thr2}
We have the following estimates such that
\begin{align}
\displaystyle
\| I_T^{CR} \varphi - \varphi \|_{L^2(T)} &\leq  C_{I}^{CR,L^2} h_T^2 |\varphi|_{H^2(T)} \quad  \forall T \in \mathbb{T}_h, \quad \forall \varphi \in H^2(T), \label{pre13}\\
| I_T^{CR} \varphi - \varphi |_{H^1(T)}  &\leq C_I^{CR,H^1} H_T    |\varphi|_{H^2(T)}  \quad  \forall T \in \mathbb{T}_h, \quad \forall \varphi \in H^2(T). \label{pre14}
\end{align}
Here,  $C_I^{CR,L^2}$ and $C_I^{CR,H^1}$ are positive constants independent of $h_T$ and $H_T$.
\end{thr}

\begin{pf*}
The proof can be found in \cite[Theorem 2]{IshKobTsu}.
\qed
\end{pf*}

\begin{rem*}
The inequality \eqref{pre14} can be improved by replacing $H_T$
with $h_T$.
\end{rem*}

The global CR interpolation error estimates are obtained as follows.
\begin{coro}  \label{pre=coro2}
Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Then, there exist $ C_{g}^{CR,L^2},  C_{g}^{CR,H^1} \> 0$, independent of $H$ and $h$, such that
\begin{align}
\displaystyle
\| I_h^{CR} \varphi - \varphi \| &\leq  C_{g}^{CR,L^2} h^2 |\varphi|_{H^2(\Omega)} \quad \forall \varphi \in H^2(\Omega), \label{pre15}\\
| I_h^{CR} \varphi - \varphi |_{H^1(\mathbb{T}_h)}  &\leq C_g^{CR,H^1} H |\varphi|_{H^2(\Omega)}  \quad \forall \varphi \in H^2(\Omega). \label{pre16}
\end{align}
\end{coro}

\subsubsection{Piecewise-constant finite element space}
We define the standard piecewise constant space as
\begin{align*}
\displaystyle
M_h^0 := \left\{ q_h \in L^2(\Omega); \ q_h|_T \in \mathcal{P}^0(T) \ \forall T \in \mathbb{T}_h \right\}.
\end{align*}
The local $L^2$-projection $\Pi_{T}^0$ from $L^2(T)$ into the space $\mathcal{P}^0(T)$ is defined by
\begin{align}
\displaystyle
\int_{T} (\Pi_{T}^0 q - q) dx = 0 \quad \forall q \in L^2(T).  \label{pre17}
\end{align}
Note that $\Pi_{T}^0 q$ is the constant function equal to $\frac{1}{|T|} \int_T q dx$. We also define the global $L^2$-projection $\Pi_h^0$ to the space $M_{h}^0$ by
\begin{align}
\displaystyle
(\Pi_h^0 q)|_T = \Pi_T^0 (q|_T) \quad \forall T \in \mathbb{T}_h, \quad \forall q \in L^2(\Omega). \label{pre18}
\end{align}
The error estimate of the $L^2$-projection is as follows.
\begin{thr}  \label{pre=thr3}
We have the error estimate of the $L^2$-projection such that
\begin{align}
\displaystyle
\| \Pi^0_T q - q \|_{L^2(T)} &\leq C_P^{L^2} h_T  |q|_{H^1(T)} \quad \forall T \in \mathbb{T}_h, \quad \forall q \in H^1(T). \label{pre19}
\end{align}
Here, $C_P^{L^2}$ is a positive constant independent of $h_T$.
\end{thr}

\begin{pf*}
The proof can be found in \cite[Theorem 2]{IshKobTsu}.
\qed
\end{pf*}

The global error estimate of the $L^2$-projection is obtained as follows.
\begin{coro}  \label{pre=coro3}
Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Then, there exists $C_{Pg}^{L^2} \> 0$, independent of $h$, such that
\begin{align}
\displaystyle
\| \Pi^0_h q - q \| &\leq  C_{Pg}^{L^2} h  |q|_{H^1(\Omega)} \quad \forall q \in H^1(\Omega). \label{pre20}
\end{align}
\end{coro}

\end{comment}

Between the RT interpolation $I_h^{RT}$ and the $L^2$-projection $\Pi_h^0$, the following relation holds:
\begin{lem} \label{pre=lem1}
For any $v \in H^1(\Omega)^d$, it holds that
\begin{align*}
\displaystyle
\mathop{\mathrm{div}} (I_h^{RT} v) = \Pi_h^0 (\mathop{\mathrm{div}} v).
\end{align*}
That is to say, the diagram
\[
\begin{CD}
     H^1(\Omega)^d @>{\mathop{\mathrm{div}}}>> L^2(\Omega) \\
  @V{I_h^{RT}}VV    @VV{\Pi_h^0}V \\
     RT^0_{h}   @>{\mathop{\mathrm{div}}}>>  M_{h}^0
\end{CD}
 \]
commutes.
\end{lem}

\begin{pf*}
The proof of this lemma is found in \cite{Bra07}.
\qed
\end{pf*}

The following relation plays an important role in the CR finite element analysis on anisotropic meshes.
\begin{lem} \label{pre=lem2}
It holds that
\begin{align}
\displaystyle
(v_h , \nabla_h \psi_h) + (\mathop{\mathrm{div}} v_h , \psi_h) = 0 \quad \forall v_h \in RT_h^0, \quad \forall \psi_h \in H_0^1(\Omega) +  CR_{h0}^1.  \label{pre21}
\end{align}
\end{lem}

\begin{pf*}
For any $v_h \in RT_h^0$ and $\psi_h \in H_0^1(\Omega) +  CR_{h0}^1$, using Green formula and the fact $v_h \cdot n_F \in \mathcal{P}^0(F)$ for any $F \in \mathcal{F}_h$, we can derive
\begin{align*}
\displaystyle
(v_h , \nabla_h \psi_h) + (\mathop{\mathrm{div}} v_h , \psi_h)
&= \sum_{T \in \mathbb{T}_h} \int_{\partial T} (v_h \cdot n_T) \psi_h ds  \notag \\
&= \sum_{F \in \mathcal{F}_h} \int_{F} [[ (v_h  \psi_h ) \cdot n_F ]] ds  \notag \\
&= \sum_{F \in \mathcal{F}_h} \int_{F} \left(  [[ v_h \cdot n_F]] \{ \! \{ \psi_h \} \!\} + \{ \! \{ v_h  \} \!\} \cdot n_F [[ \psi_h ]] \right) ds \notag \\
&=  0.
\end{align*}
\qed
\end{pf*}

\subsection{Discrete Poincar\'e Inequality on Anisotropic Meshes}
We propose the discrete Poincar\'e inequality on anisotropic meshes.
\begin{lem}[Discrete Poincar\'e inequality on anisotropic meshes] \label{pre=lem3}
Assume that $\Omega$ is convex.  If $H \leq 1$,  there exists $C(\Omega)$, independent of $h$, $H$, and the geometry of meshes, such that
\begin{align}
\displaystyle
 \| \varphi_h \| \leq C(\Omega) |\varphi_h|_{H^1(\mathbb{T}_h)} \quad \forall \varphi_h \in CR_{h0}^1. \label{pre22}
\end{align}
\end{lem}

\begin{pf*}
Let $\varphi_h \in CR_{h0}^1$. We consider the dual problem. Find $z \in H^{2}(\Omega) \cap H_0^1(\Omega)$ such that
\begin{align*}
\displaystyle
- \varDelta z = \frac{\varphi_h}{\| \varphi_h \|} \quad \text{in $\Omega$}, \quad z = 0 \quad \text{on $\partial \Omega$}.
\end{align*}
We then have a priori estimates:
\begin{align*}
\displaystyle
|z|_{H^1(\Omega)} \leq C_P, \quad |z|_{H^2(\Omega)} \leq 1,
\end{align*}
where $C_P$ is the Poincar\'e constant. We use the duality argument to show the target inequality. That is to say, we have
\begin{align*}
\displaystyle
\| \varphi_h \| &= \frac{1}{\| \varphi_h \|} (\varphi_h , \varphi_h) = (- \varDelta z , \varphi_h) = (- \mathop{\mathrm{div}} \nabla z , \varphi_h) \\
&= (- \mathop{\mathrm{div}} \nabla z , \varphi_h - \Pi_h^0 \varphi_h) - (\nabla z - I_h^{RT}(\nabla z) , \nabla_h \varphi_h) + (\nabla z , \nabla_h \varphi_h) \\
&\leq \| \varDelta z\| \|  \varphi_h - \Pi_h^0 \varphi_h \| + \| \nabla z - I_h^{RT}(\nabla z) \| |\varphi_h|_{H^1(\mathbb{T}_h)} + |z|_{H^1(\Omega)} |\varphi_h|_{H^1(\mathbb{T}_h)} \\
&\leq c \left( h+ H |\nabla z|_{H^1(\Omega)}+ C_P \right) |\varphi_h|_{H^1(\mathbb{T}_h)},
\end{align*}
which leads to
\begin{align*}
\displaystyle
\| \varphi_h \|
\leq c (2 + C_p) |\varphi_h|_{H^1(\mathbb{T}_h)} \quad \text{if $H \leq 1$}.
\end{align*}
We here used
\begin{align*}
\displaystyle
- \int_{\Omega} \mathop{\mathrm{div}} (\nabla z)  \varphi_h dx
&=  \int_{\Omega} ( \Pi_h^0  \mathop{\mathrm{div}} (\nabla z) - \mathop{\mathrm{div}} (\nabla z) ) \varphi_h dx - \int_{\Omega} ( \Pi_h^0  \mathop{\mathrm{div}} (\nabla z) ) \varphi_h dx\\
&=  \int_{\Omega} ( \Pi_h^0  \mathop{\mathrm{div}} (\nabla z) - \mathop{\mathrm{div}} (\nabla z) ) ( \varphi_h - \Pi_h^0 \varphi_h) dx \\
&\quad  - \int_{\Omega} ( \mathop{\mathrm{div}} I_h^{RT}  (\nabla z) ) \varphi_h dx\\
&= - \int_{\Omega} \mathop{\mathrm{div}} (\nabla z) \left( \varphi_h - \Pi_h^0 \varphi_h \right) dx \\
&\quad - \int_{\Omega} (\nabla z - I_h^{RT} (\nabla z) ) \cdot \nabla_h \varphi_h dx + \int_{\Omega} \nabla z \cdot \nabla_h \varphi_h dx,
\end{align*}
where
\begin{align*}
\displaystyle
 \int_{\Omega} ( \mathop{\mathrm{div}} I_h^{RT}  (\nabla z) ) \varphi_h dx
 &= \sum_{T \in \mathbb{T}_h} \int_{\partial T} n_T \cdot  I_h^{RT}  (\nabla z) \varphi_h ds - \int_{\Omega}  I_h^{RT}  (\nabla z) \cdot \nabla_h \varphi_h dx \\
 &=  \int_{\Omega} ( \nabla z -  I_h^{RT}  (\nabla z) ) \cdot \nabla_h \varphi_h dx  - \int_{\Omega}  \nabla z \cdot \nabla_h \varphi_h dx.
\end{align*}
\qed
\end{pf*}

\begin{comment}
\begin{lem}[Discrete Sobolev inequality on anisotropic meshes] \label{pre=lem3}
Assume that $\Omega$ is convex.  Let $r$ be a real number satisfying $r \in [2,\infty)$ if $d=2$ and $r \in [2,6]$ if $d=3$. If $\max_{T \in \mathbb{T}_h} ( |T|^{(2-r)/2r} H_T) \leq 1$ for any $T \in \mathbb{T}_h$, then there exists $c$, independent of $h$, for all $v_h \in CR_{h0}^1$,
\begin{align}
\displaystyle
 \| v_h \|_{L^r(\Omega)} \leq c |v_h|_{H^1(\mathbb{T}_h)}. \label{pre22}
\end{align}
\end{lem}

\begin{pf*}
Let $v_h \in CR_{h0}^1$. Let $r^{\prime}$ be a real number with $1/r + 1/r^{\prime} = 1$. We consider the dual problem. For any  $w_h \in CR_{h0}^1$, find the solution $z \in H^{2}(\Omega) \cap H_0^1(\Omega)$ such that
\begin{align*}
\displaystyle
- \varDelta z = w_h \quad \text{in $\Omega$}, \quad z = 0 \quad \text{on $\partial \Omega$},
\end{align*}
which leads to
\begin{align*}
\displaystyle
- \int_{\Omega} \varDelta z v_h dx = \langle v_h , w_h \rangle_{L^r(\Omega),L^{r^{\prime}}(\Omega)} := \int_{\Omega} v_h w_h dx.
\end{align*}
We first have the a priori estimate
\begin{align*}
\displaystyle
\| \nabla z \|_{L^2(\Omega)^d}^2
&= (\nabla z , \nabla z) = (- \varDelta z , z) = \int_{\Omega} w_h z dx
\leq \| w_h\|_{L^{r^{\prime}}(\Omega)} \| z \|_{L^r(\Omega)} \\
&\leq c  \| w_h \|_{L^{r^{\prime}}(\Omega)} \| \nabla z \|_{L^2(\Omega)^d}.
\end{align*}
We use the duality argument to show the target inequality. That is to say, we have
\begin{align*}
\displaystyle
&\| v_h \|_{L^r(\Omega)} \\
&= \sup_{w_h \in L^{r^{\prime}}(\Omega)}  \frac{\langle v_h , w \rangle_{L^r(\Omega),L^{r^{\prime}}(\Omega)}}{\| w \|_{ L^{r^{\prime}}(\Omega)}} \\
&= \sup_{w \in L^{r^{\prime}}(\Omega)} \frac{- \int_{\Omega} \varDelta z v_h dx}{\| w \|_{ L^{r^{\prime}}(\Omega)}} = \sup_{w \in L^{r^{\prime}}(\Omega)} \frac{- \int_{\Omega} \mathop{\mathrm{div}} (\nabla z)  v_h dx}{\| w \|_{ L^{r^{\prime}}(\Omega)}} \\
&\leq \sup_{w \in L^{r^{\prime}}(\Omega)} \left(  \frac{ \| \varDelta z \|_{ L^{r^{\prime}}(\Omega)} \| v_h - \Pi_h^0 v_h \|_{L^r(\Omega)} }{\| w \|_{ L^{r^{\prime}}(\Omega)}} \right.  \\
&\quad \quad \quad \quad \quad \left.+ \frac{ \| \nabla z - I_h^{RT} ( \nabla z)  \|_{L^{r^{\prime}}(\Omega)} \| \nabla _h v_h \|_{L^{r}(\Omega)} + |z|_{H^1(\Omega)} |v_h|_{H^1(\mathbb{T}_h)}}{\| w \|_{ L^{r^{\prime}}(\Omega)}} \right).
\end{align*}
Here, we used
\begin{align*}
\displaystyle
 \int_{\Omega} ( \mathop{\mathrm{div}} I_h^{RT}  (\nabla z) ) v_h d
 &= \sum_{T \in \mathbb{T}_h} \int_{\partial T} n_T \cdot  I_h^{RT}  (\nabla z) v_h ds - \int_{\Omega}  I_h^{RT}  (\nabla z) \cdot \nabla_h v_h dx \\
 &=  \int_{\Omega} ( \nabla z -  I_h^{RT}  (\nabla z) ) \cdot \nabla_h v_h dx  - \int_{\Omega}  \nabla z \cdot \nabla_h v_h dx,
\end{align*}
and
\begin{align*}
\displaystyle
- \int_{\Omega} \mathop{\mathrm{div}} (\nabla z)  v_h dx
&=  \int_{\Omega} ( \Pi_h^0  \mathop{\mathrm{div}} (\nabla z) - \mathop{\mathrm{div}} (\nabla z) ) v_h dx - \int_{\Omega} ( \Pi_h^0  \mathop{\mathrm{div}} (\nabla z) ) v_h dx\\
&=  \int_{\Omega} ( \Pi_h^0  \mathop{\mathrm{div}} (\nabla z) - \mathop{\mathrm{div}} (\nabla z) ) ( v_h - \Pi_h^0 v_h) dx \\
&\quad  - \int_{\Omega} ( \mathop{\mathrm{div}} I_h^{RT}  (\nabla z) ) v_h dx\\
&= - \int_{\Omega} \mathop{\mathrm{div}} (\nabla z) \left( v_h - \Pi_h^0 v_h \right) dx \\
&\quad + \int_{\Omega} (\nabla z - I_h^{RT} (\nabla z) ) \cdot \nabla_h v_h dx - \int_{\Omega} \nabla z \cdot \nabla_h v_h dx.
\end{align*}

Using the standard scaling argument, e.g., see \cite[Lemma 1.101]{ErnGue04}, we have
\begin{align*}
\displaystyle
 \| v_h - \Pi_T^0 v_h \|_{L^r(T)}
 &\leq c |T|^{1/r}  \| \hat{v}_h - {\Pi_{\widehat{T}}^0} \hat{v}_h \|_{L^r(\widehat{T})} \\
 &\leq c |T|^{1/r}  \| \hat{v}_h - {\Pi_{\widehat{T}}^0} \hat{v}_h \|_{L^2(\widehat{T})} \\
 &\leq c |T|^{1/r-1/2}  \| {v}_h - {\Pi_T^0} {v}_h \|_{L^2(T)} \\
 &\leq c |T|^{(2-r)/2r} h_T  \| \nabla {v}_h  \|_{L^2(T)}.
\end{align*}
By analogous argument, we have
\begin{align*}
\displaystyle
| v_h |_{W^{1,r}(T)}
&\leq c \| A^{-1} \| |T|^{1/r} | \widehat{v_h} |_{W^{1,r}(\widehat{T})}
\leq c  \| A^{-1} \| |T|^{1/r} | \widehat{v_h} |_{H^1(\widehat{T})} \\
&\leq c  \| A^{-1} \|  \| A \| |T|^{1/r - 1/2} | { v_h} |_{H^1({T})} \\
&\leq c \frac{H_T}{h_T} |T|^{1/r - 1/2} | { v_h} |_{H^1({T})}.
\end{align*}
Since $\Omega$ is convex, we obtain the global estimate, e.g. see \cite{Gri11} ???:
\begin{align*}
\displaystyle
  \| \nabla z - I_h^{RT} ( \nabla z)  \|_{L^{r^{\prime}}(\Omega)}
  &\leq c \sum_{T \in \mathbb{T}_h}  \| \nabla z - I_{T}^{RT} ( \nabla z)  \|_{L^{r^{\prime}}(T)}
  \leq cH |\nabla z|_{W^{1,r^{\prime}}(\Omega)} \\
  &\leq c H \| \varDelta z \|_{L^{r^{\prime}}(\Omega)}
  =  c H \| w \|_{L^{r^{\prime}}(\Omega)}.
\end{align*}

If $\max_{T \in \mathbb{T}_h} ( |T|^{(2-r)/2r}) H \leq 1$, it holds that
\begin{align*}
\displaystyle
&\| v_h \|_{L^r(\Omega)} \\
&\leq c \sup_{w \in L^{r^{\prime}}(\Omega)}  \frac{ \| w \|_{ L^{r^{\prime}}(\Omega)} |v_h|_{H^1(\mathbb{T}_h)} +  \| w \|_{ L^{r^{\prime}}(\Omega)} |  v_h |_{H^1(\mathbb{T}_h)} + \| w \|_{ L^{r^{\prime}}(\Omega)} |v_h|_{H^1(\mathbb{T}_h)}}{\| w \|_{ L^{r^{\prime}}(\Omega)}} \\
&\leq c  |v_h|_{H^1(\mathbb{T}_h)}.
\end{align*}
\qed
\end{pf*}
\end{comment}

\section{CR Finite Element Approximation}

\subsection{Finite Element Approximation}
The CR finite element problem is to find $u_h^{CR}  \in CR_{h0}^1$ such that
\begin{align}
\displaystyle
a_{0h}(u_h^{CR} , \varphi_h) = (f , \varphi_h) \quad \forall \varphi_h \in CR_{h0}^1, \label{CR1}
\end{align}
where $a_{0h}: ( CR_{h0}^1 + H_0^1(\Omega) ) \times ( CR_{h0}^1 + H_0^1(\Omega) ) \to \mathbb{R}$ is defined by
\begin{align*}
\displaystyle
a_{0h}(\psi_h,\varphi_h) := \sum_{T \in \mathbb{T}_h} \int_T \nabla \psi_h \cdot \nabla \varphi_h dx = (\nabla_h \psi_h , \nabla_h \varphi_h).
\end{align*}
This problem is nonconforming because $CR_{h0}^1 \not\subset H_0^1(\Omega)$.

For the CR approximate solution $u_h^{CR}  \in CR_{h0}^1$ of \eqref{CR1}, we have the a priori estimate, using \eqref{pre22},
\begin{align*}
\displaystyle
| u_h^{CR} |_{H^1(\mathbb{T}_h)}^2
\leq \| f \| \| u_h^{CR} \|
\leq C(\Omega) \| f \| | u_h^{CR} |_{H^1(\mathbb{T}_h)}.
\end{align*}
By the Lax--Milgram lemma, there exists a unique solution $u_h^{CR} \in CR_{h0}^1$ for any $f \in L^2(\Omega)$.

\subsection{Classical Error Analysis}
The starting point for error analysis is the Second Strang Lemma, e.g. see \cite[Lemma 2.25]{ErnGue04},
\begin{align}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)} \leq 2 \inf_{v_h \in CR_{h0}^1} | u -v_h |_{H^1(\mathbb{T}_h)} + \sup_{\varphi_h \in CR_{h0}^1} \frac{a_{0h}(u,\varphi_h) - (f,\varphi_h)}{| \varphi_h|_{H^1(\mathbb{T}_h)}}. \label{CR2}
\end{align}

The first term of the inequality \eqref{CR2} is estimated as follows. Using the CR interpolation error estimate \eqref{modify12}, we have, for any $u \in H^2(\Omega)$,
\begin{align}
\displaystyle
\inf_{v_h \in CR_{h0}^1} | u -v_h |_{H^1(\mathbb{T}_h)}
&\leq | u -I_h^{CR} u |_{H^1(\mathbb{T}_h)} \leq c H |u|_{H^2(\Omega)}. \label{CR3}
\end{align}

\begin{comment}
In order to estimate the nonconforming term, we first introduce the trace inequality in \cite[Lemma 1.49]{PieErn12}.
\begin{lem} \label{CR=lem4}
Let $T \in \mathbb{T}_h$ and let $F$ be a face of $T$. We then have the trace inequality
\begin{align}
\displaystyle
\| w \|^2_{L^2(F)}
\leq  \frac{h_T}{\ell_F}  \| w \|_{L^2(T)} \left( 2 \| \nabla w \|_{L^2(T)^d} + d h_T^{-1} \| w \|_{L^2(T)} \right) \quad \forall w \in H^1(\mathbb{T}_h). \label{CR3}
\end{align}
Here, $\ell_{F} := \min_{x \in F} |x - x_F|$, where $x_F$ is the vertex of $T$ opposite to $F$.
\end{lem}

\begin{pf*}
Let $w \in H^1(\mathbb{T}_h)$. let $T \in \mathbb{T}_h$ and let $F$ be a face of $T$. We set
\begin{align*}
\displaystyle
\sigma_F(x) := \frac{|F|}{d |T|} (x - x_F),
\end{align*}
The normal component of $\sigma_F$ is constant and equal to one on $F$, and it is zero on all the other faces of $T$. From the divergence theorem, we have
\begin{align*}
\displaystyle
\| w \|^2_{L^2(F)}
&= \int_F |w|^2 ds
= \int_{\partial T} |w|^2  \sigma_F(x) \cdot n_T ds
= \int_T \mathop{\mathrm{div}} (|w|^2 \sigma_F(x)) dx \\
&= \int_T 2 w \sigma_F(x) \cdot \nabla w dx + \int_T |w|^2 \mathop{\mathrm{div}} \sigma_F(x) dx.
\end{align*}
Using
\begin{align*}
\displaystyle
\| \sigma_F \|_{L^{\infty}(T)^d} \leq \frac{|F|}{d |T|} h_T, \quad \mathop{\mathrm{div}} \sigma_F =  \frac{|F|}{|T|},
\end{align*}
and the H$\rm{\ddot{o}}$lder inequality, we obtain
\begin{align*}
\displaystyle
\| w \|^2_{L^2(F)}
\leq  \frac{|F|}{d |T|} h_T \| w \|_{L^2(T)} \left( 2 \| \nabla w \|_{L^2(T)^d} + d h_T^{-1} \| w \|_{L^2(T)} \right).
\end{align*}
Since $\frac{|F|}{|T|} = \frac{d}{\ell_F}$, we have
\begin{align*}
\displaystyle
\| w \|^2_{L^2(F)}
\leq  \frac{h_T}{\ell_F}  \| w \|_{L^2(T)} \left( 2 \| \nabla w \|_{L^2(T)^d} + d h_T^{-1} \| w \|_{L^2(T)} \right).
\end{align*}
\qed
\end{pf*}

Using the trace inequality, we have the estimates,\cite[Lemma 3.32]{ErnGue04}.
\begin{lem} \label{CR=lem5}
For $T \in \mathbb{T}_h$, $\varphi \in H^1(T)$, and a face $F \in \partial T$, set $\overline{\varphi} := \frac{1}{|F|} \int_F \varphi ds$. Then, there exists $c$ such that
\begin{align}
\displaystyle
\| \varphi - \overline{\varphi} \|_{L^2(F)} \leq c   \left( \frac{d}{\ell_F} \right)^{1/2} h_T  |\varphi|_{H^1(T)} \  \forall h, \ \forall T \in \mathbb{T}_h, \ \forall F \in \partial T, \ \forall \varphi \in H^1(T). \label{CR4}
\end{align}
\end{lem}
\end{comment}

From the standard scaling argument, we have a consistency error inequality, e.g., see \cite[Lemma 3.36]{ErnGue04}.
\begin{lem}[Asymptotic Consistency]  \label{CR=lem6}
Let $u \in H^1_0(\Omega) \cap H^2(\Omega)$ be the solution of the homogeneous Dirichlet Poisson problem \eqref{intro1}. It then holds that
\begin{align}
\displaystyle
\frac{a_{0h}(u,\varphi_h) - (f,\varphi_h)}{| \varphi_h|_{H^1(\mathbb{T}_h)}} \leq c \left( \sum_{T \in \mathbb{T}_h} \frac{h_T^4}{ (\min_{F \in \partial \mathbb{T}_h} \ell_F)^2} |u|^2_{H^2(T)} \right)^{1/2} \ \forall h, \ \forall \varphi_h \in  CR_{h0}^1, \label{CR4}
\end{align}
where $\partial \mathbb{T}_h$ denotes the set of all faces $F$ of $T \in \mathbb{T}_h$. Here, $\ell_F$ denotes the distance of the vertex of $T$ opposite to $F$ to the face.
\end{lem}

\begin{pf*}
We follow \cite[Lemma 3.36]{ErnGue04}.

Let $\varphi_h \in  CR_{h0}^1$. Because $- \varDelta u = f$, we have
\begin{align*}
\displaystyle
a_{0h}(u,\varphi_h) - (f,\varphi_h)
&= \sum_{T \in \mathbb{T}_h} \int_T ( \nabla u \cdot \nabla \varphi_h - f \varphi_h ) dx \\
&= \sum_{T \in \mathbb{T}_h} \sum_{F \in \partial \mathbb{T}_h} \int_F (n_T \cdot \nabla ) u \varphi_h ds.
\end{align*}
Because each face $F$ of an element $T$ located inside $\Omega$ appears twice in the above sum, we have
\begin{align*}
\displaystyle
a_{0h}(u,\varphi_h) - (f,\varphi_h)
&= \sum_{T \in \mathbb{T}_h} \sum_{F \in \partial \mathbb{T}_h} \int_F (n_T \cdot \nabla ) u \left( \varphi_h - \overline{\varphi_h} \right) ds
\end{align*}
with the mean value
\begin{align*}
\displaystyle
\overline{\varphi_h} := \frac{1}{|F|} \int_F \varphi_h ds.
\end{align*}
Furthermore, we get
\begin{align*}
\displaystyle
a_{0h}(u,\varphi_h) - (f,\varphi_h)
&= \sum_{T \in \mathbb{T}_h} \sum_{F \in \partial \mathbb{T}_h} \int_F n_T \cdot \left(  \nabla u - \overline{\nabla u} \right)\left( \varphi_h - \overline{\varphi_h} \right) ds
\end{align*}
with the mean value
\begin{align*}
\displaystyle
n_T \cdot \overline{\nabla u} := \frac{1}{|F|} \int_F (n_T \cdot \nabla ) u ds.
\end{align*}
The Cauchy--Schwarz inequality yields
\begin{align*}
\displaystyle
a_{0h}(u,\varphi_h) - (f,\varphi_h)
&\leq  \sum_{T \in \mathbb{T}_h} \sum_{F \in \partial \mathbb{T}_h} \| \nabla u -\overline{\nabla u} \|_{L^2(F)^d} \| \varphi_h - \overline{\varphi_h} \|_{L^2(F)}.
\end{align*}
For $F \in \partial \mathbb{T}_h$, let $\widehat{T} \subset \mathbb{R}^d$ be the reference simplex and let $\Phi_T: \widehat{T} \to T$ be the corresponding affine transformation with Jacobian matrix $A_T$. Let $\widehat{F} = \Phi_T^{-1} ( F )$. Using the standard scaling argument and the trace theorem on the reference element, we have
\begin{align*}
\displaystyle
\| \varphi_h - \overline{\varphi_h} \|_{L^2(F)}
&\leq \left( \frac{|F|}{|\widehat{F}|} \right)^{1/2} \| \hat{\varphi}_h - \overline{\hat{\varphi}_h} \|_{L^2(\widehat{F})}
\leq c \left( \frac{|F|}{|\widehat{F}|} \right)^{1/2} \| \hat{\varphi}_h - \overline{\hat{\varphi}_h} \|_{H^1(\widehat{T})}.
\end{align*}
The Deny--Lions Lemma (see \cite[Lemma B.67]{ErnGue04}) implies
\begin{align*}
\displaystyle
 \| \hat{\varphi}_h - \overline{\hat{\varphi}_h} \|_{H^1(\widehat{T})}
 \leq c | \hat{\varphi}_h |_{H^1(\widehat{T})}.
\end{align*}
Using the standard scaling argument again, we obtain
\begin{align*}
\displaystyle
\| \varphi_h - \overline{\varphi_h} \|_{L^2(F)}
&\leq c \left( \frac{|F|}{|\widehat{F}|} \right)^{1/2}  | \hat{\varphi}_h |_{H^1(\widehat{T})} \\
&\leq  c \left( \frac{|F|}{|\widehat{F}|} \right)^{1/2} \| A_T \|_2 \left( \frac{|\widehat{T}|}{|{T}|} \right)^{1/2}  | {\varphi}_h |_{H^1({T})} \\
&\leq c \left( \frac{|F|}{|{T}|} \right)^{1/2} h_T | {\varphi}_h |_{H^1({T})} = c \left( \frac{d}{\ell_F} \right)^{1/2} h_T | {\varphi}_h |_{H^1({T})}.
\end{align*}
Here,  $\| A_T \|_2$ denotes the matrix $2$-norm as
\begin{align*}
\displaystyle
\| A_T \|_2 := \sup_{0 \neq x \in \mathbb{R}^d} \frac{|A_T x|}{|x|},
\end{align*}
where $| x | := ( \sum_{i=1}^d |x_i|^2)^{1/2}$ for $x \in \mathbb{R}^d$.

By analogous argument, we have
\begin{align*}
\displaystyle
 \| \nabla u -\overline{\nabla u} \|_{L^2(F)^d}
 &\leq c \left( \frac{d}{\ell_F} \right)^{1/2} h_T |u|_{H^2(T)}.
\end{align*}
We consequently get
\begin{align*}
\displaystyle
a_{0h}(u,\varphi_h) - (f,\varphi_h)
&\leq c \sum_{T \in \mathbb{T}_h} \sum_{F \in \partial \mathbb{T}_h} \frac{h_T^2}{\ell_F}  |u|_{H^2(T)} | {\varphi}_h |_{H^1({T})} \\
&\leq c \sum_{T \in \mathbb{T}_h} \frac{h_T^2}{\min_{F \in \partial \mathbb{T}_h} \ell_F}  |u|_{H^2(T)} | {\varphi}_h |_{H^1({T})} \\
&\leq c \left( \sum_{T \in \mathbb{T}_h} \frac{h_T^4}{(\min_{F \in \partial \mathbb{T}_h} \ell_F)^2}  |u|_{H^2(T)}^2  \sum_{T \in \mathbb{T}_h} | {\varphi}_h |_{H^1({T})}^2  \right)^{1/2},
\end{align*}
which leads to \eqref{CR4}.
\qed
\end{pf*}

From \eqref{CR2}, \eqref{CR3} and \eqref{CR4}, we have
\begin{align*}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
\leq c H |u|_{H^2(\Omega)} + c \left( \sum_{T \in \mathbb{T}_h} \frac{h_T^4}{ (\min_{F \in \partial \mathbb{T}_h} \ell_F)^2} |u|^2_{H^2(T)} \right)^{1/2}.
\end{align*}

Since the order of the nonconforming term does not necessary becomes the order $H$, this inequality may be overestimated.
\begin{ex*}
Let  $0 \< h_T \leq 1$. As examples, we consider two cases.
\begin{description}
  \item[(\Roman{lone})] When we use meshes including the tetrahedra $T$ with vertices $(0,0,0)^T$, $(h_T,0,0)^T$, $(0,h_T,0)^T$, and $(0,0,h_T^{\varepsilon})^T$, we have
 \begin{align*}
\displaystyle
 \frac{h_T^4}{(\min_{F \in \partial \mathbb{T}_h} \ell_F)^2} |u|^2_{H^2(T)}\leq c h_T^{2(2-\varepsilon)} |u|^2_{H^2(T)},
\end{align*}
where $1 \< \varepsilon \leq 2$. Since $H = \mathcal{O}(h)$, we get
\begin{align*}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
\leq c (h + h^{2 - \varepsilon}) |u|_{H^2(\Omega)}.
\end{align*}
  \item[(\Roman{ltwo})] When we use meshes including the tetrahedra $T$ with vertices $(0,0,0)^T$, $(h_T,0,0)^T$, $(0,h_T,0)^T$, and $(h_T^{\gamma},0,h_T^{\varepsilon})^T$, we have
 \begin{align*}
\displaystyle
 \frac{h_T^4}{(\min_{F \in \partial \mathbb{T}_h} \ell_F)^2} |u|^2_{H^2(T)} \leq c h_T^{2(2-\varepsilon)} |u|^2_{H^2(T)},
\end{align*}
where $1 \< \gamma \< \varepsilon \leq 1 + \gamma$ and  $1 \< \varepsilon \leq 2$. Since $H = \mathcal{O}(h^{1 + \gamma - \varepsilon})$, we get
\begin{align*}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
\leq c (h^{1 + \gamma - \varepsilon} + h^{2 - \varepsilon}) |u|_{H^2(\Omega)}.
\end{align*}
\end{description}
\end{ex*}

\begin{comment}

In general, the classical estimate is not valid on anisotropic meshes including, for example, the triangle $T_1$ with vertices $(0,0)^T$, $(h,0)^T$, and $(0,h^{\varepsilon})^T$, where $0 \< h \leq 1$ and $1 \leq \varepsilon$, and/or the triangle $T_2$ with vertices $(0,0)^T$, $(h,0)^T$, and $(h^{\gamma},h^{\varepsilon})^T$, where $0 \< h \leq 1$ and $1 \leq \gamma \< \varepsilon \leq 1 + \gamma$.

\begin{pf*}
Let $\varphi_h \in  CR_{h0}^1$. Since $f = - \varDelta u$,
\begin{align*}
\displaystyle
a_{0h}(u,\varphi_h) - (f,\varphi_h) = \sum_{T \in \mathbb{T}_h} \int_{\partial T} \nabla u \cdot n_T \varphi_h ds = \frac{1}{2} \sum_{F \in \mathcal{F}_h} \int_{F} \nabla u \cdot n_F [[\varphi_h]] ds.
\end{align*}
Here, each side $F \in \mathcal{F}_h$ is either side of elements $T_1,T_2 \in \mathbb{T}_h$ or part of $\partial \Omega$, The function $\nabla u \cdot n_{T_1}$ on $F \subset T_1$ has in view of the continuity of $\nabla u$ opposite sign $\nabla u \cdot n_{T_2}$ on $F \subset T_2$. On each side $F \in \mathcal{F}_h$ by the special properties of $\varphi_h \in  CR_{h0}^1$, the jump $[[ \varphi_h ]]$ has mean value zero. We consequently have
\begin{align*}
\displaystyle
\sum_{F \in \mathcal{F}_h} \int_{F} \nabla u \cdot n_F [[\varphi_h]] ds
&= \sum_{F \in \mathcal{F}_h} \int_{F} \nabla u \cdot n_F \left( [[\varphi_h]] - \overline{[[\varphi_h]]} \right) ds \\
&= \sum_{F \in \mathcal{F}_h} \int_{F} \left( \nabla u - \overline{\nabla u} \right) \cdot n_F  \left( [[\varphi_h]] - \overline{[[\varphi_h]]} \right) ds,
\end{align*}
with the mean values
\begin{align*}
\displaystyle
 \overline{\nabla u} \cdot n_F := \frac{1}{|F|} \int_F \nabla u \cdot n_F ds, \quad  \overline{[[\varphi_h]]} := \frac{1}{|F|} \int_F {[[\varphi_h]]} ds.
\end{align*}
Splitting $[[ \varphi_h ]]$ into the contributions of the two neighbouring elements $T_1$, $T_2$ results in
\begin{align*}
\displaystyle
& \int_{F} \left( \nabla u - \overline{\nabla u} \right) \cdot n_F  \left( [[\varphi_h]] - \overline{[[\varphi_h]]} \right) ds \\
 &\quad =  \int_{F} \left( \nabla u - \overline{\nabla u} \right) \cdot n_F  \left( \varphi_h|_{T_1} - \overline{\varphi_h}_F \right) ds \\
 &\quad \quad - \int_{F} \left( \nabla u - \overline{\nabla u} \right) \cdot n_F  \left( \varphi_h|_{T_2} - \overline{\varphi_h}_F \right) ds.
\end{align*}
From Lemma \ref{CR=lem5}, we have
\begin{align*}
\displaystyle
\left \| \left( \nabla u - \overline{\nabla u} \right) \cdot n_F \right \|_{L^2(F)}
&\leq  c   \left( \frac{d}{\ell_F} \right)^{1/2} h_{T_1} |u|_{H^2(T_1)}, \\
\left \| \varphi_h|_{T_1} - \overline{\varphi_h}_F \right \|_{L^2(F)}
&\leq  c   \left( \frac{d}{\ell_F} \right)^{1/2} h_{T_1} |\varphi_h|_{H^1(T_1)},
\end{align*}
and analogously for the element $T_2$. From the foregoing results, it holds that for $F = T_1 \cap T_2$,
\begin{align*}
\displaystyle
&\left| \int_{F} \nabla u \cdot n_F [[\varphi_h]] ds \right| \\
&\quad \leq \left \| \left( \nabla u - \overline{\nabla u} \right) \cdot n_F \right \|_{L^2(F)}  \left( \left \| \varphi_h|_{T_1} - \overline{\varphi_h}_F \right \|_{L^2(F)}  + \left \| \varphi_h|_{T_2} - \overline{\varphi_h}_F \right \|_{L^2(F)}  \right) \\
&\quad \leq c \frac{h_{T_1}^2}{\ell_F} |u|_{H^2(T_1)} \left( |\varphi_h|_{H^1(T_1)} +  |\varphi_h|_{H^1(T_2)} \right) \\
&\quad \leq  c \frac{h_{T_1}^2}{\ell_F} |u|_{H^2(\Omega)}  |\varphi_h|_{H^1(\mathbb{T}_h)},
\end{align*}
and correspondingly $F \subset \partial \Omega$. From this, we conclude
\begin{align*}
\displaystyle
&\left| \sum_{F \in \mathcal{F}_h} \int_{F} \nabla u \cdot n_F [[\varphi_h]] ds \right|
\leq c  |u|_{H^2(\Omega)}  |\varphi_h|_{H^1(\mathbb{T}_h)}  \frac{h^2}{\min_{F \in \mathcal{F}_h} \ell_F}.
\end{align*}
\qed
\end{pf*}
\end{comment}

\begin{comment}
The order of the nonconforming term does not necessary becomes the order $H$. In general, the classical estimate is not true on anisotropic meshes, e.g., a right-angled triangle. We give examples in the two-dimensional case.
\begin{ex*}
Let us consider the triangle $T$ with vertices $(0,0)^T$, $(2h,0)^T$, and $(h,h^{\varepsilon})^T$, where $0 \< h \leq 1$ and $1 \leq \varepsilon \leq 2$. We then have
\begin{align*}
\displaystyle
\frac{h^2}{\ell_F} = h^{2 - \varepsilon}, \quad H = \frac{2h (h^2 + h^{2 \varepsilon})}{|T|} \leq c h^{2 - \varepsilon} \quad \text{with $|T| = h^{1+\varepsilon}$}.
\end{align*}
This case is asymptotically consistent.
\end{ex*}

\begin{ex*}
Let us consider the triangle $T$ with vertices $(0,0)^T$, $(h,0)^T$, and $(0,h^{\varepsilon})^T$, where $0 \< h \leq 1$ and $1 \leq \varepsilon$. We then have
\begin{align*}
\displaystyle
\frac{h^2}{\ell_F} = h^{2 - \varepsilon}, \quad H = \frac{h^{1 + \varepsilon} \sqrt{h^2 + h^{2 \varepsilon}}}{|T|} \leq c h  \quad \text{with $|T| = \frac{1}{2} h^{1+\varepsilon}$}.
\end{align*}
This case is not asymptotically consistent.
\end{ex*}

\begin{ex*}
Let us consider the triangle $T$ with vertices $(0,0)^T$, $(h,0)^T$, and $(h^{\gamma},h^{\varepsilon})^T$, where $0 \< h \leq 1$ and $1 \leq \gamma \< \varepsilon \leq 1 + \gamma$. We then have
\begin{align*}
\displaystyle
\frac{h^2}{\ell_F} = h^{2 - \varepsilon}, \ H = \frac{h \sqrt{h^{2 \gamma} + h^{2 \varepsilon}} \sqrt{( h - h^{\gamma})^2 + h^{2 \varepsilon}}}{|T|} \leq c h^{1 + \gamma - \varepsilon}  \ \text{with $|T| = \frac{1}{2} h^{1+\varepsilon}$}.
\end{align*}
This case is not asymptotically consistent.
\end{ex*}
\end{comment}

\subsection{Argument via the RT Interpolation Error}
To overcome the difficulty, we use the relation \eqref{pre21} in Lemma \ref{pre=lem2}, e.g., see also \cite{AcoDur99,LiuKik18}.

\begin{lem}[Asymptotic Consistency]  \label{CR=lem7}
We assume that $\Omega$ is convex. Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Let $u \in H^1_0(\Omega) \cap H^2(\Omega)$ be the solution of the homogeneous Dirichlet Poisson problem \eqref{intro1}. Then, there exists $c$, independent of $H$, such that
\begin{align}
\displaystyle
& \sup_{\varphi_h \in CR_{h0}^1} \frac{a_{0h}(u,\varphi_h) - (f,\varphi_h)}{| \varphi_h|_{H^1(\mathbb{T}_h)}}
\leq c H \| f \|.  \label{CR5}
\end{align}
\end{lem}

\begin{pf*}
Using \eqref{pre21}, we have, for any $w_h \in RT_h^0$,
\begin{align*}
\displaystyle
\sup_{\varphi_h \in CR_{h0}^1} \frac{a_{0h}(u,\varphi_h) - (f,\varphi_h)}{| \varphi_h|_{H^1(\mathbb{T}_h)}} =
 \sup_{\varphi_h \in CR_{h0}^1} \frac{ (\nabla u - w_h , \nabla_h \varphi_h) - ( \mathop{\mathrm{div}} w_h + f, \varphi_h  ) }{|\varphi_h|_{H^1(\mathbb{T}_h)}}.
\end{align*}
We set $w_h := I_h^{RT} \nabla u$. From Lemma \ref{pre=lem1}, we get
\begin{align*}
\displaystyle
 \mathop{\mathrm{div}} ( I_h^{RT} \nabla u ) = \Pi_h^0 \mathop{\mathrm{div}} (\nabla u) = - \Pi_h^0 f.
\end{align*}
Furthermore, we have, for any $\varphi_h \in CR_{h0}^1$,
\begin{align*}
\displaystyle
  ( - \Pi_h^0 f + f,  \Pi_h^0 \varphi_h) &= 0.
\end{align*}
We thus obtain
\begin{align*}
\displaystyle
&(\nabla u - I_h^{RT} \nabla u, \nabla_h \varphi_h) - ( - \Pi_h^0 f + f, \varphi_h  ) \\
&= (\nabla u - I_h^{RT} \nabla u, \nabla_h \varphi_h) - ( - \Pi_h^0 f + f, \varphi_h  - \Pi_h^0 \varphi_h) \\
&\leq \| \nabla u - I_h^{RT} \nabla u \|_{L^2(\Omega)^d} |\varphi_h|_{H^1(\mathbb{T}_h)} + \| f -  \Pi_h^0 f \|  \| \varphi_h  - \Pi_h^0 \varphi_h \| \\
&\leq c H |u|_{H^2(\Omega)}  |\varphi_h|_{H^1(\mathbb{T}_h)} + c h \| f \| |\varphi_h|_{H^1(\mathbb{T}_h)}.
\end{align*}
\qed
\end{pf*}

We consequently obtain the error estimate of the CR finite element method on anisotropic meshes.
\begin{thr} \label{CR=thr4}
We assume that $\Omega$ is convex. Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Let $u \in H^1_0(\Omega) \cap H^2(\Omega)$ be the solution of the homogeneous Dirichlet Poisson problem \eqref{intro1} with data $f \in L^2(\Omega)$. Let $u_h^{CR} \in CR_{h0}^1$ be the approximate solution of \eqref{CR1}. Then, there exists $c$, independent of $H$, such that
\begin{align}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)} \leq c H \| f \|. \label{CR6}
\end{align}
\end{thr}

\begin{pf*}
Using \eqref{CR2}, \eqref{modify12} and \eqref{CR5}, we have
\begin{align*}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
&\leq 2 \inf_{v_h \in CR_{h0}^1} | u -v_h |_{H^1(\mathbb{T}_h)} + \sup_{\varphi_h \in CR_{h0}^1} \frac{a_{0h}(u,\varphi_h) - (f,\varphi_h)}{| \varphi_h|_{H^1(\mathbb{T}_h)}} \\
&\leq 2  | u -I_h^{CR} u |_{H^1(\mathbb{T}_h)} + c H \| f \|  \leq c H \| f \|,
\end{align*}
which leads to the estimate \eqref{CR6}.
\qed
\end{pf*}

We next give the $L^2$ error estimate of the CR finite element method on anisotropic meshes, see also \cite{LiuKik07,LiuKik18,BreSco08}.

\begin{thr} \label{CR=thr5}
We assume that $\Omega$ is convex. Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Let $u \in H^1_0(\Omega) \cap H^2(\Omega)$ be the solution of the homogeneous Dirichlet Poisson problem \eqref{intro1} with data $f \in L^2(\Omega)$. Let $u_h^{CR} \in CR_{h0}^1$ be the approximate solution of \eqref{CR1}. Then, there exists $c$, independent of $H$, such that
\begin{align}
\displaystyle
\| u - u_h^{CR} \| \leq c H^2 \| f \|. \label{CR7}
\end{align}
\end{thr}

\begin{pf*}
We set $e_h :=  u - u_h^{CR} $. Let  $z \in H^2(\Omega) \cap H_0^1(\Omega)$ satisfy
\begin{align}
\displaystyle
a_0(\varphi , z ) = ( \varphi , e_h) \quad \forall \varphi \in H_0^1(\Omega) \label{CR8}
\end{align}
and $z_h^{CR} \in CR_{h0}^1$ satisfy
\begin{align}
\displaystyle
a_{0h}(\varphi_h , z_h^{CR}) = (\varphi_h , e_h) \quad \forall \varphi_h \in CR_{h0}^1. \label{CR9}
\end{align}
We then have
\begin{align}
\displaystyle
\| e_h \|^2 &= (e_h , e_h ) = a_{0h}(u , z )  - a_{0h}(u_h^{CR} , z_h^{CR}) \notag \\
&= a_{0h}(u -  u_h^{CR} , z - z_h^{CR} ) + a_{0h}( u -  u_h^{CR} ,  z_h^{CR}) + a_{0h}( u_h^{CR} ,  z - z_h^{CR} ) \notag \\
&= a_{0h}(u -  u_h^{CR} , z - z_h^{CR} ) \notag \\
&\quad + a_{0h}( u -  u_h^{CR} ,  z_h^{CR} - I_h^{CR} z) +   a_{0h}( u -  u_h^{CR} ,   I_h^{CR} z) \notag\\
&\quad + a_{0h}( u_h^{CR} - I_h^{CR} u ,  z - z_h^{CR} ) +a_{0h}( I_h^{CR} u ,  z - z_h^{CR} ).  \label{CR10}
\end{align}

Using Theorem \ref{CR=thr4}, the first term on the right hand side of \eqref{CR10} can be estimated as
\begin{align}
\displaystyle
 a_{0h}( u - u_h^{CR} ,  z- z_h^{CR})
 &\leq | u - u_h^{CR} |_{H^1(\mathbb{T}_h)}   | z- z_h^{CR} |_{H^1(\mathbb{T}_h)} \notag \\
 &\leq c H^2  \| f \| \|  e_h\|.  \label{CR11}
\end{align}
For the second and fourth terms on the right hand side of \eqref{CR10}, we have
\begin{align}
\displaystyle
& a_{0h}( u -  u_h^{CR} ,  z_h^{CR} - I_h^{CR} z) \notag \\
 &\quad =  a_{0h}( u -  u_h^{CR} ,  z_h^{CR} -  z) +  a_{0h}( u -  u_h^{CR} ,  z - I_h^{CR} z) \notag \\
 &\quad \leq  | u - u_h^{CR} |_{H^1(\mathbb{T}_h)} \left(  | z_h^{CR} - z |_{H^1(\mathbb{T}_h)} +  | z- I_h^{CR} z |_{H^1(\mathbb{T}_h)} \right) \notag \\
 &\quad \leq  c H^2  \| f \| \|  e_h\|, \label{CR12}
\end{align}
and, analogously,
\begin{align}
\displaystyle
a_{0h}( u_h^{CR} - I_h^{CR} u ,  z - z_h^{CR} )
&\leq c H^2 \| f \| \|  e_h\|. \label{CR13}
\end{align}
From \eqref{CR8}, \eqref{CR9} and \eqref{pre21}, we have
\begin{align*}
\displaystyle
&a_{0h}( u -  u_h^{CR} ,   I_h^{CR} z) \\
&\quad = a_{0h}( u ,   I_h^{CR} z) - a_{0h}( u_h^{CR} ,   I_h^{CR} z)  = (\nabla u , \nabla_h I_h^{CR} z) - (f ,  I_h^{CR} z) \\
&\quad = (\nabla u , \nabla_h I_h^{CR} z - \nabla z) - (f ,  I_h^{CR} z - z) + (\nabla u , \nabla z) - (f , z) \\
&\quad = (\nabla u - I_h^{RT} \nabla u , \nabla_h I_h^{CR} z - \nabla z) - (f + \mathop{\mathrm{div}} ( I_h^{RT} \nabla u ),  I_h^{CR} z - z).
\end{align*}
From Lemma \ref{pre=lem1} and $ \mathop{\mathrm{div}} ( I_h^{RT} \nabla u ) = - \Pi_h^0 f$, we have
\begin{align}
\displaystyle
&a_{0h}( u -  u_h^{CR} ,   I_h^{CR} z) \notag \\
&\quad =  (\nabla u - I_h^{RT} \nabla u , \nabla_h I_h^{CR} z - \nabla z) - (f - \Pi_h^0 f,  I_h^{CR} z - z) \notag\\
&\quad \leq \| \nabla u - I_h^{RT} \nabla u \|_{L^2(\Omega)^d} | I_h^{CR} z - z |_{H^1(\mathbb{T}_h)} + \| f - \Pi_h^0 f \| \|  I_h^{CR} z - z \|  \notag\\
&\quad \leq c H^2 \| f \| \|  e_h\|.  \label{CR14}
\end{align}

Analogously, from $ \mathop{\mathrm{div}} ( I_h^{RT} \nabla z ) = - \Pi_h^0 e_h$, we have
\begin{align}
\displaystyle
&a_{0h}( I_h^{CR} u ,  z - z_h^{CR} ) \notag \\
&\quad=  (\nabla_h I_h^{CR} u - \nabla u , \nabla z - I_h^{RT} \nabla z ) - (I_h^{CR} u - u , e_h + \mathop{\mathrm{div}} ( I_h^{RT} \nabla z)) \notag\\
&\quad \leq | I_h^{CR} u - u |_{H^1(\mathbb{T}_h)} \|  \nabla z - I_h^{RT} \nabla z \|_{L^2(\Omega)^d} + \|  I_h^{CR} u - u \| \| e_h  - \Pi_h^0 e_h \| \notag \\
&\quad \leq c H^2 \| f \| \|  e_h\|. \label{CR15}
\end{align}
Combining \eqref{CR10}, \eqref{CR11}, \eqref{CR12}, \eqref{CR13}, \eqref{CR14}, and \eqref{CR15}, we finally get
\begin{align*}
\displaystyle
\| e_h \|^2 \leq c H^2 \| f \| \|  e_h\|,
\end{align*}
which leads to the target estimate.
\qed
\end{pf*}

\section{RT Finite Element Error Estimates}

\subsection{Dual mixed formulation of the Poisson problem}
The Poisson equation \eqref{intro1} $- \varDelta u =  - \mathop{\mathrm{div}} \nabla u = f$ can be written as the following system. Find $(\sigma,u):\Omega \to \mathbb{R}^{d} \times \mathbb{R}$ such that
\begin{subequations} \label{mix1}
\begin{align}
\displaystyle
\sigma - \nabla u &= 0  \quad \text{in $\Omega$}, \label{mix1a}\\
\mathop{\mathrm{div}} \sigma &=  - f  \quad \text{in $\Omega$}, \label{mix1b}\\
u &= 0 \quad \text{on $\partial \Omega$}. \label{mix1c}
\end{align}
\end{subequations}
We consider the following dual mixed formulation: Find $(\sigma , u) \in H(\mathop{\mathrm{div}};\Omega) \times L^2(\Omega)$ such that
\begin{subequations} \label{mix2}
\begin{align}
\displaystyle
a(\sigma,v) + b(v,u) &= 0 \quad \forall v \in  H(\mathop{\mathrm{div}};\Omega), \label{mix2a}\\
b(\sigma,q) &= - (f ,q) \quad \forall q \in L^2(\Omega), \label{mix2b}
\end{align}
\end{subequations}
where bilinear forms $a: H(\mathop{\mathrm{div}};\Omega) \times H(\mathop{\mathrm{div}};\Omega) \to \mathbb{R}$ and $b:H(\mathop{\mathrm{div}};\Omega) \times L^2(\Omega) \to \mathbb{R}$ are defined by
\begin{align*}
\displaystyle
&a(\sigma,v) := (\sigma , v), \quad b(v,q) := (\mathop{\mathrm{div}} v , q).
\end{align*}
We set $X_0 := \{ v \in H(\mathop{\mathrm{div}};\Omega)  ; \ b(v,q) = 0 \ \forall q \in L^2(\Omega) \}$. Because there exists a constant $c \> 0$ such that
\begin{align*}
\displaystyle
a(v,v) \geq c \| v \|^2_{H(\mathop{\mathrm{div}};\Omega) } \quad \forall v \in X_0
\end{align*}
and the bilinear form $b(_\cdot , _ \cdot)$ satisfies the inf--sup condition
\begin{align}
\displaystyle
\inf_{0 \neq q \in L^2(\Omega)} \sup_{0 \neq v \in H(\mathop{\mathrm{div}};\Omega)  } \frac{b(v,q)}{\| v \|_{H(\mathop{\mathrm{div}};\Omega) } \| q \|}
\geq \beta_* \> 0, \label{mix3}
\end{align}
\eqref{mix2} is uniquely solvable; e.g., see \cite{GirRav86,BofBreFor13}.

\subsection{RT Approximate Problem}
We consider the following RT approximate problem. Find $(\sigma_h^{RT} , u_h^{RT}) \in RT^0_{h} \times M_{h}^0$ such that
\begin{subequations} \label{mix4}
\begin{align}
\displaystyle
a(\sigma_h^{RT},v_h) + b(v_h,u_h^{RT}) &= 0, \quad \forall v_h \in  RT^0_{h}, \label{mix4a}\\
b(\sigma_h^{RT},q_h) &= - (f ,q_h), \quad \forall q_h \in M_{h}^0.\label{mix4b}
\end{align}
\end{subequations}
This setting is conforming because $RT^0_{h} \times M_{h}^0 \subset H(\mathop{\mathrm{div}};\Omega) \times L^2(\Omega)$. It is given later that the discrete inf--sup condition
\begin{align*}
\displaystyle
\inf_{q_h \in M_{h}^0} \sup_{v_h \in RT^0_{h}} \frac{b(v_h , q_h)}{\| v_h \|_{H(\mathop{\mathrm{div}};\Omega)} \| q_h \| } \geq c_* \> 0
\end{align*}
holds, where $c_*$ is a constant independent of $h$.

\subsection{Error Estimates of the RT Finite Element Approximation}
This section gives error estimates of the mixed finite element approximation \eqref{mix4}. We  emphasise that we do not impose the shape regularity condition and the maximum-angle condition for the mesh partition. That is, we assume that $\{ \mathbb{T}_h \}$ is a family of conformal meshes satisfying Assumption \ref{pre=ass1}.

\begin{lem} \label{mix=lem8}
Let $D \subset \mathbb{R}^d$ be a bounded domain. For any $g \in L^2(D)$, there exists $v \in H^1(D)^d$ such that
\begin{align}
\displaystyle
\mathop{\mathrm{div}} v = g \quad \text{in $D$} \label{mix5}
\end{align}
and
\begin{align}
\displaystyle
| v |_{H^1(D)^d} \leq  \| g \|_{L^2(D)}, \quad \| v \|_{L^2(\Omega)^d} \leq  C_P(D) \| g \|_{L^2(D)},  \label{mix6}
\end{align}
where $C_P(D)$ is the Poincar\'e constant.
\end{lem}

\begin{pf*}
The proof can be found in  \cite[Lemma 2.2]{BofBre08}.
\qed
\end{pf*}

\begin{comment}
\begin{pf*}
Let $B \subset \mathbb{R}^d$ be a ball containing $D$. We set
\begin{align*}
\displaystyle
\tilde{g} :=
\begin{cases}
g \quad \text{in $D$},\\
0 \quad \text{in $B \backslash D$}.
\end{cases}
\end{align*}
There then exists a unique solution ${u} \in H_0^1(B) \cap H^2(B)$ such that
\begin{align*}
\displaystyle
\varDelta {u} = \tilde{g} \quad \text{in $B$}, \quad {u} = 0 \quad \text{on $\partial B$}.
\end{align*}
 It is known that $u$ satisfies the a priori estimate
\begin{align*}
\displaystyle
| {u} |_{H^2(D)} \leq | {u} |_{H^2(B)} \leq \| \varDelta {u} \|_{L^2(B)} = \| \tilde{g} \|_{L^2(B)} = \| {g} \|_{L^2(D)}.
\end{align*}
We also get the a priori estimate
\begin{align*}
\displaystyle
| {u} |_{H^1(D)} \leq | {u} |_{H^1(B)} \leq C_P(D) \| \tilde{g} \|_{L^2(B)} =  C_P(D) \| {g} \|_{L^2(D)}.
\end{align*}
Therefore, setting $v := \nabla {u} \in H^1(D)^d$, \eqref{mix5} and \eqref{mix6} hold.
\qed
\end{pf*}
\end{comment}

We next give the discrete inf--sup condition.
\begin{lem}[Discrete inf--sup condition] \label{mix=lem9}
If $ C_{G}^{RT} H \leq 1$, there exists a constant $c_*$, depending only on the Poincar\'e constant, such that
\begin{align}
\displaystyle
\inf_{q_h \in M_{h}^0} \sup_{v_h \in RT^0_{h}} \frac{b(v_h , q_h)}{\| v_h \|_{H(\mathop{\mathrm{div}};\Omega)} \| q_h \|} \geq c_* \> 0,\label{mix7}
\end{align}
where $ C_{G}^{RT}$ is the constant appearing in Corollary \ref{modify=coro3}.
\end{lem}

\begin{pf*}
Let $q_h \in  M_{h}^0$. From Lemma \ref{mix=lem8}, there exists $v \in H^1(\Omega)^d$ such that $\mathop{\mathrm{div}} v = q_h$ in $\Omega$, $| v |_{H^1(\Omega)^d} \leq   \| q_h \|$, and $\| v \|_{L^2(\Omega)^d} \leq C_P(\Omega) \| q_h \|$.

By the Gauss theorem, we have
\begin{align*}
\displaystyle
\sum_{T \in \mathbb{T}_h} \int_{\partial T} v \cdot n_T ds = \sum_{T \in \mathbb{T}_h} \int_{T} \mathop{\mathrm{div}} v dx = \int_{\Omega} q_h dx.
\end{align*}
From the definition of the Raviart--Thomas interpolation, we conclude that
\begin{align*}
\displaystyle
 \int_{\Omega} \mathop{\mathrm{div}} ( I_T^{RT} v) p_h dx &= \sum_{T \in \mathbb{T}_h} p_h \int_T \mathop{\mathrm{div}} ( I_T^{RT} v) dx = \sum_{T \in \mathbb{T}_h} p_h \int_{\partial T} n_T \cdot (I_T^{RT} v) ds \\
&= \sum_{T \in \mathbb{T}_h} p_h \int_{\partial T} v \cdot n_T ds = \int_{\Omega} q_h p_h dx \quad \forall p_h \in M_h^0.
\end{align*}
Therefore,  it follows that $\mathop{\mathrm{div}} (I_h^{RT} v) = q_h$.

From the definitions, we have
\begin{align*}
\displaystyle
\| I_h^{RT} v \|^2_{H(\mathop{\mathrm{div}};\Omega)}
&= \| I_h^{RT} v \|_{L^2(\Omega)^d}^2 + \| \mathop{\mathrm{div}} (I_h^{RT} v) \|^2 \\
&\leq 2 \| I_h^{RT} v - v \|_{L^2(\Omega)^d}^2 + 2 \| v \|_{L^2(\Omega)^d}^2 + \| q_h \|^2 \\
&\leq 2 ( C_{G}^{RT})^2 H^2 |v|^2_{H^1(\Omega)^d} + 2 C_P(\Omega)^2 \| q_h \|^2 + \| q_h \|^2 \\
&\leq \left( 3 + 2 C_P(\Omega)^2 \right) \| q_h \|^2.
\end{align*}
We thus have
\begin{align*}
\displaystyle
 \sup_{v_h \in RT^0_{h}} \frac{b(v_h , q_h)}{\| v_h \|_{H(\mathop{\mathrm{div}};\Omega)}}
& \geq \frac{b(I_h^{RT} v , q_h)}{\| I_h^{RT} v \|_{H(\mathop{\mathrm{div}};\Omega)}}
 \geq \frac{1}{\left( 3 + 2 C_P(\Omega)^2 \right)^{1/2}} \frac{(q_h,q_h)}{\| q_h \|},
\end{align*}
and the proof of \eqref{mix7} is completed with $c_* := \left( 3 + 2 C_P(\Omega)^2 \right)^{- 1/2}$.
\qed
\end{pf*}

From the discrete equations \eqref{mix4} and their continuous counterpart \eqref{mix2}, we obtain the Galerkin orthogonality
\begin{subequations} \label{mix8}
\begin{align}
\displaystyle
a(\sigma - \sigma_h^{RT} , v_h) + b(v_h , u - u_h^{RT}) &= 0 \quad \forall v_h \in RT^0_{h}, \label{mix8a} \\
b(\sigma - \sigma_h^{RT} , q_h) &= 0 \quad \forall q_h \in M_{h}^0. \label{mix8b}
\end{align}
\end{subequations}
We then get the following C\'ea-lemma-type estimates with the help of \eqref{mix8} and the inf--sup condition \eqref{mix7}.

\begin{thr} \label{mix=thr6}
Let $\sigma \in H^1(\Omega)^d$ and $\sigma_h^{RT} \in RT_h^0$ be the solutions of \eqref{mix1} and  \eqref{mix4}, respectively. We then have
\begin{align}
\displaystyle
\| \sigma - \sigma_h^{RT}\|_{L^2(\Omega)^d} \leq  \| \sigma - I_h^{RT} \sigma \|_{L^2(\Omega)^d}. \label{mix9}
\end{align}
Furthermore, let $(\sigma,u) \in H^1(\Omega)^d \times L^2(\Omega)$ and $(\sigma_h^{RT},u_h^{RT}) \in RT_h^0 \times M_h^0$ be the solutions of \eqref{mix1} and  \eqref{mix4}, respectively. Then, if $C_{G}^{RT} H \leq 1$,  it holds that
\begin{align}
\displaystyle
\| u - u_h^{RT} \|
\leq \| u - \Pi_h^0 u \| + c_*^{-1} \|  \sigma - \sigma_h^{RT}  \|_{L^2(\Omega)^d}. \label{mix11}
\end{align}
Here, $C_{G}^{RT}$ and $c_*$ are respectively the constants appearing in Corollary \ref{modify=coro3} and Lemma \ref{mix=lem9}.
\end{thr}

\begin{pf*}
The proof can be found in \cite[Lemma 3.7, Lemma 3.9]{BofBre08}.
\qed
\end{pf*}

Using Theorem \ref{mix=thr6} and the interpolation error estimates of Corollary \ref{modify=coro1} and \ref{modify=coro3}, we thus have the error estimates of the mixed finite element approximation \eqref{mix4} on anisotropic meshes violating the maximum-angle condition.

\begin{thr} \label{mix=thr7}
let $(\sigma,u) \in H^1(\Omega)^d \times H^1(\Omega)$ and $(\sigma_h^{RT},u_h^{RT}) \in RT_h^0 \times M_h^0$ be the solutions of \eqref{mix1} and  \eqref{mix4}, respectively. Then, there exists a constant $c_1 \> 0$, independnt of  $\sigma$, $H$, and the geometric properties of $\mathbb{T}_h$,  such that
\begin{align}
\displaystyle
\| \sigma - \sigma_h^{RT}\|_{L^2(\Omega)^d}
&\leq c_1 H | \sigma |_{H^1(\Omega)^d}.  \label{mix10}
\end{align}
Furthermore, if $C_{G}^{RT} H \leq 1$, there exists a constant $c_2 \> 0$, depending on the discrete inf--sup condition but independent of $\sigma$, $u$, $h$, $H$, and the geometric properties of $\mathbb{T}_h$
\begin{align}
\displaystyle
\| u - u_h^{RT} \|
&\leq c_2 \left(h |u|_{H^1(\Omega)} + H |\sigma|_{H^1(\Omega)^d} \right). \label{mix12}
\end{align}
Here, $C_{G}^{RT}$ is the constant appearing in Corollary \ref{modify=coro3}.
\end{thr}

\begin{comment}
\begin{pf*}
From \eqref{mix8a} and the definition of the $L^2$-projection, we have
\begin{align*}
\displaystyle
a(\sigma - \sigma_h^{RT} , v_h) + b(v_h , \Pi_h^0 u - u_h^{RT}) &= 0 \quad \forall v_h \in RT^0_{h}.
\end{align*}
With help of the discrete inf--sup stability \eqref{mix7} and H$\rm{\ddot{o}}$lder's inequality,
\begin{align*}
\displaystyle
\| u - u_h^{RT} \|
&\leq \| u - \Pi_h^0 u \| + \| \Pi_h^0 u - u_h^{RT} \| \\
&\leq \| u - \Pi_h^0 u \| + c_*^{-1}  \sup_{v_h \in RT^0_{h}} \frac{b(v_h , \Pi_h^0 u - u_h^{RT})}{\| v_h \|_{H(\mathop{\mathrm{div}};\Omega)} } \\
&=   \| u - \Pi_h^0 u \| + c_*^{-1}  \sup_{v_h \in RT^0_{h}} \frac{a(\sigma_h^{RT} - \sigma, v_h)}{\| v_h \|_{H(\mathop{\mathrm{div}};\Omega(} } \\
&\leq  \| u - \Pi_h^0 u \| + c_*^{-1} \| \sigma_h^{RT} - \sigma \|_{L^2(\Omega)^d},
\end{align*}
which concludes \eqref{mix11}.

Furthermore, using the inequality \eqref{pre14} in Corollary \ref{pre=coro1} and \eqref{mix10}, we have \eqref{mix12}.
\qed
\end{pf*}
\end{comment}

\section{Relationship between the RT and CR Finite Element Approximation}
This section shows the relationship between the RT and CR problems. Find $(\bar{\sigma}_h^{RT} , \bar{u}_h^{RT}) \in RT^0_{h} \times M_{h}^0$ such that
\begin{subequations} \label{RTCR1}
\begin{align}
\displaystyle
a(\bar{\sigma}_h^{RT},v_h) + b(v_h,\bar{u}_h^{RT}) &= 0 \quad \forall v_h \in  RT^0_{h}, \label{RTCR1a}\\
b(\bar{\sigma}_h^{RT},q_h) &= - (\Pi_h^0 f ,q_h) \quad \forall q_h \in M_{h}^0 \label{RTCR1b}
\end{align}
\end{subequations}
and find $\bar{u}_h^{CR} \in CR_{h0}^1$ such that
\begin{align}
\displaystyle
a_{0h}(\bar{u}_h^{CR} , \varphi_h) = (\Pi_h^0 f , \varphi_h) \quad \forall \varphi_h \in CR_{h0}^1. \label{RTCR2}
\end{align}
Here, \eqref{RTCR2} is the CR approximation of the Poisson equation
\begin{align}
\displaystyle
- \varDelta \bar{u} = \Pi_h^0 f \quad \text{in $\Omega$}, \quad \bar{u} = 0 \quad \text{on $\partial \Omega$}. \label{RTCR3}
\end{align}

In the case of $d=2$, it is well known that there exists a relationship between $(\bar{\sigma}_h^{RT} , \bar{u}_h^{RT}) $ and $\bar{u}_h^{CR}$ introduced by Marini; for example, \cite{Mar85}. See also \cite{LiuKik07,KikSai16,LiuKik18}. We here show the relation in the three dimensional case.

Let us consider a tetrahedron $T \subset \mathbb{R}^3$ such as that in Figure \ref{fig1}. Let $x_i$ ($i=1,2,3,4$) be the vertices and $m_{i,j}$ the midpoints of edges of the tetrahedron; that is, $m_{i,j} := \frac{1}{2} (x_i + x_j)$. Furthermore, for $1 \leq i \leq 4$, let $F_i$ be the face of the tetrahedron opposite $x_i$. Then, by simple calculation, we find the equality
\begin{align*}
\displaystyle
L := \sum_{i=1}^4 |x_i - x_T|^2 =  |m_{1,4} - m_{2,3}|^2 +  |m_{1,3} - m_{2,4}|^2 +  |m_{1,2} - m_{3,4}|^2,
\end{align*}
holds, where $x_T$ is the barycentre of $T$ such that $x_T := \frac{1}{4} \sum_{i=1}^4 x_i$.

\begin{figure}
\includegraphics[width=7cm]{tetrahedron2.png}
\caption{Tetrahedron}
\label{fig1}
\end{figure}

We present a quadrature scheme over a simplex $T \subset \mathbb{R}^3$ (e.g., \cite[p.307]{Str71}) that is easily conformed.
\begin{lem} \label{RTCR=lem10}
For any $f \in \mathcal{C}^0(T)$, the quadrature scheme
\begin{align*}
\displaystyle
\int_T f(x) dx &\sim - \frac{|T|}{20} \sum_{i=1}^4 f(x_i) + \frac{|T|}{5} \sum_{1 \leq i \textless j \leq 4} f(m_{i,j})
\end{align*}
is exact for polynomials of degree less than or equal to $2$;
\begin{align}
\displaystyle
\int_T f(x) dx +  \frac{|T|}{20} \sum_{i=1}^4 f(x_i) - \frac{|T|}{5} \sum_{1 \leq i \< j \leq 4}f(m_{i,j}) = 0 \quad \forall f \in \mathcal{P}^2(T). \label{RTCR4}
\end{align}
\end{lem}

Define the function $\varphi_T$ by
\begin{align}
\displaystyle
\varphi_T(x) &:=
\begin{cases} \label{RTCR5}
L - 12 |x - x_T|^2, \quad \text{on $T$},\\
0, \quad \text{otherwise}.
\end{cases}
\end{align}
We then have the following lemma.

\begin{lem}  \label{RTCR=lem11}
It holds that
\begin{align}
\displaystyle
\frac{1}{|F_i|} \int_{F_i} \varphi_T (x) ds &= 0, \quad i=1,2,3,4, \label{RTCR6} \\
\frac{1}{|T|} \int_T \varphi_T (x) dx &= \frac{2}{5} L, \label{RTCR7} \\
\frac{1}{|T|} \int_T | \nabla \varphi_T (x)|^2 dx &=  \frac{144}{5} L.\label{RTCR8}
\end{align}
\end{lem}

\begin{pf*}
From second-order three-point numerical integration over $F_1$,
\begin{align*}
\displaystyle
\int_{F_1} f(x) ds &= \frac{|F_1|}{3} \left( f (m_{2,3}) + f (m_{3,4}) + f (m_{2,4}) \right)  \quad \forall f \in \mathcal{P}^2(T),
\end{align*}
we have
\begin{align*}
\displaystyle
&\frac{1}{|F_1|} \int_{F_1} \varphi_T (x) ds \\
&= \frac{1}{3} \left( \varphi_T (m_{2,3}) + \varphi_T (m_{3,4}) + \varphi_T (m_{2,4}) \right) \\
&= \frac{1}{3} \left( 3 L - 12 \left(  |m_{2,3} - x_T|^2 + |m_{3,4} - x_T|^2 + |m_{2,4} - x_T|^2 \right) \right) \\
&= \frac{1}{3} \left( 3 L - \frac{12}{4} \left(  |m_{2,3} - m_{1,4}|^2 + |m_{3,4} - m_{1,2}|^2 + |m_{2,4} - m_{1,3}|^2 \right) \right) = 0,
\end{align*}
which leads to \eqref{RTCR6}.

Next, using \eqref{RTCR4}, we have
\begin{align*}
\displaystyle
&\frac{1}{|T|} \int_T \varphi_T (x) dx \\
&= - \frac{1}{20} \sum_{i=1}^4 \varphi_T(x_i) + \frac{1}{5} \sum_{1 \leq i \< j \leq 4}\varphi_T(m_{i,j}) \\
&= - \frac{1}{20} \left( 4 L - 12 \sum_{i=1}^4 |x_i - x_T|^2 \right) + \frac{1}{5} \left( 6L - 12 \sum_{1 \leq i \< j \leq 4} |m_{i,j} - x_T|^2 \right) \\
&= \frac{2}{5} L,
\end{align*}
which leads to \eqref{RTCR7}. We here used
\begin{align*}
\displaystyle
&\sum_{1 \leq i \< j \leq 4} |m_{i,j} - x_T|^2  \\
&=  |m_{1,2} - x_T|^2 +  |m_{1,3} - x_T|^2 +  |m_{1,4} - x_T|^2 \\
&\quad +  |m_{2,3} - x_T|^2 +  |m_{2,4} - x_T|^2 +  |m_{3,4} - x_T|^2 \\
&= \frac{1}{4} \left( 2 |m_{1,2} - m_{3,4}|^2 + 2 |m_{1,3} - m_{2,4}|^2 +2  |m_{1,4} - m_{2,3}|^2  \right) = \frac{L}{2}.
\end{align*}
We similarly obtain
\begin{align*}
\displaystyle
&\frac{1}{|T|} \int_T | \nabla \varphi_T (x)|^2 dx \\
&= \frac{24^2}{|T|} \int_T |x - x_T|^2 dx \\
&= - \frac{24^2}{20} \sum_{i=1}^4 |x_i - x_T|^2 + \frac{24^2}{5} \sum_{1 \leq i \< j \leq 4} |m_{i,j} - x_T|^2 = \frac{144}{5} L,
\end{align*}
which leads to \eqref{RTCR8}.
\qed
\end{pf*}

We set the bubble space $B_h$ by
\begin{align}
\displaystyle
B_h := \{ b_h \in L^2(\Omega); \ b_h |_T \in \mathop{\mathrm{span}} \{ \varphi_T\}, \ \forall T \in \mathbb{T}_h \}. \label{RTCR9}
\end{align}
Then, for any $\psi_h \in CR_{h0}^1$ and $b_h \in B_h$, because one writes $b_h|_T = c_b\varphi_T$ for $c_b \in \mathbb{R}$, it holds that
\begin{align*}
\displaystyle
(\nabla_h \psi_h, \nabla_h b_h)  &= \sum_{T \in \mathbb{T}_h} c_b \int_T \nabla \psi_h \cdot \nabla \varphi_T dx \\
&= \sum_{T \in \mathbb{T}_h} c_b \left\{ \sum_{F \subset \partial T} (n_F \cdot \nabla \psi_h ) \int_F \varphi_T ds - \int_T \varDelta \psi_h \varphi_T dx \right\} = 0.
\end{align*}
We here used the facts that \eqref{RTCR6}, $n_F \cdot \nabla \psi_h$ is constant on $F$, and $\varDelta \psi_h = 0$ on $T$. That is to say, two finite element spaces $CR_{h0}^1$ and $B_h$ are orthogonal to  each other.

Furthermore, we define the finite element space $X_h^{bCR}$ by
\begin{align}
\displaystyle
X_h^{bCR} := CR_{h0}^1 + B_h = \{ \psi_h + b_h ; \ \psi_h \in CR_{h0}^1, \ b_h \in B_h  \}. \label{RTCR10}
\end{align}
We consider the following finite element problem. Find $u_h^{bCR} \in X_h^{bCR}$ such that
\begin{align}
\displaystyle
a_{0h}(u_h^{bCR} , \varphi_h) = (\nabla_h u_h^{bCR} , \nabla_h \varphi_h) = (\Pi_h^0 f , \varphi_h) \quad \forall \varphi_h \in X_{h}^{bCR}.\label{RTCR11}
\end{align}
The solution $u_h^{bCR} \in X_h^{bCR}$ is then decomposed as $u_h^{bCR} = \bar{u}_h^{CR} + b_h$ with $\bar{u}_h^{CR} \in  CR_{h0}^1$ and $b_h\in  B_h$. Note that $\bar{u}_h^{CR}$ and $b_h$  respectively satisfy \eqref{RTCR2} and the equation
\begin{align}
\displaystyle
a_{0h}(b_h , c_h) = (\nabla_h b_h , \nabla_h c_h) = (\Pi_h^0 f , c_h) \quad \forall c_h \in  B_h. \label{RTCR12}
\end{align}
On each element $T \in \mathbb{T}_h$, \eqref{RTCR12} has the form
\begin{align*}
\displaystyle
\gamma_T \int_T \nabla \varphi_T \cdot \nabla \varphi_T dx = \int_T \Pi_T^0 f \varphi_T dx, \quad \gamma_T \in \mathbb{R}.
\end{align*}
From \eqref{RTCR7} and \eqref{RTCR8}, we have
\begin{align}
\displaystyle
\gamma_T = \frac{1}{72} \Pi_T^0 f \quad \forall T \in \mathbb{T}_h. \label{RTCR13}
\end{align}
\begin{thr}  \label{RTCR=thr8}
Let ${u}_h^{bCR} \in X_{h}^{bCR}$ be the solution of \eqref{RTCR11} and $(\bar{\sigma}_h^{RT},\bar{u}_h^{RT}) \in RT^0_{h} \times M_{h}^0$ the solution of \eqref{RTCR1}. We then have $\nabla_h u_h^{bCR} \in RT^0_{h}$ and
\begin{align}
\displaystyle
\bar{\sigma}_h^{RT} &= \nabla u_h^{bCR} \quad\forall T \in \mathbb{T}_h,  \label{RTCR14}\\
\bar{u}_h^{RT} &= \Pi_T^0  u_h^{bCR}  \quad\forall T \in \mathbb{T}_h.\label{RTCR15}
\end{align}
\end{thr}

\begin{pf*}
The proof can be found in \cite{HuMa15}.
\qed
\end{pf*}

From Theorem \ref{RTCR=thr8}, for $d=3$, the following lemma holds.
\begin{lem}  \label{RTCR=lem12}
Let $\overline{u}_h^{CR} \in CR_{h0}^1$ be the solution of \eqref{RTCR2} and $(\bar{\sigma}_h^{RT},\bar{u}_h^{RT}) \in RT^0_{h} \times M_{h}^0$ be the solution of \eqref{RTCR1}. We then have the relationships
\begin{align}
\displaystyle
\bar{\sigma}_h^{RT}|_T &= \nabla \bar{u}_h^{CR} - \frac{1}{3} \Pi_T^0 f (x - x_T) \quad \forall T \in \mathbb{T}_h, \label{RTCR16} \\
\bar{u}_h^{RT}|_T &= \Pi_T^0  \bar{u}_h^{CR} + \frac{1}{180} \Pi_T^0 f  \sum_{i=1}^4 |x_i - x_T|^2  \quad \forall T \in \mathbb{T}_h. \label{RTCR17}
\end{align}
\end{lem}

Using relationship between the RT and CR finite element methods, we have the error estimate of the CR finite element approximation with the bubble function.

\begin{lem} \label{RTCR=lem13}
We assume that $\Omega$ is convex. Let $\{ \mathbb{T}_h \}$ be a family of conformal meshes satisfying Assumption \ref{pre=ass1}. Let $\bar{u} \in  H_0^1(\Omega) \cap H^2(\Omega)$ be the solution of \eqref{RTCR3} and $u_h^{bCR} \in X_{h}^{bCR}$ be the solution of the CR problem \eqref{RTCR11}. There then exists  a constant $c \> 0$ independent of $\bar{u}$, $h$, $H$ and the geometric properties of $\mathbb{T}_h$ such that
\begin{align}
\displaystyle
| \bar{u} - u_h^{bCR} |_{H^1(\mathbb{T}_h)} \leq  c H \| \Pi_h^0 f \|. \label{RTCR18}
\end{align}
\end{lem}

\begin{pf*}
Let $(\bar{\sigma}_h^{RT},\bar{u}_h^{RT}) \in RT^0_{h} \times M_{h}^0$ be the solution of \eqref{RTCR1}. From Theorem \ref{RTCR=thr8}, it holds that $\nabla_h u_h^{bCR} \in RT^0_{h}$ and $\bar{\sigma}_h^{RT} = \nabla_h u_h^{bCR}$. Setting $\bar{\sigma} := \nabla \bar{u} \in H^1(\Omega)^d$, we then have, using inequality \eqref{mix10}, that
\begin{align*}
\displaystyle
| \bar{u} - u_h^{bCR} |_{H^1(\mathbb{T}_h)}
&=  \left( \sum_{T \in \mathbb{T}_h} \| \bar{\sigma} - \bar{\sigma}_h^{RT} \|^2_{L^2(T)^d} \right)^{1/2} \\
&\leq c H |\bar{\sigma}|_{H^1(\Omega)^d} = c H |\bar{u}|_{H^2(\Omega)} \leq c H \| \Pi_h^0 f \|.
\end{align*}
\qed
\end{pf*}

\begin{comment}

From Lemma \ref{RTCR=lem13}, the equivalence of the RT and CR problems, we have the following theorem.
\begin{thr} \label{RTCR=thr10}
Let $u \in  H_0^1(\Omega) \cap H^2(\Omega)$ be the solution of \eqref{intro1} and $u_h^{CR} \in CR_{h0}^1$ be the solution of the Crouzeix--Raviart problem \eqref{CR1}. It then holds that
\begin{align}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
&\leq  c_1 \left(  \| \nabla u - \bar{\sigma}_h^{RT}  \|_{L^2(\Omega)^d} + h \| f \|  \right) \leq c_2 H \| f \|, \label{RTCR19} \\
\| \nabla u - \bar{\sigma}_h^{RT}  \|_{L^2(\Omega)^d} &\leq | u - u_h^{CR} |_{H^1(\mathbb{T}_h)} + c_3 h \| f \|, \label{RTCR20}
 \end{align}
 where $c_1$, $c_2$ and $c_3$ are independent of $h$, $H$ and the geometric properties of $\mathbb{T}_h$
\end{thr}

\begin{pf*}
Let $(\bar{\sigma}_h^{RT},\bar{u}_h^{RT}) \in RT^0_{h} \times M_{h}^0$ be the solution of \eqref{RTCR1}, $\bar{u}_h^{CR} \in CR_{h0}^1$ be the solution of the Crouzeix--Raviart problem \eqref{RTCR2}, and $u_h^{bCR} \in X_{h}^{bCR}$ be the solution of the Crouzeix--Raviart problem \eqref{RTCR11}. Then, from Theorem \ref{RTCR=thr8}, it holds that $\nabla_h u_h^{bCR} \in RT^0_{h}$ and $\bar{\sigma}_h^{RT} = \nabla_h u_h^{bCR}$.

From the triangle inequality, we have
\begin{align*}
\displaystyle
&| u - u_h^{CR} |_{H^1(\mathbb{T}_h)} \\
&= \| \nabla u - \nabla_h u_h^{CR} \|_{L^2(\Omega)^d} \\
&\leq  \| \nabla u - \bar{\sigma}_h^{RT} \|_{L^2(\Omega)^d} +  \|  \bar{\sigma}_h^{RT} - \nabla_h \bar{u}_h^{CR} \|_{L^2(\Omega)^d} + \| \nabla_h  \bar{u}_h^{CR} - \nabla_h  u_h^{CR} \|_{L^2(\Omega)^d}.
\end{align*}

From \eqref{intro3} in the two-dimensional case and \eqref{RTCR16} in the three-dimensional case, we get
\begin{align*}
\displaystyle
\|  \bar{\sigma}_h^{RT} - \nabla_h \bar{u}_h^{CR} \|_{L^2(\Omega)^d}^2
&= \sum_{T \in \mathbb{T}_h} \int_T |\bar{\sigma}_h^{RT} - \nabla \bar{u}_h^{CR}|^2 dx \\
&= \frac{1}{d^2}\sum_{T \in \mathbb{T}_h} (\Pi_T^0 f)^2 \int_T |x - x_T|^2 dx \\
&\leq \frac{1}{d^2} \sum_{T \in \mathbb{T}_h} |T|^{-1} \| f \|^2_{0,2,T} h_T^2 |T|
= \frac{1}{d^2} \sum_{T \in \mathbb{T}_h} h_T^2 \| f \|^2_{0,2,T}.
\end{align*}
We here used H$\rm{\ddot{o}}$lder's inequality to obtain
\begin{align*}
\displaystyle
(\Pi_T^0 f)^2 = \frac{1}{|T|^2} \left( \int_T f dx \right)^2 \leq  \frac{1}{|T|^2} |T| \| f \|^2_{0,2,T} = |T|^{-1} \| f \|^2_{0,2,T}.
\end{align*}

For any $w_h \in CR_{h0}^1$, from \eqref{CR1} and \eqref{RTCR2}, we have
\begin{align*}
\displaystyle
a_{0h}(\bar{u}_h^{CR} - u_h^{CR} , w_h) &= (f - \Pi_h^0 f , w_h) = (f - \Pi_h^0 f , w_h - \Pi_h^0 w_h).
\end{align*}
Substituting $\bar{u}_h^{CR} - u_h^{CR} \in CR_{h0}^1$ into $w_h$ in the above equation, we have, using the error estimate of the $L^2$-projection \eqref{pre12},
\begin{align*}
\displaystyle
\| \nabla_h  \bar{u}_h^{CR} - \nabla_h  u_h^{CR} \|_{L^2(\Omega)^d}^2
&= (f - \Pi_h^0 f , \bar{u}_h^{CR} - u_h^{CR} - \Pi_h^0 (\bar{u}_h^{CR} - u_h^{CR} )) \\
&\leq  C_P^{L^2} h \| f - \Pi_h^0 f \| \| \nabla_h  \bar{u}_h^{CR} - \nabla_h  u_h^{CR} \|_{L^2(\Omega)^d} \\
&\leq  C_P^{L^2} h \| f \| \| \nabla_h  \bar{u}_h^{CR} - \nabla_h  u_h^{CR} \|_{L^2(\Omega)^d}.
\end{align*}
We here used $\| f - \Pi_h^0 f \| \leq \| f \|$ which is derived from the definition of the $L^2$-projection in the first inequality.

From Lemma \ref{RTCR=lem13}, we have
\begin{align*}
\displaystyle
 \| \nabla u - \bar{\sigma}_h^{RT}  \|_{L^2(\Omega)^d}
 &\leq  \| \nabla u -  \nabla \bar{u
 }  \|_{L^2(\Omega)^d} +  \| \nabla \bar{u} -  \nabla_h u_h^{bCR}  \|_{L^2(\Omega)^d} \\
 &\leq c_{-1}\| f - \Pi_h^0 f \|_{H^{-1}(\Omega)} +  \| \nabla \bar{u} -  \nabla_h u_h^{bCR}\|_{L^2(\Omega)^d} \\
 &\leq c_{-1} \| f - \Pi_h^0 f \|_{H^{-1}(\Omega)} + c H \| \Pi_h^0 f \|,
\end{align*}
where $c_{-1} := \| - \varDelta^{-1} \|_{\mathcal{L}(H^{-1}(\Omega) ; H^1_0(\Omega))}$. For any $\varphi \in H_0^1(\Omega)$, considering that $f - \Pi_h^0 f$ is $L^2$-orthogonal to $M_{h}^0$ and the error estimate of the $L^2$-projection \eqref{pre12}, we have
\begin{align*}
\displaystyle
(f - \Pi_h^0 f , \varphi)
&= (f - \Pi_h^0 f , \varphi - \Pi_h^0 \varphi) \\
&\leq \| f - \Pi_h^0 f \| \| \varphi - \Pi_h^0 \varphi \| \\
&\leq C_P^{L^2} h \| f \| | \varphi |_{H^1(\Omega)}
\end{align*}
and therefore
\begin{align*}
\displaystyle
\|f - \Pi_h^0 f \|_{H^{-1}(\Omega)}
&\leq C_P^{L^2} h \| f \|.
\end{align*}

Gathering the above estimates, we have
\begin{align*}
\displaystyle
| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}
 &\leq c_{-1} \| f - \Pi_h^0 f \|_{H^{-1}(\Omega)} + c H \| \Pi_h^0 f \| + \frac{1}{d} h \| f \| +  C_P^{L^2} h \| f \| \\
 &\leq C_P^{L^2} h \| f \| + c H \| f \| + \frac{1}{d} h \| f \| +  C_P^{L^2} h \| f \| \\
&\leq c H \| f \|.
\end{align*}
We again use \eqref{intro3} in the two-dimensional case and \eqref{RTCR16} in the three-dimensional case for the proof of \eqref{RTCR20}. That is, we have
\begin{align*}
\displaystyle
& \| \nabla u - \bar{\sigma}_h^{RT}  \|_{L^2(\Omega)^d} \\

 &\leq  \| \nabla u -  \nabla_h {u}_h^{CR} \|_{L^2(\Omega)^d} +  \| \nabla_h {u}_h^{CR} -  \nabla_h \bar{u}_h^{CR} \|_{L^2(\Omega)^d} +  \|  \nabla_h \bar{u}_h^{CR}  - \bar{\sigma}_h^{RT} \|_{L^2(\Omega)^d} \\
 &\leq \| \nabla u -  \nabla_h {u}_h^{CR} \|_{L^2(\Omega)^d} + c h \| f \|.
\end{align*}
\qed
\end{pf*}

\end{comment}

\section{Numerical Results}
This section presents results of numerical examples. Let $\Omega := (0,1)^3$. Let $u_h^L$ and $u_h^{CR}$ be the $\mathcal{P}^1$-Lagrange and $\mathcal{P}^1$-CR finite element solutions, respectively, for the model problem
\begin{align*}
\displaystyle
- \varDelta u &= 2y(1-y)z(1-z) + 2x(1-x)z(1-z) + 2x(1-x)y(1-y) \quad \text{in $\Omega$}, \\
u &= 0 \quad \text{on $\partial \Omega $},
\end{align*}
which is the exact solution $u = x(1-x)y(1-y)z(1-z)$.

\begin{figure}[htbp]

  \includegraphics[width=7cm]{M8L22.png}
  \caption{Mesh: $M=8$, $N=22$}
   \label{fig2}
   \vspace{0.5cm}
   \includegraphics[width=8cm]{3dmesh2.png}
  \caption{Elements}
   \label{fig3}
\end{figure}

Let $M$ be the division number of each side of the bottom face and $N$ the division number of the height of $\Omega$ with $N \sim M^{\gamma}$ (see Fig. \ref{fig2}). There are two elements as shown in Fig. \ref{fig3}.

If an exact solution $u$ is known, the error $e_h := u - u_h$ and $e_{h/2} := u - u_{h/2}$ are computed numerically for two mesh sizes $h$ and $h/2$. The convergence indicator $r$ is defined by
\begin{align*}
\displaystyle
r = \frac{1}{\log(2)} \log \left( \frac{\| e_h \|_X}{\| e_{h/2} \|_X} \right).
\end{align*}

We set $h := \frac{1}{M}$. The parameter $H$ is then $H = \mathcal{O}( h^{2 - \gamma} )$. We compute the convergence order with respect to $H_0^1$ and $L^2$ norms defined by
\begin{align*}
\displaystyle
&Err_h^L(H^1) := \frac{| u - u_h^L |_{H^1(\Omega)}}{\| \varDelta u \|}, \quad  Err_h^{L}(L^2) := \frac{\| u - u_h^L \|}{\| \varDelta u \|}, \\
&Err_h^{CR}(H^1) := \frac{| u - u_h^{CR} |_{H^1(\mathbb{T}_h)}}{\| \varDelta u \|}, \quad  Err_h^{CR}(L^2) := \frac{\| u - u_h^{CR} \|}{\| \varDelta u \|},
\end{align*}
for three cases: $\gamma = 1.5$, $\gamma = 1.9$ and $\gamma = 2.0$. In order to compute the above norms, we use the five-order fifteen-point numerical integration introduced in \cite{Kea86}. The results are give in Table \ref{table1}, Table \ref{table2} when $\gamma = 1.5$, Table \ref{table3}, Table \ref{table4} when $\gamma = 1.9$, and Table \ref{table5}, Table \ref{table6} when $\gamma = 2.0$. Further, $N_p^L$ and  $N_p^{CR}$ denote respectively the degrees of freedom for the $\mathcal{P}^1$-Lagrange finite element and the $\mathcal{P}^1$-CR finite element.

\begin{table}[htb]
\caption{Error of the $\mathcal{P}^1$-Lagrange finite element solution ($\gamma = 1.5$)}
\centering
\begin{tabular}{l|l|l|l|l|l|l|l|l} \hline
$M$& $N$ &  $h$ & $H$ & $N_p^L$ & $Err_h^L(H^1)$ & $r$& $Err_h^L(L^2)$  & $r$  \\ \hline \hline
4& 8 & 2.50e-01 & 5.00e-01  &225 &  1.2043e-01&   & 9.5321e-03&  \\
8& 22 & 1.25e-01 & 3.54e-01 & 1,863  & 7.0318e-02   & 0.78   & 3.1646e-03    & 1.59   \\
16& 64 & 6.25e-02 & 2.50e-01 & 18,785 & 4.4662e-02  &  0.65 & 1.2570e-03   & 1.33 \\
32& 182 &  3.13e-02 & 1.77e-01 & 199,287 & 2.9479e-02   & 0.60  & 5.4477e-04   & 1.21 \\
\hline
\end{tabular}
\label{table1}
\end{table}

\begin{table}[htb]
\caption{Error of the $\mathcal{P}^1$-CR finite element solution ($\gamma = 1.5$)}
\centering
\begin{tabular}{l|l|l|l|l|l|l|l|l} \hline
$M$& $N$ &  $h $ & $H$ & $N_p^{CR}$ & $Err_h^{CR}(H^1)$ & $r$ & $Err_h^{CR}(L^2)$  & $r$  \\ \hline \hline
4& 8 &2.50e-01  &5.00e-01 & 1,440 &8.2569e-02   &   & 3.8242e-03    &  \\
8&  22 & 1.25e-01 & 3.54e-01 & 14,912 & 4.0629e-02  & 1.02   & 8.8356e-04    & 2.11 \\
16& 64 & 6.25e-02 & 2.50e-01 & 168,448 &  2.0042e-02 & 1.02  & 2.0485e-04    & 2.11 \\
32& 182 &   3.13e-02  &1.77e-01 & 1,889,024  &  9.9579e-03 & 1.01  & 4.8960e-05   & 2.07 \\
\hline
\end{tabular}
\label{table2}
\end{table}

\begin{table}[htb]
\caption{Error of the $\mathcal{P}^1$-Lagrange finite element solution ($\gamma = 1.9$)}
\centering
\begin{tabular}{l|l|l|l|l|l|l|l|l} \hline
$M$& $N$ &  $h $ & $H$ & $N_p^L$ & $Err_h^L(H^1)$ & $r$ & $Err_h^L(L^2)$  & $r$  \\ \hline \hline
4& 14  & 2.50e-01 & 8.71e-01  & 345 &  1.4873e-01 &   & 1.4032e-02   &  \\
8& 52 & 1.25e-01 & 8.12e-01 & 4,293 & 1.2167e-01  &  0.29 &  9.3061e-03  & 0.59 \\
16& 194 & 6.25e-02 &  7.58e-01 & 56,355 & 1.0919e-01  & 0.16  &7.4989e-03    & 0.31\\
32& 724 & 3.13e-02 & 7.07e-01 & 789,525  & 1.0128e-01  &  0.11 & 6.4558e-03    & 0.22 \\
\hline
\end{tabular}
\label{table3}
\end{table}

\begin{table}[htb]
\caption{Error of the $\mathcal{P}^1$-CR finite element solution ($\gamma = 1.9$)}
\centering
\begin{tabular}{l|l|l|l|l|l|l|l|l} \hline
$M$& $N$ &  $h $ & $H$ & $N_p^{CR}$ & $Err_h^{CR}(H^1)$ & $r$ & $Err_h^{CR}(L^2)$  & $r$ \\ \hline \hline
4& 14 & 2.50e-01 & 8.71e-01 & 2,496 & 7.9756e-02 &   & 3.2993e-03   &  \\
8& 52 & 1.25e-01 & 8.12e-01 & 35,072  & 3.9708e-02   & 1.01   & 7.7177e-04  & 2.10 \\
16& 194 & 6.25e-02 & 7.58e-01 & 509,568  & 1.9814e-02  & 1.00  &   1.8781e-04 &2.04 \\
32& 724 &  3.13e-02 & 7.07e-01 & 7,508,480 &  9.9003e-03  &  1.00 & 4.6546e-05    & 2.01 \\
\hline
\end{tabular}
\label{table4}
\end{table}

\begin{table}[htb]
\caption{Error of the $\mathcal{P}^1$-Lagrange finite element solution ($\gamma = 2.0$)}
\centering
\begin{tabular}{l|l|l|l|l|l|l|l|l} \hline
$M$& $N$ &  $h$ & $H$ & $N_p^L$ & $Err_h^L(H^1)$ & $r$ & $Err_h^L(L^2)$  & $r$ \\ \hline \hline
4& 16  & 2.50e-01 & 1.00   & 425 & 1.5862e-01  &   & 1.5909e-02   &  \\
8& 64 & 1.25e-01 & 1.00  & 5,265 & 1.4079e-01  &  0.17 & 1.2472e-02   & 0.35 \\
16& 256 & 6.25e-02  & 1.00 & 74,273 & 1.3597e-01  & 0.05  & 1.1646e-02    & 0.10 \\
32& 1,024 & 3.13e-02 & 1.00 & 1,116,225  & 1.3474e-01   & 0.01  & 1.1442e-02   & 0.03 \\
\hline
\end{tabular}
\label{table5}
\end{table}

\begin{table}[htb]
\caption{Error of the $\mathcal{P}^1$-CR finite element solution ($\gamma = 2.0$)}
\centering
\begin{tabular}{l|l|l|l|l|l|l|l|l} \hline
$M$& $N$ &  $h $ & $H$& $N_p^{CR}$ & $Err_h^{CR}(H^1)$ & $r$ & $Err_h^{CR}(L^2)$  & $r$  \\ \hline \hline
4& 16  & 2.50e-01 & 1.00 & 2,848 & 7.9473e-02  &   &   3.2264e-03 &  \\
8& 64 & 1.25e-01 & 1.00  & 43,136  & 3.9647e-02  & 1.00  &  7.6153e-04  & 2.08  \\
16& 256  & 6.25e-02  & 1.00  & 672,256 & 1.9803e-02   & 1.00  & 1.8680e-04    & 2.03  \\
32& 1,024  & 3.13e-02 & 1.00 & 10,618,880  &  9.8984e-03 & 1.00  &  4.6458e-05  & 2.01 \\
\hline
\end{tabular}
\label{table6}
\end{table}

Observing the numerical results, the convergence indicators $r$ in each norms are respectively
\begin{align*}
\displaystyle
&| u - u_h^L |_{H^1(\Omega)} = \mathcal{O}(H), \quad \| u - u_h^L \| = \mathcal{O}(H^2), \\
&| u - u_h^{CR} |_{H^1(\mathbb{T}_h)} = \mathcal{O}(h), \quad \| u - u_h^{CR} \| = \mathcal{O}(h^2),
\end{align*}
where $H = \mathcal{O}( h^{2 - \gamma} )$. Meanwhile, the theoretical results are as follows:
\begin{align*}
\displaystyle
&| u - u_h^L |_{H^1(\Omega)} = \mathcal{O}(H), \quad \| u - u_h^L \| = \mathcal{O}(H^2), \\
&| u - u_h^{CR} |_{H^1(\mathbb{T}_h)} = \mathcal{O}(H), \quad \| u - u_h^{CR} \| = \mathcal{O}(H^2),
\end{align*}
if $\Omega$ is convex and $u \in H^2(\Omega) \cap H_0^1(\Omega)$. In this numerical examples, the CR finite element approximation is superior to the Lagrange finite element approximation on this anisotropic meshes. The theoretical explanation of this point is still open.

\begin{acknowledgements}
This work was supported by JSPS KAKENHI Grant Number \\ JP16H03950. We would like to thank the anonymous referee for the valuable comments.
\end{acknowledgements}