Gradient recovery for elliptic interface problem: III. Nitsche's method

TL;DR

This paper introduces a Nitsche's method-based gradient recovery technique for elliptic interface problems that maintains high accuracy near interfaces, overcoming issues caused by skinny triangles in unfitted meshes.

Contribution

It proposes a novel gradient recovery method using Nitsche's method, ensuring superconvergence and improved gradient accuracy near interfaces in unfitted mesh scenarios.

Findings

Achieves superconvergence of gradient recovery

Maintains accuracy near interfaces with skinny triangles

Validated by numerical examples

Abstract

This is the third paper on the study of gradient recovery for elliptic interface problem. In our previous works [H. Guo and X. Yang, 2016, arXiv:1607.05898 and {\it J. Comput. Phys.}, 338 (2017), 606--619], we developed {gradient recovery methods} for elliptic interface problem based on body-fitted meshes and immersed finite element methods. Despite the efficiency and accuracy that these methods bring to recover the gradient, there are still some cases in unfitted meshes where skinny triangles appear in the generated local body-fitted triangulation that destroy the accuracy of recovered gradient near the interface. In this paper, we propose a gradient recovery technique based on Nitsche's method for elliptic interface problem, which avoids the loss of accuracy of gradient near the interface caused by skinny triangles. We analyze the supercloseness between the gradient of the numerical solution by the Nitsche's method and the gradient of interpolation of the exact solution, which leads to the superconvergence of the proposed gradient recovery method. We also present several numerical examples to validate the theoretical results.

Figures

Tables

Table 1. Table 1 : Numerical results for Example 5.1 with β 1 = 10 , β 2 = 1 formulae-sequence subscript 𝛽 1 10 subscript 𝛽 2 1 \beta_{1}=10,\beta_{2}=1 .

$h$	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
1/16	4.61e-02	–	2.37e-02	–	1.82e-02	–
1/32	2.34e-02	0.98	9.34e-03	1.34	7.70e-03	1.25
1/64	1.17e-02	1.00	3.28e-03	1.51	2.75e-03	1.48
1/128	5.88e-03	1.00	1.17e-03	1.48	9.95e-04	1.47
1/256	2.94e-03	1.00	4.08e-04	1.52	3.36e-04	1.56
1/512	1.47e-03	1.00	1.43e-04	1.51	1.17e-04	1.53
1/1024	7.35e-04	1.00	5.08e-05	1.49	4.17e-05	1.48

Table 2. Table 2 : Numerical results for Example 5.1 with β 1 = 1000 , β 2 = 1 formulae-sequence subscript 𝛽 1 1000 subscript 𝛽 2 1 \beta_{1}=1000,\beta_{2}=1 .

$h$	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
1/16	4.19e-02	–	2.62e-02	–	2.15e-02	–
1/32	2.13e-02	0.98	9.98e-03	1.39	8.52e-03	1.33
1/64	1.06e-02	1.00	3.53e-03	1.50	3.09e-03	1.46
1/128	5.33e-03	1.00	1.25e-03	1.50	1.12e-03	1.47
1/256	2.66e-03	1.00	4.33e-04	1.52	3.75e-04	1.57
1/512	1.33e-03	1.00	1.52e-04	1.51	1.29e-04	1.54
1/1024	6.66e-04	1.00	5.41e-05	1.49	4.59e-05	1.49

Table 3. Table 3 : Numerical results for Example 5.1 with β 1 = 1 , β 2 = 100000 formulae-sequence subscript 𝛽 1 1 subscript 𝛽 2 100000 \beta_{1}=1,\beta_{2}=100000 .

$h$	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
1/16	1.99e-01	–	2.95e-02	–	3.23e-02	–
1/32	9.97e-02	1.00	9.94e-03	1.57	1.06e-02	1.61
1/64	4.98e-02	1.00	3.53e-03	1.50	3.08e-03	1.78
1/128	2.49e-02	1.00	1.19e-03	1.56	1.05e-03	1.55
1/256	1.25e-02	1.00	4.33e-04	1.46	3.85e-04	1.45
1/512	6.23e-03	1.00	1.56e-04	1.47	1.38e-04	1.48
1/1024	3.12e-03	1.00	5.51e-05	1.50	4.85e-05	1.51

Table 4. Table 4 : Numerical results for Example 5.1 with β 1 = 100000 , β 2 = 1 formulae-sequence subscript 𝛽 1 100000 subscript 𝛽 2 1 \beta_{1}=100000,\beta_{2}=1 .

$h$	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
1/16	4.19e-02	–	2.62e-02	–	2.15e-02	–
1/32	2.13e-02	0.98	9.99e-03	1.39	8.54e-03	1.33
1/64	1.06e-02	1.00	3.54e-03	1.50	3.10e-03	1.46
1/128	5.33e-03	1.00	1.25e-03	1.50	1.12e-03	1.46
1/256	2.66e-03	1.00	4.33e-04	1.52	3.84e-04	1.55
1/512	1.33e-03	1.00	1.52e-04	1.51	1.37e-04	1.48
1/1024	6.66e-04	1.00	5.41e-05	1.49	4.65e-05	1.56

Table 5. Table 5 : Numerical results for Example 5.2.

$h$	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
1/16	8.86e-02	–	5.81e-02	–	3.74e-02	–
1/32	3.90e-02	1.19	1.50e-02	1.95	1.19e-02	1.65
1/64	1.90e-02	1.04	4.37e-03	1.78	3.57e-03	1.74
1/128	9.48e-03	1.00	1.57e-03	1.48	1.29e-03	1.47
1/256	4.74e-03	1.00	5.63e-04	1.48	4.72e-04	1.45
1/512	2.37e-03	1.00	2.00e-04	1.50	1.71e-04	1.47
1/1024	1.18e-03	1.00	7.06e-05	1.50	6.62e-05	1.37

Table 6. Table 6 : Numerical results for Example 5.3.

DOF	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
2573	2.49e-02	–	1.32e-02	–	8.68e-03	–
10265	1.29e-02	0.48	4.42e-03	0.79	3.20e-03	0.72
41009	6.57e-03	0.49	1.56e-03	0.75	1.25e-03	0.68
163937	3.30e-03	0.50	5.09e-04	0.81	4.46e-04	0.74
655553	1.66e-03	0.50	1.74e-04	0.77	1.59e-04	0.74
2621825	8.28e-04	0.50	5.92e-05	0.78	5.48e-05	0.77

Table 7. Table 7 : Numerical results for Example 5.4.

DOF	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
1381	2.14e-01	–	9.04e-02	–	8.01e-02	–
5489	1.12e-01	0.47	2.51e-02	0.93	2.02e-02	1.00
21889	5.67e-02	0.49	7.60e-03	0.86	6.65e-03	0.80
87425	2.85e-02	0.50	2.32e-03	0.86	2.11e-03	0.83
349441	1.42e-02	0.50	7.29e-04	0.83	6.99e-04	0.80
1397249	7.12e-03	0.50	2.43e-04	0.79	2.39e-04	0.78

Table 8. Table 8 : Numerical results for Example 5.5.

DOF	$De$	order	$D^{i}e$	order	$D_r^{r}e$	order
6514	1.30e-01	–	5.36e-02	–	4.39e-02	–
26027	6.70e-02	0.48	1.53e-02	0.91	1.58e-02	0.74
104053	3.39e-02	0.49	4.23e-03	0.93	4.74e-03	0.87
416105	1.70e-02	0.50	1.17e-03	0.93	1.32e-03	0.92
1664209	8.51e-03	0.50	3.26e-04	0.92	3.72e-04	0.91
6656417	4.25e-03	0.50	9.26e-05	0.91	9.97e-05	0.95

Equations

$$-\nabla \cdot (\beta (z)\nabla u(z))$$

$$u$$

$$[[u]]$$

$$[[\beta {\partial }_{n}u]]$$

$$[[w]]=w_1-w_2,$$

$$\beta(z)=\left\{\begin{array}[]{ccc}\beta_{1}(z)&\text{if }z=(x,y)\in\Omega_{1},\\ \beta_{2}(z)&\text{if }z=(x,y)\in\Omega_{2},\\ \end{array}\right.$$

$$\Vert w{\Vert }_{k,p,D_1\cup D_2}={\left(\Vert w{\Vert }_{k,p,D_1}^p+\Vert w{\Vert }_{k,p,D_2}^p\right)}^{1/p},$$

$$|w{|}_{k,p,D_1\cup D_2}={\left(|w{|}_{k,p,D_1}^p+|w{|}_{k,p,D_2}^p\right)}^{1/p}.$$

$$\frac{h_T}{{\rho }_T}\le \sigma ,$$

$${\mathcal{T}}_{\Gamma ,h}=\left\{T\in {\mathcal{T}}_h:\Gamma \cap \overline{T}\ne \varnothing \right\},$$

$${\Omega }_{\Gamma ,h}=\bigcup _{T\in {\mathcal{T}}_{\Gamma ,h}}T.$$

$${\mathcal{T}}_{i,h}=\left\{T\in {\mathcal{T}}_h:\overline{{\Omega }_i}\cap \overline{T}\ne \varnothing \right\},i=1,2;$$

$${\Omega }_{i,h}=\bigcup _{T\in {\mathcal{T}}_{i,h}}T,{\omega }_{i,h}=\bigcup _{T\in {\mathcal{T}}_{i,h}\setminus {\mathcal{T}}_{\Gamma ,h}}T,i=1,2.$$

$$V_{i,h}=\left\{v\in C^0({\Omega }_{i,h}):v{|}_T\in {\mathbb{P}}_1(T)\text{ for any }T\in {\mathcal{T}}_{i,h}\right\},i=1,2,$$

$$V_h=\left\{v_h=(v_{1,h},v_{2,h}):v_{i,h}\in V_{i,h},i=1,2\right\},$$

$$V_{h,0}=\left\{v_h\in V_h:v_h{|}_{\partial \Omega }=0\right\}.$$

$${\kappa }_1{|}_T=\frac{{\beta }_2|T_1|}{{\beta }_2|T_1|+{\beta }_1|T_2|},{\kappa }_2{|}_T=\frac{{\beta }_1|T_2|}{{\beta }_2|T_1|+{\beta }_1|T_2|},$$

$$\{\{w\}\}={\kappa }_1{v}_{1,h}+{\kappa }_2{v}_{2,h},\{\{w\}{\}}^{\ast }={\kappa }_2{v}_{1,h}+{\kappa }_1{v}_{2,h}.$$

$$a_h(u_h,v_h)=L_h(v_h),\forall v_h\in V_{h,0};$$

$$\begin{array}{rl}{a}_h(u_h,v_h)=& \sum _{i=1}^2{\left(\beta \nabla u_{i,h},\nabla v_{i,h}\right)}_{{\Omega }_i}-{⟨[[u_h]],\{\{\beta {\partial }_n{v}_h\}\}⟩}_{\Gamma }\\ & -{⟨[[v_h]],\{\{\beta {\partial }_n{u}_h\}\}⟩}_{\Gamma }+h^{-1}{⟨\gamma [[u_h]],[[v_h]]⟩}_{\Gamma },\end{array}$$

$$L_h(v_h)=\sum _{i=1}^2(f,\nabla v_{i,h}{)}_{{\Omega }_i}-⟨q,\{\{\beta {\partial }_n{v}_h\}\}{⟩}_{\Gamma }+{⟨\gamma q,[[v_h]],⟩}_{\Gamma }+{⟨g,\{\{w\}{\}}^{\ast }⟩}_{\Gamma },$$

$$\gamma {|}_T=\frac{2h_T|{\Gamma }_T|}{|T_1|/{\beta }_1+|T_2|/{\beta }_2}.$$

$$a_h(\tilde{u},v_h)=L(v_h),\forall v_h\in V_{h,0},$$

$$a_h(\tilde{u}-u_h,v_h)=0,\forall v_h\in V_{h,0},$$

$$|||v||{|}_h^2=\Vert \nabla u{\Vert }_{0,{\Omega }_1\cup {\Omega }_2}^2+\sum _{T\in {\mathcal{T}}_{\Gamma ,h}}{h}_T\Vert \{\{{\partial }_{n}v\}\}{\Vert }_{0,{\Gamma }_T}^2+\sum _{T\in {\mathcal{T}}_{\Gamma ,h}}{h}_T^{-1}\Vert [[v]]{\Vert }_{0,{\Gamma }_T}^2.$$

$$C|||v_h||{|}_h^2\le a_h(v_h,v_h),\forall v_h\in V_h.$$

$$|||\tilde{u}-u_h||{|}_h\le Ch\Vert u{\Vert }_{2,{\Omega }_1\cup {\Omega }_2},$$

$$||\tilde{u}-u_h|{|}_{0,\Omega }\le C{h}^2\Vert u{\Vert }_{2,{\Omega }_1\cup {\Omega }_2},$$

$$\sum _{T\in {\mathcal{T}}_h^2}|T|=\mathcal{O}(h^{\sigma }).$$

$$(E_{i}w){|}_{{\Omega }_i}=w,$$

$$\Vert E_{i}w{\Vert }_{s,\Omega }\le C\Vert w{\Vert }_{s,{\Omega }_i},\forall w\in H^s({\Omega }_i),s=0,1,2,3.$$

$$I_h^{\ast }v=(I_{1,h}^{\ast }{v}_1,I_{1,h}^{\ast }{v}_2),$$

Full text

1 Introduction

Elliptic interface problems arise in many applications such as fluid dynamics and materials science, where the background consists of rather different materials on the subdomains separated by smooth curves called interface. The numerical challenge of interface problems comes from the fact that the solution, in general, has low global regularity due to the discontinuity of parameter (e.g. dielectric constant) at the interface. Standard finite element methods have been studied for elliptic interface problems by aligning the triangulation along the interface (body-fitted meshes), and are proven to achieve optimal convergence rates in both $L_{2}$ and energy norms [8, 3, 14, 54]. However, when the interface leads to subdomains of complex geometry, it is non-trivial and time-consuming to generate body-fitted meshes.

To overcome the difficulty of mesh generation in standard finite element method, tremendous effort has been input to develop numerical methods using unfitted (Cartesian) meshes. [49] is the first to propose the immersed boundary method (IBM) to simulate blood flow using Cartesian meshes. The idea of IBM is to use a Dirac $\delta$-function to model the discontinuity and discretize it to distribute a singular source to nearest grid points [49, 50]. But IBM only achieves the first-order accuracy. To improve the accuracy, Leveque and Li characterized the discontinuity as jump conditions and proposed the immersed interface method (IIM) [35]. IIM constructs special finite difference schemes to incorporate the jump conditions near the interface. High order unfitted finite difference methods including matched interface and boundary (MIB) method are also proposed in [60, 59]. We refer to [38] for a review on IIM and other unfitted finite difference methods.

In the meantime, unfitted numerical methods using finite element formulation are also developed for elliptic interface problems. The extended finite element method [7, 43, 19] enriches the standard continuous finite element space by adding some special basis functions to capture the discontinuity. The immersed finite element methods [36, 39, 40] modify the basis functions to satisfy the homogeneous jump conditions on interface elements. For the nonconforming immersed finite element method(IFEM) in [39], the numerical solution is continuous inside each element but can be discontinuous on the boundary of each element. Recently, there are also improved versions of IFEM such as the Petrov-Galerkin IFEM [30, 33, 31], symmetric and consistent IFEM [34], and partially penalized IFEM [42].

The Nitsche’s method [5, 11, 12, 10, 28, 25, 26, 27, 29], also called the cut finite element method, is firstly proposed by Hansbo and Hansbo in [25] to solve elliptic interface problems using unfitted meshes. It is further extended to deal with elastic problems with strong and weak discontinuities [26]. The study of the Nitsche’s method for Stokes interface problems can be found in [29]. The key idea of the Nitsche’s method is to construct an approximate solution on each fictitious domain and use Nitsche’s technique [48] to patch them together. A similar idea was used to develop the fictitious domain method [10, 11]. The robust forms of the unfitted Nitsche’s method were given in [2, 52]. The recent development of the cut finite element method is referred to the review paper [12].

For elliptic interface problems, computation of gradient plays an important role in many practical problems as discussed in [41], which demands numerical methods of high order accuracy. For standard elliptic problems, it is well known that the gradient recovery techniques [61, 63, 62, 57, 23, 1, 47, 46, 13, 55, 4] can reconstruct a highly accurate approximate gradient from the primarily computed data with reasonable cost. But for elliptic interface problems, only a few works have been done on the gradient recovery and associated superconvergence theory. For example in [53], a supercloseness result between the gradient of the linear finite element solution and the gradient of the linear interpolation is proved for a two-dimensional interface problem with a body-fitted mesh. For IFEM, Chou et al. introduced two special interpolation formulae to recover flux with high order accuracy for the one-dimensional linear and quadratic IFEM [15, 16]. Moreover, Li and his collaborators recently proposed an augmented immersed interface method [41] and a new finite element method [51] to accurately compute the gradient of the solution to elliptic interface problems. In our recent work [20], we proposed an improved polynomial preserving recovery for elliptic interface problems based on a body-fitted mesh and proved the superconvergence on both mildly unstructured meshes and adaptively refined meshes. Later in the two-dimensional case [21], we proposed gradient recovery methods based on symmetric and consistent IFEM [34] and Petrov-Galerkin IFEM [30, 32, 33] and numerically verified its superconvergence. In [24], we also provided a supercloseness result for the partially penalized IFEMs and proved that the recovered gradient using the gradient recovery method in [21] is superconvergent to the exact gradient.

Despite the efficiency and accuracy that the methods mentioned above bring to recover the gradient, there are still some cases in unfitted meshes where skinny triangles appear in the generated local body-fitted triangulation that destroy the accuracy of recovered gradient near the interface. In this paper, we propose an unfitted polynomial preserving recovery (UPPR) based on the Nitsche’s method. The key idea is to decompose the domain into two overlapping subdomains, named fictitious domains, by the interface and the triangulation proposed in the Nitsche’s method. On each fictitious domain, the standard linear finite element space will be used, and thus the classical polynomial preserving recovery (PPR) can be applied in each fictitious domain. Compared to previous gradient recovery methods [20, 21], the new method does not require generating a local body-fitted mesh and therefore avoids the drawback caused by skinny triangles. In general, the exact solutions of the interface problems are piecewise smooth on each subdomain. It implies that the extension of the exact solution on each subdomain to the whole domain is smooth, based on which, the recovered gradient using the interpolation of the exact solution is proven to be superconvergent to the exact gradient at rate of $\mathcal{O}(h^{2})$, and this is similar to the classical PPR for standard elliptic problems. In addition, we prove $\mathcal{O}(h^{1.5})$ supercloseness between the gradient given by the Nitsche’s method and the gradient of the interpolation of the exact solution by a sharp argument. This enables us to establish the complete superconvergence theory for the proposed UPPR.

The rest of the paper is organized as follows: We introduce briefly the elliptic interface problem and the Nitsche’s method in Section 2. In Section 3, we analyze the supercloseness for the Nitsche’s method and prove the $\mathcal{O}(h^{1.5})$ supercloseness between the gradient of the finite element solution and the gradient of the interpolation of the exact solution. In section 4, we describe the UPPR for the Nitsche’s method and establish its superconvergence theory. In Section 5, we present several numerical examples to confirm our theoretical results.

2 Nitsche’s method for elliptic interface problem

In this section, we first introduce the elliptic interface problem and associated notations, and then summarize the unfitted finite element discretization based on Nitsche’s method proposed in [2, 25] as a preparation for the unfitted polynomial preserving recovery (UPPR) method introduced later.

2.1 Elliptic interface problem

Let $\Omega$ be a bounded polygonal domain with Lipschitz boundary $\partial\Omega$ in $\mathbb{R}^{2}$. A $C^{2}$-curve $\Gamma$ divides $\Omega$ into two disjoint subdomains $\Omega_{1}$ and $\Omega_{2}$ as in Figure 1. We consider the following elliptic interface problem

[TABLE]

where $\partial_{n}u=(\nabla u)\cdot n$ with $n$ being the unit outward normal vector of $\Gamma$ and the jump $\llbracket w\rrbracket$ on $\Gamma$ is defined as

[TABLE]

with $w_{i}=w|_{\Omega_{i}}$ being the restriction of $w$ on $\Omega_{i}$. The diffusion coefficient $\beta(z)\geq\beta_{0}$ is a piecewise smooth function, i.e.

[TABLE]

which has a finite jump of function value at the interface $\Gamma$.

In this paper, we use the standard notations for Sobolev spaces and their associated norms as in [9, 17, 18]. For any bounded domain $D\subset\Omega$, the Sobolev space with norm $\|\cdot\|_{k,p,D}$ and seminorm $|\cdot|_{k,p,D}$ is denoted by $W^{k,p}(D)$. When $p=2$, $W^{k,2}(D)$ is simply denoted by $H^{k}(D)$ and the subscript $p$ is omitted in its associate norm and seminorm. Similar notations are applied to subdomains of $\Gamma$. Let $(\cdot,\cdot)_{D}$ and $\langle\cdot,\cdot\rangle_{\Gamma}$ denote the standard $L_{2}$ inner products of $L_{2}(D)$ and $L_{2}(\Gamma)$, respectively. For a bounded domain $D=D_{1}\cup D_{2}$ with $D_{1}\cap D_{2}=\emptyset$, let $W^{k,p}(D_{1}\cup D_{2})$ be the function space consisting of piecewise Sobolev functions $w$ such that $w|_{D_{1}}\in W^{k,p}(D_{1})$ and $w|_{D_{2}}\in W^{k,p}(D_{2})$, whose norm is defined as

[TABLE]

and seminorm is defined as

[TABLE]

In this paper, we denote $C$ as a generic positive constant which can be different at different occurrences. In addition, it is independent of mesh size and the location of the interface.

2.2 Nitsche’s method

Let $\mathcal{T}_{h}$ be a triangulation of $\Omega$ independent of the location of the interface $\Gamma$. For any element $T\in\mathcal{T}_{h}$, let $h_{T}$ be the diameter of $T$ and $\rho_{T}$ be the diameter of the circle inscribed in $T$. In addition, we make the following assumptions on the triangulation.

Assumption 2.1**.**

The triangulation $\mathcal{T}_{h}$ is shape regular in the sense that there is a constant $\sigma$ such that

[TABLE]

for any $T\in\mathcal{T}_{h}$.

Assumption 2.2**.**

The interface $\Gamma$ intersects each interface element boundary $\partial T$ exactly twice, and each open edge at most once.

To define the finite element space, denote the set of all elements that intersect the interface $\Gamma$ by

[TABLE]

and denote the union of all such type elements by

[TABLE]

Denote the set of all elements covering subdomain $\Omega_{i}$ to be

[TABLE]

and let

[TABLE]

Figure 2 gives an illustration of $\Omega_{i,h}$ and $\omega_{i,h}$. We remark that $\Omega_{1,h}$ and $\Omega_{2,h}$ overlap on $\Omega_{\Gamma,h}$, which is shown as the shaded part in Figures 2b and 2c.

Let $V_{i,h}$ be the standard continuous linear finite element space on $\Omega_{i,h}$, i.e.

[TABLE]

where $\mathbb{P}_{1}(T)$ is the space of polynomials with degree $\leq 1$ on $T$. Then, we define the finite element space $V_{h}$ as

[TABLE]

and $V_{h,0}$ as

[TABLE]

Note that a function in $V_{h}$ is a vector-valued function from $\mathbb{R}^{2}\mapsto\mathbb{R}^{2}$, which has a zero component in $\omega_{1,h}\bigcup\omega_{2,h}$ but in general two non-zero components in $\mathcal{T}_{\Gamma,h}$. It means that one will have two sets of basis functions for any element $T$ in $\mathcal{T}_{\Gamma,h}$: one for $V_{1,h}$ and the other for $V_{2,h}$.

For any $T$ in $\mathcal{T}_{\Gamma,h}$, denote $T_{i}=T\cap\Omega_{i}$ as the part of $T$ in $\Omega_{i}$, with $|T_{i}|$ being the measure of $T_{i}$ in $\mathbb{R}^{2}$. Denote $\Gamma_{T}=\Gamma\cap T$ as the part of $\Gamma$ in $T$, with $|\Gamma_{T}|$ being the measure of $\Gamma_{T}$ in $\mathbb{R}^{1}$. To increase the robustness of the Nitsche’s method, we introduce two weights as used in [2]

[TABLE]

which satisfies that $\kappa_{1}+\kappa_{2}=1$. Then, we define the weighted averaging of a function $v_{h}$ in $V_{h}$ on the interface $\Gamma$ as

[TABLE]

Using the notations introduced above, the Nitsche’s method [12, 25, 2] for the elliptic interface problem (2.1) is to find $u_{h}\in V_{h,0}$ such that

[TABLE]

where the bilinear form $a_{h}$ is defined as

[TABLE]

and the linear functional $L_{h}$ is defined as

[TABLE]

with the stability parameter

[TABLE]

In [2], the discrete variational form is shown to be consistent as the following theorem:

Theorem 2.3**.**

Let $u$ be the solution of the interface problem (2.1). Then we have

[TABLE]

where $\tilde{u}=(u_{1},u_{2})$ with the short-hand notation $u_{i}=u|_{\Omega_{i}}$ as used in (2.2).

Theorem 2.3 implies the following Galerkin orthogonality:

Corollary 2.4**.**

Let $u$ be the solution of (2.1) and $u_{h}$ be the solution of the discrete problem (2.16). Then we have

[TABLE]

with $\tilde{u}=(u_{1},u_{2})$.

To analyze the stability of the bilinear form $a_{h}(\cdot,\cdot)$, we introduce the following mesh-dependent norm [12, 25]

[TABLE]

In [25], it is shown that the bilinear form $a_{h}(\cdot,\cdot)$ is coercive with respect to the above mesh-dependent norm in the following sense

Theorem 2.5**.**

There is a constant C such that

[TABLE]

Based on the above coercivity, Hansbo et al. proved the following optimal convergence result [25]:

Theorem 2.6**.**

Let $u$ be the solution to the interface problem (2.1) and $u_{h}$ be the finite element solution to (2.16). If $u\in H^{2}(\Omega_{1}\cup\Omega_{2})$, then

[TABLE]

and

[TABLE]

with $\tilde{u}=(u_{1},u_{2})$.

3 Supercloseness Analysis

In this section, we establish the supercloseness result between the gradient of the finite element solution and the gradient of the interpolation of the exact solution as a preparation for the superconvergence anaysis of the proposed gradient recovery technique based on Nitsche’s method. For that propose, we need the triangulation $\mathcal{T}_{h}$ to satisfy Condition $(\sigma,\alpha)$ as explained below.

Two adjacent triangles are said to form an $\mathcal{O}(h^{1+\alpha})$ approximate parallelogram if the lengths of any two opposite edges differ only by $\mathcal{O}(h^{1+\alpha})$.

Definition 3.1**.**

The triangulation $\mathcal{T}_{h}$ is called to satisfy Condition $(\sigma,\alpha)$ if there exist a partition $\mathcal{T}_{h}^{1}\cup\mathcal{T}_{h}^{2}$ of $\mathcal{T}_{h}$ and positive constants $\alpha$ and $\sigma$ such that every two adjacent triangles in $\mathcal{T}_{h}^{1}$ form an $\mathcal{O}(h^{1+\alpha})$ parallelogram and

[TABLE]

Remark 3.1**.**

For the Nitsche’s method, we usually use a Cartesian mesh which is independent of the location of the interface. Therefore any Cartesian mesh satisfies Condition $(\sigma,\alpha)$ with $\sigma=\infty$ and $\alpha=1$.

To define the interpolation operator, we need to extend the function defined on the subdomain $\Omega_{i}$ to the whole domain $\Omega$. Let $E_{i}$, $i=1,2$, be the $H^{3}$-extension operator from $H^{3}(\Omega_{i})$ to $H^{3}(\Omega)$ such that

[TABLE]

and

[TABLE]

Let $I_{i,h}$ be the standard nodal interpolation operator from $C(\overline{\Omega})$ to $V_{i,h}$. Define the interpolation operator for the finite element space $V_{h}$ as

[TABLE]

where

[TABLE]

Optimal approximation capability of $I_{h}^{\ast}$ is proved in [25]. Assume $\mathcal{T}_{h}$ satisfies Condition $(\sigma,\alpha)$, and then we can prove the following theorem:

Theorem 3.2**.**

Suppose the triangulation $\mathcal{T}_{h}$ satisfies Condition $(\sigma,\alpha)$. Let $u$ be the solution of the interface problem (2.1) and $I_{h}^{\ast}u$ be the interpolation of $u$ in the finite element space $V_{h,0}$. If $u\in H^{1}(\Omega)\cap H^{3}(\Omega^{-}\cup\Omega^{+})\cap W^{2,\infty}(\Omega^{-}\cup\Omega^{+})$, then for $\tilde{u}=(u_{1},u_{2})$ and all $v_{h}\in V_{h,0}$,

[TABLE]

where $\rho=\min(\alpha,\frac{\sigma}{2},\frac{1}{2})$.

Proof.

By (2.17), we have

[TABLE]

Since $\mathcal{T}_{h}$ satisfies Condition $(\sigma,\alpha)$, it follows that $\mathcal{T}_{1,h}\setminus\mathcal{T}_{\Gamma,h}$ and $\mathcal{T}_{2,h}\setminus\mathcal{T}_{\Gamma,h}$also satisfy Condition $(\sigma,\alpha)$. Notice that the restriction of the finite element space $V_{i,h}$ on $\omega_{i,h}$ is just the standard continuous linear finite element space for $i=1,2$. Then by Lemma $2.1$ in [55], we have

[TABLE]

where $\rho=\min(\alpha,\frac{\sigma}{2},\frac{1}{2})$. To estimate $I_{3}$, the Cauchy-Schwartz inequality implies

[TABLE]

[TABLE]

where we have used the fact $\sum_{T\in\mathcal{T}_{\Gamma,h}}1\approx\mathcal{O}(h^{-1})$. By the Cauchy-Schwartz inequality and the trace inequality in [25, 52], we have

[TABLE]

Similarly, we can estimate $F_{5}$ and $F_{6}$ as

[TABLE]

Combing all the above estimations, we get (3.5). ∎

Now we state our main supercloseness result as follows.

Theorem 3.3**.**

Assume the same hypothesis as in Theorem 3.2 and let $u_{h}$ be the finite element solution of the discrete variational problem (2.16), then

[TABLE]

where $\rho=\min(\alpha,\frac{\sigma}{2},\frac{1}{2})$.

Proof.

By Corollary 2.4, we have the following Galerkin orthogonality

[TABLE]

where $\tilde{u}=(u_{1},u_{2})$. Then we have

[TABLE]

Taking $v_{h}=u_{h}-I_{h}^{\ast}u$ and using Theorem 3.2 and Theorem 2.5, we prove (3.11). ∎

Remark 3.2**.**

Theorem 3.3 implies that we can have the same supercloseness result as the partially penalized IFEM [24] .

4 Superconvergent gradient recovery

In this section, we first propose the unfitted polynomial preserving recovery (UPPR) technique based on the Nitsche’s method, then prove that the recovered gradient by UPPR is superconvergent to the exact gradient on mildly unstructured meshes.

4.1 Unfitted polynomial preserving recovery

To accurately recover the gradient, we notice that the finite element solution $u_{h}$ of the Nitsche’s method (2.16) consists of two part: $u_{1,h}$ and $u_{2,h}$. Also, by the fact that we can smoothly extend the exact solution $u|_{\Omega_{i}}$ ($i=1,2$) to the whole domain, it is safe to assume $u|_{\Omega_{i}}$ and its extension $E_{i}u|_{\Omega_{i}}$ is smooth in general. For each $i\in\left\{1,2\right\}$, $u_{i,h}\in V_{h}$ is a continuous piecewise polynomial on fictitious domain $\Omega_{i,h}$ but its gradient $\nabla u_{i,h}$ is only a piecewise constant function. This motivates us to naturally use some smoothing operators such as superconvergent patch recovery (SPR) and polynomial preserving recovery (PPR) to smooth the discontinuous gradient into a continuous one on each fictitious domain $\Omega_{i,h}$.

To this end, let $G_{h}^{i}$ be the PPR gradient recovery operator [57, 47] on the fictitious domain $\Omega_{i,h}$ for $i=1,2$. Then $G_{h}^{i}$ is a linear operator from $V_{i,h}$ to $V_{i,h}\times V_{i,h}$ whose value at each nodal point is obtained by the local least squares fitting using sampling points only located in $\Omega_{i,h}$. According to [22, 57, 47], the gradient recovery operator $G_{h}^{i}$ is bounded in the sense that

[TABLE]

and is consistent in following sense that

[TABLE]

for $i=1,2$.

Let $u_{h}$ be the finite element solution of the discrete variational problem (2.16). We define the recovered gradient of $u_{h}$ as

[TABLE]

The linearity of $G_{h}^{i}$ implies $R_{h}$ is a linear operator from $V_{h}$ to $V_{h}\times V_{h}$. $R_{h}$ is called the unfitted polynomial preserving recovery (UPPR).

Remark 4.1**.**

The definition of the gradient recovery operator can be presented in a more general form. In fact, $G_{h}^{i}$ can be chosen as any local gradient recovery operators [56] like simple averaging, weight averaging, SPR and PPR. For simplicity and efficiency, we only consider $G_{h}^{i}$ as PPR here.

Remark 4.2**.**

The main idea is to use the standard PPR on each fictitious domain $\Omega_{i,h}$, which is similar to the improved polynomial preserving recovery for the finite element method based on body-fitted meshes [20]. But the proposed method (4.3) does not require the mesh fitting the interface, and that is why we call it unfitted polynomial preserving recovery method.

Remark 4.3**.**

We considered the gradient recovery technique for immersed finite element methods in [21], which is also based on unfitted meshes. The gradient recovery technique in [21] needs to generate a local body-fitted mesh by dividing every interface triangle into three sub-triangles, which can lead to skinny triangles and therefore a loss of accuracy. The proposed gradient recovery operator (4.3) overcomes this drawback.

Note that as a function in $V_{h}\times V_{h}$, $R_{h}u_{h}$ is continuous on each subdomain $\Omega_{i}$ and is discontinuous in the whole domain $\Omega$ which approximates the exact gradient $\nabla u$. Also, similar to the finite element solution $u_{h}$, both $G_{h}^{1}u_{1,h}$ and $G_{h}^{2}u_{2,h}$ in (4.3) are, in general, non-zero on interface triangles $T\in\mathcal{T}_{h,\Gamma}$. For the gradient recovery operator $R_{h}$ (4.3), we can show that it is consistent as follows:

Theorem 4.1**.**

Let $R_{h}:V_{h}\rightarrow V_{h}\times V_{h}$ be the unfitted polynomial preserving recovery operator defined in (4.3) and $I_{h}^{\ast}$ be the interpolation of $u$ into the finite element space $V_{h}$ as defined in (3.3). If $u\in H^{3}(\Omega_{1}\cup\Omega_{2})$, then we have

[TABLE]

with $\tilde{\nabla}u=(\nabla u_{1},\nabla u_{2})$.

Proof.

By (2.4), (3.3), (3.4), and (4.2), we have

[TABLE]

Taking square root on both sides of (4.5) completes our proof. ∎

Remark 4.4**.**

Theorem 4.1 means the recovered gradient using the interpolation of the exact solution is superconvergent to the exact gradient at a rate of $\mathcal{O}(h^{2})$. It is similar to the classical PPR operator for regular elliptic problems.

4.2 Superconvergence analysis

In the following, we shall show the superconvergence property of the proposed UPPR. Our main superconvergent tool is the supercloseness result provided in Section 3.

Theorem 4.2**.**

Under the same hypothesis as in Theorem 3.2. We further assume that $u_{h}$ is the finite element solution of the discrete variational problem (2.16). Then we have

[TABLE]

where $\tilde{\nabla}u=(\nabla u_{1},\nabla u_{2})$ and $\rho=\min(\alpha,\frac{\sigma}{2},\frac{1}{2})$.

Proof.

By the triangle inequality, we have

[TABLE]

According to Theorem 4.1, we have

[TABLE]

For $F_{2}$, we have

[TABLE]

[TABLE]

Theorem 3.3 implies that

[TABLE]

Then, we estimate $F_{4}$ as

[TABLE]

where we have used the fact $\sum_{T\in\mathcal{T}_{\Gamma,h}}1\approx\mathcal{O}(h^{-1})$. Similarly, we have

[TABLE]

Combining the estimates for $F_{3}$, $F_{4}$, and $F_{5}$, we have

[TABLE]

which completes the proof. ∎

By the above superconvergence result, we naturally define a local a posteriori error estimator on an element $T\in\mathcal{T}_{h}$ :

[TABLE]

and the corresponding global error estimator

[TABLE]

Theorem 4.2 implies the error estimator (4.6) (or (4.7)) is asymptotically exact for the Nitsche’s method:

Theorem 4.3**.**

Assume the same hypothesis in Theorem 3.2 and let $u_{h}$ be the finite element solution of the discrete variational problem (2.16). Further assume that there is a constant $C(u)>0$ such that

[TABLE]

then it holds that

[TABLE]

where $\tilde{\nabla}u=(\nabla u_{1},\nabla u_{2})$ and $\rho=\min(\alpha,\frac{\sigma}{2},\frac{1}{2})$.

Proof.

By Theorem 4.2 and (4.8), we have

[TABLE]

∎

Remark 4.5**.**

For interface problems, there are two types of errors: the error introduced by geometric discretization and the error introduced by the singularity of the solution. The first type of error can be predicted by the curvature of the interface [53, 58]. The error estimator (4.6) or (4.7) can be used to estimate the second type of error.

5 Numerical examples

In this section, we show the performance of proposed unfitted polynomial preserving recovery (UPPR) method by several numerical examples with both simple and complex interface geometries. The computational domains of all examples are chosen as $\Omega=(-1,1)\times(-1,1)$. For the first two numerical examples, the uniform triangulations of $\Omega$ are obtained by dividing $\Omega$ into $N^{2}$ sub-squares and then dividing each sub-square into two right triangles. The resulted uniform mesh size is $h=\frac{2}{N}$. For convenience, we use the following errors in all the examples:

[TABLE]

with $\tilde{\nabla}u=(\nabla u_{1},\nabla u_{2})$.

Example 5.1. In this example, we consider the interface problem (2.1) with homogeneous jump condition as in [39]. The interface is a circular interface of radius $r_{0}=0.5$. The exact solution is

[TABLE]

where $r=\sqrt{x^{2}+y^{2}}$.

We consider the following four typical different jump ratios: $\beta_{1}/\beta_{2}=1/10$ (moderate jump), $\beta_{1}/\beta_{2}=1/1000$ (large jump), $\beta_{1}/\beta_{2}=1/100000$ (huge jump), and $\beta_{1}/\beta_{2}=100000$ (huge jump). The numerical errors are displayed in Tables 1-4. We observe an optimal convergence in the $H^{1}$-seminorm as predicted by Theorem 2.6. The observed $\mathcal{O}(h^{1.5})$ supercloseness and superconvergence confirm our theoretical results. In addition, we observe the same superconvergence results in all different cases. It means that the superconvergence results are independent of the jump ratio of the coefficient. In Figure 3, we plot the recovered gradient on the initial mesh.

Example 5.2. In this example, we consider the flower-shape interface problem with non-homogeneous jump conditions as studied in [44, 59]. The interface curve $\Gamma$ in polar coordinates is given by

[TABLE]

It contains both convex and concave parts, as demonstrated in Figure 4a. The diffusion coefficient is piecewise constant with $\beta_{1}=1$ and $\beta_{2}=10$. The right-hand side function $f$ in (2.1a) is chosen to match the exact solution

[TABLE]

and the jump conditions (2.1c)-(2.1d) are provided by the exact solution.

In Figure 4b, we plot the numerical solution on the initial mesh which clearly indicates the non-homogeneous jump in function value. We show the numerical results in Table 5. As expected, we observe the first-order convergence for the gradient of finite element solution. For the recovered gradient, $\mathcal{O}(h^{1.5})$ convergence is observed, which is in agreement with Theorem 4.2. The recovered gradient on the initial mesh is visualized in Figure 5.

Example 5.3. In this example, we consider the interface problem with complex geometrical structure as in [45]. The interface in polar coordinates is given by

[TABLE]

The interface and subdomains are plotted in Figure 4a. The coefficient function is

[TABLE]

and the exact function is

[TABLE]

As plotted in Figure 6a, the interface contains complex geometrical structure. To guarantee the Assumption 2.2, we need an extremely fine mesh. It would increase the computational cost. To reduce the computational cost, we propose an adaptive strategy to generate an initial unfitted mesh. Here we use the curvature-based a posterior estimator to guide the refinement of the mesh as in [53]. Different from the mesh generated in [53], the resulted mesh is an unfitted mesh and all triangles are perfect right triangles.

Figure 6b plots the generated initial unfitted mesh. It is easy to see that the mesh is refined around the part of the interface with high curvature. The other four levels of unfitted meshes are obtained by uniform refinement. The numerical results are summarized in Table 6. Note that in Table 6, convergence rates are listed with respect to the degree of freedom (DOF). The corresponding convergent rates with respect to the mesh size $h$ are double of what we present in Table 6. The gradient of finite element solution converges to the exact gradient at the rate of $\mathcal{O}(h)$ while the recovered gradient superconverges at the rate of $\mathcal{O}(h^{1.5})$. Additionally, the predicted supercloseness is observed in the numerical experiment.

Example 5.4. In this example, we consider the interface problem as in [6, 53]. The interface in parametric form is given by

[TABLE]

where

[TABLE]

The coefficient function is

[TABLE]

and the exact solution is

[TABLE]

The interface $\Gamma$, shown in Figure 7a, contains complex geometrical structure. We use the same algorithm as in Example 5.3 to generate an initial unfitted mesh which is plotted in Figure 7b. Table 7 lists the numerical results. Clearly, we observe the desired optimal convergence and superconvergence rates.

Example 5.5. In this example, we consider the interface problem as in [37, 53]. The interface $\Gamma$ in parametric form is defined by

[TABLE]

where $r(\theta)=r_{0}+r_{1}\sin(\omega\theta)$, $0\leq\theta<2\pi$.

In this test, we take $r_{0}=0.4$, $r_{1}=0.2$, $\omega=20$, and $x_{c}=y_{c}=0.02\sqrt{5}$. The coefficient $\beta$ is a piecewise constant with $\beta_{1}=1$ and $\beta_{2}=10$. The exact function is

[TABLE]

The interface $\Gamma$ is plotted in Figure 8a and the adaptively refined initial mesh is shown in Figure 8b. The numerical results are given in Table 8. The observed results confirm the first-order convergence rate as predicted by Theorem 2.6. For the errors $De^{i}$ and $De^{r}$, $\mathcal{O}(h^{1.8})$ order decaying rates can be observed which are better than our theoretical results. Compared to the numerical results using a body-fitted mesh in [53], we achieve the same accuracy by using an unfitted mesh with about one sixth of the total mesh grid points.

6 Conclusion

In this paper, we propose a new gradient recovery technique based on the Nitsche’s method. Compared to our previous works [20, 21, 24], it avoids the loss of accuracy of gradient near the interface caused by skinny triangles. By proving the supercloseness result for the Nitsche’s method, we are able to show that the recovered gradient is superconvergent to the exact gradient. As a byproduct, we propose a curvature estimator based adaptive algorithm to generate initial unfitted triangulations for the elliptic interface problems with complex geometry, which greatly reduces the computational cost as illustrated in Examples 5.3, 5.4 and 5.5. The future work is planned in several different directions: firstly, we will extend the study into three dimension problems; secondly, we will consider other type equations like elastic interface problems and wave propagation problems in heterogeneous media; thirdly, we will combine the curvature estimator and the recovery-based a posterior error estimator to derive adaptive algorithms for the complex interface problems.

Acknowledgement

This work was partially supported by the NSF grant DMS-1418936 and DMS-1107291.