A Discontinuous Galerkin Method by Patch Reconstruction for Elliptic Interface Problem on Unfitted Mesh

TL;DR

This paper introduces a novel discontinuous Galerkin method with patch reconstruction for elliptic interface problems on unfitted meshes, achieving optimal error estimates and improved efficiency over conforming finite element methods.

Contribution

The paper presents a new DG method using patch reconstruction with one degree of freedom per element, removing typical interface intersection constraints and enhancing efficiency.

Findings

Achieves optimal error estimates in L2 and DG energy norms.

Does not require constraints on interface intersection with mesh elements.

Numerical examples demonstrate improved efficiency over conforming finite element methods.

Abstract

We propose a discontinuous Galerkin(DG) method to approximate the elliptic interface problem on unfitted mesh using a new approximation space. The approximation space is constructed by patch reconstruction with one degree of freedom per element. The optimal error estimates in both L2 norm and DG energy norm are obtained, without the typical constraints for DG method on how the interface intersects to the elements in the mesh. Other than enjoying the advantages of DG method, our method may achieve even better efficiency than the conforming finite element method, as illustrated by numerical examples.

Figures

Tables

Table 1. Table 1. The uniform # ​ S i ​ ( K ) # superscript 𝑆 𝑖 𝐾 \#S^{i}(K) for 1 ≤ m ≤ 3 1 𝑚 3 1\leq m\leq 3 .

	$m$	1	2	3
$\#S^i\left(K\right)$	$d=2$	5	9	15
	$d=3$	9	18	38

Table 2. Table 2. The convergence orders under L 2 superscript 𝐿 2 L^{2} norm and DG energy norm for Example 6.

order $m$	$h$	$L^2$ error	order	DG error	order
$m=1$	2.00e-1	7.661e-3	-	1.835e-1	-
	1.00e-1	2.515e-3	1.61	5.022e-2	1.00
	5.00e-2	6.498e-4	1.95	2.445e-2	1.03
	2.50e-2	1.653e-4	1.97	1.199e-2	1.02
	1.25e-2	4.202e-5	1.98	1.156e-2	1.00
$m=2$	2.00e-1	4.727e-4	-	9.283e-3	-
	1.00e-1	6.423e-5	2.85	2.393e-3	1.95
	5.00e-2	7.249e-6	3.16	5.872e-4	2.01
	2.50e-2	9.171e-7	2.98	1.505e-4	1.97
	1.25e-2	1.126e-7	3.02	6.401e-5	2.02
$m=3$	2.00e-1	1.229e-4	-	3.145e-3	-
	1.00e-1	1.126e-5	3.45	3.361e-4	2.41
	5.00e-2	9.603e-7	3.55	5.721e-5	2.53
	2.50e-2	8.249e-8	3.55	1.816e-5	2.55
	1.25e-2	6.999e-9	3.56	3.108e-6	2.55

Table 3. Table 3. The L 2 superscript 𝐿 2 L^{2} errors ‖ u − u h 1 ‖ L 2 ​ ( Ω 0 ∪ Ω 1 ) subscript norm 𝑢 superscript subscript 𝑢 ℎ 1 superscript 𝐿 2 subscript Ω 0 subscript Ω 1 \|u-u_{h}^{1}\|_{L^{2}(\Omega_{0}\cup\Omega_{1})} and ‖ u − u h 2 ‖ L 2 ​ ( Ω ~ 0 ∪ Ω ~ 1 ) subscript norm 𝑢 superscript subscript 𝑢 ℎ 2 superscript 𝐿 2 subscript ~ Ω 0 subscript ~ Ω 1 \|u-u_{h}^{2}\|_{L^{2}(\widetilde{\Omega}_{0}\cup\widetilde{\Omega}_{1})} .

		$h=1/4$	$h=1/8$	$h=1/16$	$h=1/32$
$m=1$	${\Vert u-u_h^1\Vert }_{L^2\left({\Omega }_0\cup {\Omega }_1\right)}$	7.1514e-2	3.0977e-2	1.0696e-2	2.8099e-3
	${\Vert u-u_h^2\Vert }_{L^2\left({\tilde{\Omega }}_0\cup {\tilde{\Omega }}_1\right)}$	1.0216e-1	3.3384e-2	1.1173e-2	2.8274e-3
$m=2$	${\Vert u-u_h^1\Vert }_{L^2\left({\Omega }_0\cup {\Omega }_1\right)}$	6.1523e-2	2.0271e-3	1.8898e-4	2.3557e-5
	${\Vert u-u_h^2\Vert }_{L^2\left({\tilde{\Omega }}_0\cup {\tilde{\Omega }}_1\right)}$	3.4460e-2	2.4149e-3	2.1523e-4	2.3799e-5
$m=3$	${\Vert u-u_h^1\Vert }_{L^2\left({\Omega }_0\cup {\Omega }_1\right)}$	5.1865e-2	3.9632e-4	1.7698e-5	9.0215e-7
	${\Vert u-u_h^2\Vert }_{L^2\left({\tilde{\Omega }}_0\cup {\tilde{\Omega }}_1\right)}$	1.7197e-2	3.7255e-4	1.5553e-5	8.7805e-7

Equations

$$h_K=\text{diam}(K),\forall K\in {\mathcal{T}}_h,h_e=|e|,\forall e\in {\mathcal{E}}_h,$$

$${\rho }_s\underset{\tilde{K}\in \Delta (K)}{\max }{h}_{\tilde{K}}\le h_K,$$

$$\Vert v{\Vert }_{L^2(\partial K)}^2\le C\left(h_K^{-1}\Vert v{\Vert }_{L^2(K)}^2+h_K\Vert \nabla v{\Vert }_{L^2(K)}^2\right),\forall v\in H^1(K).$$

$$\Vert \nabla v{\Vert }_{L^2(K)}$$

$$e^0$$

$$K^0$$

$$(\partial K{)}^0$$

$${\mathcal{T}}_h^0$$

$${\mathcal{E}}_h^0$$

$$S_t^i(K)=\bigcup _{\begin{array}{rl}\tilde{K}\in {\mathcal{T}}_h^i,& \hat{K}\in S_{t-1}^i(K)\\ \tilde{K}\cap \hat{K}& =e\in {\mathcal{E}}_h\end{array}}\tilde{K},t=1,2,\cdots$$

$${\mathcal{I}}_K^i=\left\{{\mathbf{x}}_{\tilde{K}}|\forall \tilde{K}\in S^i(K)\right\}.$$

$${\mathcal{R}}_K^{i}g={\mathrm{argmin}}_{p\in {\mathbb{P}}_m(S^i(K))}\sum _{\mathbf{x}\in {\mathcal{I}}_K^i}|p(\mathbf{x})-g(\mathbf{x}){|}^2.$$

$$p{|}_{{\mathcal{I}}_K^i}=0\text{implies}p{|}_{S^i(K)}\equiv 0,i=0,1.$$

$$({\mathcal{R}}^{0}g){|}_K$$

$$({\mathcal{R}}^{1}g){|}_K$$

$$w_K^i(\mathbf{x})=\{\begin{array}{l}1,\mathbf{x}={\mathbf{x}}_K,\\ 0,\mathbf{x}\in \tilde{K},\tilde{K}\ne K.\end{array}$$

$${\mathcal{R}}^{i}g=\sum _{K\in {\mathcal{T}}_h^i}g({\mathbf{x}}_K){\lambda }_K^i(\mathbf{x}),\forall g\in C^0(\Omega ).$$

$$\Vert E^{i}w{\Vert }_{H^s(\Omega )}\le C\Vert w{\Vert }_{H^s({\Omega }_i)},0\le s\le t.$$

$$\Lambda (m,S^i(K))=\underset{p\in {\mathbb{P}}_m(S^i(K))}{\max }\frac{{\max }_{\mathbf{x}\in S^i(K)}|p(\mathbf{x})|}{{\max }_{\mathbf{x}\in {\mathcal{I}}_K^i}|p(\mathbf{x})|}.$$

$$\Vert E^{i}g-{\mathcal{R}}^i(E^{i}g){\Vert }_{H^q(K)}$$

$$\Vert {\nabla }^q(E^{i}g-{\mathcal{R}}^i(E^{i}g)){\Vert }_{L^2(\partial K)}$$

$$(\mathcal{R}w){|}_K\triangleq \{\begin{array}{l}({\mathcal{R}}_K^0{E}^{0}w){|}_K,\text{for }K\in {\mathcal{T}}_h^0\backslash {\mathcal{T}}_h^{\Gamma },\\ ({\mathcal{R}}_K^1{E}^{1}w){|}_K,\text{for }K\in {\mathcal{T}}_h^1\backslash {\mathcal{T}}_h^{\Gamma },\\ ({\mathcal{R}}_K^i{E}^{i}w){|}_{K^i},\text{for }K\in {\mathcal{T}}_h^{\Gamma },i=0,1.\end{array}$$

$$-\nabla \cdot \beta \nabla u$$

$$u$$

$$[[u]]$$

$$[[\beta \nabla u\cdot {\mathrm{n}}_{\Gamma }]]$$

$$[[\mathbf{q}]]=\{\begin{array}{l}\mathbf{q}{|}_{K_{+}}\cdot {\mathrm{n}}_{K_{+}}+\mathbf{q}{|}_{K_{-}}\cdot {\mathrm{n}}_{K_{-}}\\ \mathbf{q}{|}_{K_{+}^i}\cdot {\mathrm{n}}_{K_{+}}+\mathbf{q}{|}_{K_{-}^i}\cdot {\mathrm{n}}_{K_{-}}\\ \mathbf{q}{|}_K\cdot {\mathrm{n}}_K\\ (\mathbf{q}{|}_{K^1}-\mathbf{q}{|}_{K^0})\cdot {\mathrm{n}}_{\Gamma }\end{array}[[v]]$$

$$\{\mathbf{q}\}=\{\begin{array}{l}\frac{1}{2}(\mathbf{q}{|}_{K_{+}}+\mathbf{q}{|}_{K_{-}})\\ \frac{1}{2}(\mathbf{q}{|}_{K_{+}^i}+\mathbf{q}{|}_{K_{-}^i})\\ \mathbf{q}{|}_K\\ \frac{1}{2}(\mathbf{q}{|}_{K^1}+\mathbf{q}{|}_{K^0})\end{array}\{v\}$$

$$b_h(u_h,v_h)$$

$$\displaystyle-\left[\sum_{e\in\mathcal{E}_{h}}\int_{e^{0}\cup e^{1}}+\sum_{K\in\mathcal{T}_{h}^{\Gamma}}\int_{\Gamma_{K}}\right]\Big{(}[[u_{h}]]\cdot\{\beta\nabla v_{h}\}+[[v_{h}]]\cdot\{\beta\nabla u_{h}\}\Big{)}\mathrm{d}\boldsymbol{s}$$

$$+\sum _{e\in {\mathcal{E}}_h}{\int }_{e^0\cup e^1}\frac{\eta }{h_e}[[u_h]]\cdot [[v_h]]\mathrm{d}\mathbf{s}+\sum _{K\in {\mathcal{T}}_h^{\Gamma }}{\int }_{{\Gamma }_K}\frac{\eta }{h_K}[[u_h]]\cdot [[v_h]]\mathrm{d}\mathbf{s},$$

$$l_h(v_h)$$

$$+\sum _{K\in {\mathcal{T}}_h^{\Gamma }}{\int }_{{\Gamma }_K}b\{v_h\}\mathrm{d}\mathbf{s}-\sum _{K\in {\mathcal{T}}_h^{\Gamma }}{\int }_{{\Gamma }_K}a{\mathrm{n}}_{\Gamma }\cdot \{\beta \nabla v_h\}\mathrm{d}\mathbf{s}$$

$$+\sum _{e\in {\mathcal{E}}_h^b}{\int }_e\frac{\eta }{h_e}g{v}_h\mathrm{d}\mathbf{s}+\sum _{K\in {\mathcal{T}}_h^{\Gamma }}{\int }_{{\Gamma }_K}\frac{\eta }{h_K}a{\mathrm{n}}_{\Gamma }\cdot [[v_h]]\mathrm{d}\mathbf{s},$$

$$\displaystyle W_{h}=\Big{\{}v\in L^{2}(\Omega)\ \Big{|}$$

$$\displaystyle v|_{K^{i}}\in H^{2}(K^{i}),\ i=0,1,\ \text{for }K\in\mathcal{T}_{h}^{\Gamma}\Big{\}}.$$

$$|||v_h||{|}^2=\Vert \nabla v_h{\Vert }_{L^2({\mathcal{T}}_h^0\cup {\mathcal{T}}_h^1)}^2$$

Full text

The optimal error estimates in both $L^{2}$ norm and DG energy norm are obtained, without restrictions on how the interface intersects the elements in the mesh. The stability near the interface is ensured by the patch reconstruction and no special numerical flux is required. The convergence order by numerical results in both 2D and 3D agrees with the error estimates perfectly. More than enjoying the advantages of DG method, the new method may achieve even better efficiency in number of degree of freedom than the conforming finite element method as illustrated by our numerical examples.

keywords: Elliptic interface problems, Patch reconstructed, Discontinuous Galerkin method, Unfitted mesh.

1. Introduction

In the last decades, numerical methods for the elliptic interface problem have attracted pervasive attention since the pioneering work of Peskin [41], for example, the immersed interface method by LeVeque and Li [28, 33], Mayo’s method on irregular regions [38], the method in [53] with second-order accuracy in the $L^{\infty}$ norm. In the finite difference fold, we also refer to [35, 21, 22, 16, 11, 40] for some other interesting methods. Meanwhile, finite element (FE) method is also popular for solving the interface problem. Based on the geometrical relationship between the grid and the interface, FE methods could be classified into two categories: interface-fitted method and interface-unfitted method. The body-fitted grid enforces the mesh to align with the interface to render a high-order accurate approximation [12, 6]. However, generating a fitted mesh with satisfied quality is sometimes a nontrivial and time-consuming task [50, 51]. Therefore, there are some techniques for FE methods based on unfitted grid. The unfitted FE method can date back to the [5], which introduced a penalty term to weakly enforce the jump on the interface. Li proposed the immersed FE method in [32], which processes a better approximate solution by modifying the basis functions near interface to capture the jump of the solution. We refer to [3, 34, 2, 48, 10, 17] for some recent works. Let us note that the extended FE method is also a popular discretization method [7].

In 2002, A. Hansbo and P. Hansbo proposed an unfitted FE method with the piecewise linear space and proved an optimal order of convergence [20]. The numerical solution comes from separate solutions defined on each subdomain and the jump conditions are imposed weakly by Nitsche’s method. Wadbro et al. [46] developed a uniformly well-conditioned FE method based on Nitsche’s method. Wu and Xiao [23, 50] presented a $hp$ unfitted FE method, which is extended to the three dimensional case. To achieve high-order accuracy and enjoy additional flexibility, some authors tried to apply DG method to the elliptic interface problem, for example the local DG method in [18], the hybridizable DG method in [25] on fitted mesh, and the $hp$ DG method in [37] on unfitted mesh.

Though high-order accuracy can be obtained, solid difficulties remain for DG methods in solving problems with complex interfaces. To fit curved interfaces, Cangiani et al. [9] introduced elements with curved faces to give an adaptive DG method recently. As one of the latest work on unfitted mesh, Burman and Ern [8] proposed a hybrid high-order method, while an extra assumption on the meshes are required to ensure the mesh cells are cut favorably by the interface [37]. In this paper, we are trying to propose a DG method on unfitted mesh for the interface problem still using Nitsche’s method. The novel point is that we adopt a new approximation space by patch reconstruction with one degree of freedom (DOF) per element following the methodology in [30, 29]. The new space may be regarded as a subspace of the approximation space used in [37]. Thanks to the flexibility in choosing reconstruction patches, we may allow the interface to intersect elements in a very general manner, in comparison to the methods in [8, 37]. Following the standard DG discretization, the elliptic interface problem is approximated by using a symmetric interior penalty bilinear form with a Nitsche-type penalization at the interface. The optimal error estimate is then derived in both DG energy norm and $L^{2}$ norm. The patch reconstruction can provide the stability near the interface and no cut-dependent numerical flux is used. We note that the idea of using the patch of interface elements to improve the numerical stability can also be found in [19, 13, 23, 17]. In addition, the classical DG methods for elliptic problems were challenged [24, 55] since it may use more DOFs than traditional conforming FE methods. As a new observation, we demonstrate by numerical examples that using our new approximation space, one needs much less DOFs than classical DG methods. For high-order approximations, number of DOFs can be even less than conforming FE methods to achieve the same numerical error.

The rest of this paper is organized as follows. In Section 2, we introduce the reconstruction operator and the new approximation space, and we also give the basic properties of the approximation space. In Section 3, the approximation to the elliptic interface problem is proposed and we derive the optimal error estimate in DG energy norm and $L^{2}$ norm. In Section 4, we present a lot of numerical examples to verify the error estimate in Section 3. To show the performance of our method in efficiency, we make a comparison of number of DOFs with respect to the numerical error between different methods. We also solve a problem that admits solutions with low regularities to illustrate the robustness of our method.

2. Approximation Space

Let $\Omega\subset\mathbb{R}^{d}$, $d=2$ or $3$, be a convex and polygonal (polyhedral) domain with boundary $\partial\Omega$ and let $\Gamma$ be a $C^{2}$-smooth interface which divides $\Omega$ into two open sets $\Omega_{0}$ and $\Omega_{1}$ satisfying $\Omega_{0}\cap\Omega_{1}=\varnothing,\ \overline{\Omega}=\overline{\Omega}_{0}\cup\overline{\Omega}_{1}$ and $\Gamma=\overline{\Omega}_{0}\cap\overline{\Omega}_{1}$. We denote by $\mathcal{T}_{h}$ a partition of $\Omega$ into polygonal (polyhedral) elements. Here we do not require the faces of elements in $\mathcal{T}_{h}$ align with the interface (see Fig 1).

Let $\mathcal{E}^{\circ}_{h}$ be the set of all interior faces of $\mathcal{T}_{h}$, $\mathcal{E}^{b}_{h}$ the set of the faces on $\partial\Omega$ and then $\mathcal{E}_{h}=\mathcal{E}_{h}^{\circ}\cup\mathcal{E}_{h}^{b}$. We set

[TABLE]

and we denote by $h$ the biggest one among the diameters of all elements in $\mathcal{T}_{h}$. We assume that $\mathcal{T}_{h}$ is share-regular in the sense of satisfying the conditions introduced in [4], which are: there exist

• •

two positive numbers $N$ and $\sigma$ which are independent of mesh size $h$;

• •

a compatible sub-decomposition $\widetilde{\mathcal{T}}_{h}$ into shape-regular triangles (tetrahedrons);

such that

• •

any polygon (polyhedron) $K\in\mathcal{T}_{h}$ admits a decomposition $\widetilde{\mathcal{T}}_{h|K}$ which has less than $N$ shape-regular triangles (tetrahedrons);

• •

the share-regularity of $\widetilde{K}\in\widetilde{\mathcal{T}}_{h}$ follows [14]: the ratio between $h_{\widetilde{K}}$ and $\rho_{\widetilde{K}}$ is bounded by $\sigma$: $h_{\widetilde{K}}/\rho_{\widetilde{K}}\leq\sigma$ where $\rho_{\widetilde{K}}$ is the radius of the largest ball inscribed in $\widetilde{K}$.

The above regularity requirements could bring some useful consequences which are trivial to verify [4]:

1. M1

there exists a constant $\rho_{v}$ that only depends on $N$ and $\sigma$ such that $\rho_{v}h_{K}\leq h_{e}$ for every element $K$ and every edge $e$ of $K$. 2. M2

there exists a constant $\rho_{s}$ that only depends on $N$ and $\sigma$ such that for every element $K$ the following holds true

[TABLE]

where $\Delta(K)=\left\{K^{\prime}\in\mathcal{T}_{h}\ |\ K^{\prime}\cap K\neq\varnothing\right\}$ is the collection of the elements touching $K$. 3. M3

there exists a constant $\tau$ that only depends on $N$ and $\sigma$ such that for every element $K$, there is a disk (ball) inscribed in $K$ with center at the point $\boldsymbol{z}_{K}\in K$ and the radius $\tau h_{K}$. 4. M4

[Trace inequality] there exists a constant $C$ such that

[TABLE] 5. M5

[Inverse inequality] there exists a constant $C$ such that

[TABLE]

where $\mathbb{P}_{m}(\cdot)$ denotes the polynomial space of degree less than $m$.

Let us note that throughout the paper, $C$ and $C$ with a subscript are generic constants that may be different from line to line but are independent of the mesh size $h$ and how the interface cuts the mesh. Given a bounded domain $D\subset\mathbb{R}^{d}$ and an integer $r\geq 0$, we would use the standard notations and definitions for the spaces $H^{r}(D)$, $L^{r}(D)$ and their corresponding inner products and norms. Then we will use the following notations related to the partition:

[TABLE]

[TABLE]

Furthermore, we denote by $\mathcal{T}_{h}^{\Gamma}=\left\{K\in\mathcal{T}_{h}\ |\ K\cap\Gamma\neq\varnothing\right\}$ the set of the elements that are divided by $\Gamma$ and by $\mathcal{E}_{h}^{\Gamma}=\left\{e\in\mathcal{E}_{h}\ |\ e\cap\Gamma\neq\varnothing\right\}$ the set of the faces that are divided by $\Gamma$. We set $\mathcal{T}_{h}^{\backslash\Gamma}=\mathcal{T}_{h}\backslash\mathcal{T}_{h}^{\Gamma}$ and $\mathcal{E}_{h}^{\backslash\Gamma}=\mathcal{E}_{h}\backslash\mathcal{E}_{h}^{\Gamma}$. For an element $K\in\mathcal{T}_{h}^{\Gamma}$ we denote $\Gamma_{K}=K\cap\Gamma$.

We make the following assumptions about the mesh, which are actually easy to be fulfilled.

Assumption 1**.**

For any face $e\in\mathcal{E}_{h}^{\Gamma}$, the intersection $e\cap\Gamma$ is simply connected; that is, $\Gamma$ does not intersect an interior face multiple times.

Assumption 2**.**

For any element $K\in\mathcal{T}_{h}^{\Gamma}$, there exists a line (plane) $\widetilde{\Gamma}_{K}$ and a smooth function $\psi$ that maps $\widetilde{\Gamma}_{K}$ onto $\Gamma_{K}$.

Assumption 3**.**

For any element $K\in\mathcal{T}_{h}^{\Gamma}$, there exist two elements $K^{0}_{\circ},K^{1}_{\circ}\in\Delta(K)$ such that $K^{0}_{\circ}\subset\Omega^{0}$ and $K^{1}_{\circ}\subset\Omega^{1}$.

We note that Assumption 1 and Assumption 2 ensure the interface $\Gamma$ is well resolved by the mesh [36] and such similar geometric assumptions are commonly used in numerically solving interface problems [37, 20, 50, 46, 8]. In Fig. 2, we present some examples of cut elements to illustrate the assumptions.

For the given partition $\mathcal{T}_{h}$, we follow the idea in [30, 29] to define the reconstruction operator for solving the elliptic interface problem. First, for every element $K\in\mathcal{T}_{h}$, we specify its barycenter $x_{K}$ as a sampling point. Second, for each element $K\in\mathcal{T}_{h}^{i}\backslash\mathcal{T}_{h}^{\Gamma}(i=0,1)$, we will construct an element patch $S^{i}(K)$ for $K$. The element patch $S^{i}(K)$ is a set of elements and consists of $K$ and some elements around $K$. We start from assigning a threshold value $\#S(K)^{i}$ that is used to control the size of $S^{i}(K)$, and we setup the element patch $S^{i}(K)$ in a recursive manner. Let $S_{0}^{i}(K)=\left\{K\right\}$, then we define $S^{i}_{t}(K)$ as

[TABLE]

Once $S^{i}_{t}(K)$ has collected $\#S^{i}(K)$ elements, we stop the procedure and let $S^{i}(K)=S^{i}_{t}(K)$. Clearly, the cardinality of $S^{i}(K)$ is just the value $\#S^{i}(K)$. For any element $K\in\mathcal{T}_{h}^{\backslash\Gamma}$, we only construct one element patch which satisfies that if $K\in\mathcal{T}_{h}^{i}\backslash\mathcal{T}_{h}^{\Gamma}$ then $S^{i}(K)\subset\mathcal{T}_{h}^{i}$. For any element $K\in\mathcal{T}_{h}^{\Gamma}$, we assume that $K\in S^{0}(K_{\circ}^{0})$ and $K\in S^{1}(K_{\circ}^{1})$ where $K^{0}_{\circ}$ and $K^{1}_{\circ}$ are defined in Assumption 2. With $\#S^{0}(K_{\circ}^{0})$ and $\#S^{1}(K_{\circ}^{1})$ to be mildly greater than that in [30, 29], the assumption can be fulfilled according to the method of constructing the element patch. Consequently, for each element $K\in\mathcal{T}_{h}^{\Gamma}$, we have two element patches $S^{0}(K)=S^{0}(K_{\circ}^{0})$ and $S^{1}(K)=S^{1}(K_{\circ}^{1})$. In Appendix A, we present the detailed algorithm and give some examples to illustrate the construction of the element patch.

For any element $K\in\mathcal{T}_{h}$, we denote by $\mathcal{I}_{K}^{i}(i=0,1)$ the set of sampling points located inside $S^{i}(K)$,

[TABLE]

For any function $g\in C^{0}(\Omega)$ and an element $K\in\mathcal{T}_{h}$, we seek a polynomial $\mathcal{R}_{K}^{i}g$ defined on $S^{i}(K)$ of degree $m$ by solving the following least squares problem:

[TABLE]

The existence and uniqueness of the solution to (3) are decided by the position of the sampling nodes in $\mathcal{I}_{K}^{i}$. Here we follow [31] to make the following assumption:

Assumption 4**.**

For any element $K\in\mathcal{T}_{h}$ and $p\in\mathbb{P}_{m}(S^{i}(K))$,

[TABLE]

This assumption actually rules out the situation that all the points in $\mathcal{I}^{i}_{K}$ are located on an algebraic curve of degree $m$. Definitely, this assumption requires the cardinality $\#S^{i}(K)$ shall be greater than $\text{dim}(\mathbb{P}_{m})$. Hereafter, we always require this assumption holds.

Since the solution to (3) is linearly dependent on $g$, we define two interpolation operators $\mathcal{R}^{i}$ for $g$:

[TABLE]

Given $\mathcal{R}^{i}(i=0,1)$ and $g\in C^{0}(\Omega)$, the function $g$ is mapped to a piecewise polynomial function of degree $m$ on $\mathcal{T}_{h}^{i}$. We denote by $V_{h}^{i}$ the image of the operator $\mathcal{R}^{i}$. For any element $K$, We pick up a function $w_{K}^{i}(\boldsymbol{x})\in C^{0}(\Omega)$ such that

[TABLE]

It should be noted that in element $K$ we do not care about the values of $w_{K}^{i}(\boldsymbol{x})$ at $\boldsymbol{x}\neq\boldsymbol{x}_{K}$ and such continuous functions obviously exist. Then it is easy to check that $V_{h}^{i}=\mathrm{span}\left\{\lambda_{K}^{i}\ |\ \lambda^{i}_{K}=\mathcal{R}^{i}w_{K}^{i}\right\}$, and one can write the operator $\mathcal{R}^{i}$ in an explicit way:

[TABLE]

In Appendix B, we present a one-dimensional example to show more details of construction of $\lambda_{K}^{i}$ and its computer implementation.

Remark 1**.**

The computational cost of constructing the approximation spaces $V_{h}^{0}$ and $V_{h}^{1}$ mainly consists of two parts. The first is the construction of element patches. We adopt a recursive strategy on every element for construction as we illustrate in Appendix A. The number of recursive steps is related to the order $m$ and in numerical experiments we take $1\leq m\leq 3$. We also list the values $\#S^{i}(K)$ that are used in numerical experiments in Section 4. In this case at most $6$ recursive steps are required on each element. Hence this part is very cheap. The second part is to solve the function $\lambda_{K}^{i}$ on each element. In this part, the main step is to solve an inverse of a $\#S^{i}(K)\times\#S^{i}(K)$ matrix, as we demonstrate in Appendix B. Thus, the computational cost of the second part is still small.

The operators $\mathcal{R}^{i}(i=0,1)$ are defined for functions in $C^{0}(\Omega)$, while we only concern the case for the functions in $H^{t}(\Omega_{0}\cup\Omega_{1})(t\geq 2)$. Hence, we choose two extension operators to extend the functions in $H^{t}(\Omega_{0}\cup\Omega_{1})$ to be defined in $H^{t}(\Omega)$ [1]. For any function $w\in H^{t}(\Omega_{0}\cup\Omega_{1})$, there exist two operators $E^{i}:H^{t}(\Omega_{i})\rightarrow H^{t}(\Omega)$ such that $(E^{i}w)|_{\Omega_{i}}=w$ and

[TABLE]

Now let us study the approximation property of the operator $\mathcal{R}^{i}$. We define $\Lambda(m,S^{i}(K))$ for all element patches as

[TABLE]

We note that under some mild conditions on $S^{i}(K)$, $\Lambda(m,S^{i}(K))$ admits a uniform upper bound $\Lambda_{m}$, which is crucial in the convergence analysis. We refer to [30, Assumption A] for the geometrical conditions on element patches. These conditions in fact exclude the case that the points in $\mathcal{I}_{K}^{i}$ are very close to an algebraic curve of degree $m$. We also proved that if the size of the element patch $S^{i}(K)$ is greater than a certain number, then the geometrical conditions will be satisfied, see [30, Lemma 6] and [31, Lemma 3.4]. We note that that this number is usually too great and we prefer not to adopt it in the implementation. In numerical tests, we observe that out method can still work very well under the case that the value $\#S^{i}(K)$ is far less than the theoretical value. In Section 4, we list the values of $\#S^{i}(K)$ that are used in numerical examples. In addition, we refer to [31] for some numerical experiments about the size of the element patch and the upper bound $\Lambda_{m}$.

Remark 2**.**

For a special case when $\#S^{i}(K)=\text{dim}(\mathbb{P}_{m})$, we may replace the constant $\Lambda(m,S^{i}(K))$ by the Lebesgue constant [43, p.24]. In this case, the solution to the problem (3) is the Lagrange interpolation polynomial. Unfortunately, we have little knowledge of the Lebesgue constant in two or three dimensions.

With $\Lambda_{m}$, we have the local approximation error estimates.

Theorem 1**.**

Let $g\in H^{t}(\Omega_{0}\cup\Omega_{1})(t\geq 2)$, there exist constants $C$ such that for any $K\in\mathcal{T}_{h}^{i}(i=0,1)$ the following estimates hold true:

[TABLE]

where $s=\min(t+1,m)$.

Proof.

It is a direct consequence of [30, lemma 2.4] or [29, lemma 2.5]. ∎

Finally, we give the definition of our approximation space $V_{h}$ by concatenaing the two spaces $V_{h}^{0}$ and $V_{h}^{1}$. Let us define a global interpolation operator $\mathcal{R}$: for any function $w\in H^{t}(\Omega_{0}\cup\Omega_{1})$, $\mathcal{R}w$ is piecewise defined by

[TABLE]

The image of $\mathcal{R}$ is actually our new approximation space $V_{h}$. We notice that for any function $w\in H^{t}(\Omega_{0}\cup\Omega_{1})$, $\mathcal{R}w$ is a combination of $\mathcal{R}^{0}w$ and $\mathcal{R}^{1}w$ that $(\mathcal{R}w)|_{K^{i}}=(\mathcal{R}^{i}w)_{K^{i}}(i=0,1)$, and the approximation error estimates of $\mathcal{R}$ are the direct consequence from (5).

3. Approximation to Elliptic Interface Problem

We consider the standard elliptic interface problem: find $u$ in $H^{2}(\Omega_{0}\cup\Omega_{1})$ such that

[TABLE]

where $\beta$ is a positive constant function on $\Omega_{i}(i=0,1)$ but may be discontinuous across the interface $\Gamma$, and $\boldsymbol{\mathrm{n}}_{\Gamma}$ denotes the the unit normal of $\Gamma$ pointing to $\Omega_{0}$ (see Fig 1). The source term $f$, the Dirichlet data $g$ and the jump term $a$, $b$ are assumed to be in $L^{2}(\Omega)$, $H^{3/2}(\partial\Omega)$, $H^{3/2}(\Gamma)$, $H^{1/2}(\Gamma)$, respectively, to ensure (6) has a unique solution. We refer to [44, 27, 42, 26] for more details. In (6), the jump operator $[[\cdot]]$ takes the standard sense in DG framework. More precisely, we define the jump operator $[[\cdot]]$ and average operator $\{\cdot\}$ as below,

[TABLE]

where $v$ is a scalar-valued function and $\boldsymbol{q}$ is a vector-valued function. For $e\in\mathcal{E}_{h}^{\circ}$, we let $K_{+}$ and $K_{-}$ be two neighbouring elements that share a common face $e$. $\boldsymbol{\mathrm{n}}_{K_{+}}$ and $\boldsymbol{\mathrm{n}}_{K_{-}}$ are the unit outer normal on $e$ corresponding to $\partial K_{+}$ and $\partial K_{-}$, respectively. In the case $e\in\mathcal{E}_{h}^{b}$, we let $e$ be a face of the element $K$.

Now we define the bilinear form $b_{h}(\cdot,\cdot)$ and the linear form $l_{h}(\cdot)$:

[TABLE]

for $\forall u_{h},v_{h}\in W_{h}$, and

[TABLE]

for $\forall v_{h}\in W_{h}$, where $W_{h}$ denotes the following broken Sobolev space

[TABLE]

The penalty parameter $\eta$ is nonnegative and will be specified later on. For any $v_{h}\in W_{h}$, let us define a DG energy norm $|\!|\!|\cdot|\!|\!|$ as

[TABLE]

where

[TABLE]

The approximation problem to the elliptic interface problem (6) is then defined as: find $u_{h}\in V_{h}$ such that

[TABLE]

An immediate consequence from the definitions of the bilinear form $b_{h}(\cdot,\cdot)$ and the linear form $l_{h}(\cdot)$ is the validity of the Galerkin orthogonality, which plays a key role in the error estimate later on.

Lemma 1**.**

Let $u\in H^{2}(\Omega_{0}\cup\Omega_{1})$ be the exact solution and let $u_{h}\in V_{h}$ be the solution to (8), the Galerkin orthogonality holds true:

[TABLE]

Proof.

By $[[u]]=0$ on any $e\in\mathcal{E}_{h}^{\circ}$ and $[[u]]=g\boldsymbol{\mathrm{n}}$ on any $e\in\mathcal{E}_{h}^{b}$, we observe that

[TABLE]

Applying integration by parts, we have that

[TABLE]

Combining above two equations implies $b_{h}(u_{h},v_{h})=b_{h}(u,v_{h})$, which completes the proof. ∎

Next we verify the boundedness and coercivity of the bilinear form $b_{h}(\cdot,\cdot)$ with respect to the energy norm $|\!|\!|\cdot|\!|\!|$. For this purpose, we need to estimate the error on the interface. Here we first give the discrete trace inequality, which is crucial in the error estimate.

Lemma 2**.**

For any $K\in\mathcal{T}_{h}^{\Gamma}$, there exists a constant $C$ such that

[TABLE]

where $\partial K^{i}=(\partial K)^{i}\cup\Gamma_{K}$.

Proof.

Since $K\in\mathcal{T}_{h}^{\Gamma}$, we have that the patch $S^{i}(K)$ is the same as the patch $S^{i}(K_{\circ}^{i})$. From the definition of the least squares problem (3), it is clear that the solution to (3) on $S^{i}(K)$ is the same as the solution to (3) on $S^{i}(K_{\circ}^{i})$. Particularly, $\nabla^{\alpha}v_{h}|_{K^{i}}$ and $\nabla^{\alpha}v_{h}|_{K_{\circ}^{i}}$ are exactly the same polynomial which is denoted as $\tilde{p}$. Based on M3, there exists a constant $\hat{\tau}$ such that $B(\boldsymbol{z}_{K_{\circ}^{i}},\hat{\tau}h_{K_{\circ}^{i}})\subset K_{\circ}^{i}$, where $B(\boldsymbol{z},r)$ is a ball with center at $\boldsymbol{z}$ and radius $r$. From Assumption 2, we have that $K\in\Delta(K_{\circ}^{i})$. By the mesh regularity M2, there exists a constant $\tilde{\tau}$ such that $\partial K^{i}\subset B(\boldsymbol{z}_{K_{\circ}^{i}},\tilde{\tau}h_{K_{\circ}^{i}})$ and there exists a constant $C$ such that $h_{K}\leq Ch_{K_{\circ}^{i}}$. We note that here the constants $\hat{\tau}$, $\tilde{\tau}$ and $C$ only depend on $N$ and $\sigma$. We further deduce that

[TABLE]

The third inequality follows from the inverse inequality $\|\hat{p}\|_{L^{\infty}(B(0,1))}\leq C\|\hat{p}\|_{L^{2}(B(0,\hat{\tau}/\tilde{\tau}))}$ for any $\hat{p}\in\mathbb{P}_{m}(B(0,1))$ and the pullback using the bijective affine map from $B(\boldsymbol{z}_{K_{\circ}^{i}},\tilde{\tau}h_{K_{\circ}^{i}})$ to $B(0,1)$. As $\Gamma$ is of class $C^{2}$, it is easy to show (cf. [12, 50]) $|\Gamma_{K}|\leq Ch_{K}^{d-1}$. We complete the proof by observing $|\partial K^{i}|\leq h_{K}^{d-1}$ and $|B(\boldsymbol{z}_{K_{\circ}^{i}},\tilde{\tau}h_{K_{\circ}^{i}})|\leq Ch_{K_{\circ}^{i}}^{d}$. ∎

Lemma 3**.**

There exists a positive constant $h_{0}$ independent of $h$ and the location of the interface such that for all $h\leq h_{0}$ and any element $K\in\mathcal{T}_{h}^{\Gamma}$, the following trace inequality holds true:

[TABLE]

See the proof of this lemma in [50, 20, 51].

Now we are ready to claim the continuity and coercivity of the bilinear form $b_{h}(\cdot,\cdot)$.

Theorem 2**.**

Let $b_{h}(\cdot,\cdot)$ be the bilinear form defined in (7) with sufficiently large $\eta$. Then there exist positive constants $C$ such that

[TABLE]

Proof.

By Cauchy-Schwartz inequality, for $\forall u,v\in W_{h}$ we directly obtain that

[TABLE]

which directly gives us the continuity result (12).

To obtain (13), we first define a weaker norm $\|\cdot\|_{*}$ which is a more natural one for analyzing coercivity. For any $w_{h}\in V_{h}$, $\|\cdot\|_{*}$ is given by

[TABLE]

From the trace estimate (1) and the inverse inequality (2), we immediately obtain that

[TABLE]

By the trace estimate (10) and the mesh regularity M1, we have that

[TABLE]

and

[TABLE]

The above inequalities give us

[TABLE]

which actually indicates $|\!|\!|w_{h}|\!|\!|\leq C\|w_{h}\|_{*}$ and the equivalence of $|\!|\!|\cdot|\!|\!|$ and $\|\cdot\|_{*}$ restricted on $V_{h}$.

Then we consider to bound the trace terms in the bilinear form with respect to the norm $\|\cdot\|_{*}$. For the face $e\in\mathcal{E}_{h}^{\circ}$, we let $e$ be shared by two neighbouring elements $K^{-}$ and $K^{+}$. For any $e\in\mathcal{E}_{h}^{\circ}\cap\mathcal{E}_{h}^{\backslash\Gamma}$, we apply the Cauchy-Schwartz inequality to get that

[TABLE]

for any $\varepsilon>0$. For any $e\in\mathcal{E}_{h}^{\circ}\cap\mathcal{E}_{h}^{\Gamma}$ and $i=0,1$, we deduce that

[TABLE]

By trace inequality (10) and (1), for any $e\in\mathcal{E}_{h}^{\circ}$ we have

[TABLE]

Together with (14) - (16), we obtain that

[TABLE]

and

[TABLE]

For any $e\in\mathcal{E}_{h}^{b}$, it is similar to derive that

[TABLE]

Further for any $K\in\mathcal{T}_{h}^{\Gamma}$, we again apply the trace estimate (10) to obtain that

[TABLE]

The inequality (17) yields that

[TABLE]

Combining all above inequalities, we conclude that there exists a constant $C$ such that

[TABLE]

for any $\varepsilon>0$. We can directly let $\varepsilon=\beta/(2C)$ and select a sufficiently large $\eta$ to ensure $b_{h}(v_{h},v_{h})\geq C\|v_{h}\|_{*}^{2}$, which completes the proof. ∎

Remark 3**.**

To ensure the stability near the interface, some unfitted methods [51, 37, 20] may require a weighted average $\{v\}=\kappa_{0}v|_{\Omega_{0}}+\kappa_{1}v|_{\Omega_{1}}$ where $\kappa_{0}$ and $\kappa_{1}$ are the cut-dependent parameters like $\kappa_{i}=|K_{i}|/|K|(i=0,1)$ for elements in $\mathcal{T}_{h}^{\Gamma}$. In our method, another advantage is just taking the arithmetic one could also guarantee the stability and we note that this advantage is brought by the patch reconstruction. In addition, the analysis can be adapted to their choices without any difficulty.

Now let us give the approximation error in the DG energy norm $|\!|\!|\cdot|\!|\!|$.

Lemma 4**.**

Let $u\in H^{t}(\Omega_{0}\cup\Omega_{1})$ with $t\geq 2$, there exists a constant $C$ such that

[TABLE]

where $s=\min(m+1,t)$.

Proof.

From (5), it is trivial to obtain

[TABLE]

Then using trace inequality (1) and (5), for any $K\in\mathcal{T}_{h}^{i}(i=0,1)$ we have

[TABLE]

From the above two inequalities and (4), we could conclude

[TABLE]

Finally we use (11) to bound the error on the interface. For any $K\in\mathcal{T}_{h}^{\Gamma}$, we obtain

[TABLE]

A summation over all $K\in\mathcal{T}_{h}^{\Gamma}$ gives us

[TABLE]

Similarly, we could yield

[TABLE]

Combining all the inequalities above gives the error estimate (18), which completes the proof. ∎

We are now ready to prove a priori error estimates.

Theorem 3**.**

Let $u\in H^{t}(\Omega_{0}\cup\Omega_{1})$ with $t\geq 2$ be the exact solution to (6) and let $u_{h}\in V_{h}$ be the solution to (8), then there exist constants $C$ such that the following error estimates hold true:

[TABLE]

and

[TABLE]

where $s=\min(m+1,t)$.

Proof.

Together with the Galerkin orthogonality (9), boundedness (12) and coercivity (13) of the bilinear form $b_{h}(\cdot,\cdot)$ we could have a bound of $|\!|\!|u-u_{h}|\!|\!|$. For any $v_{h}\in V_{h}$, we obtain that

[TABLE]

Hence,

[TABLE]

Combining (18) immediately gives us the estimate (19).

Finally we obtain the optimal convergence order in $L^{2}$ norm with the standard duality argument. Let $\phi\in H^{2}(\Omega_{0}\cup\Omega_{1})$ be the solution of

[TABLE]

and satisfy [5]

[TABLE]

We denote by $\phi_{I}=\mathcal{R}\phi$ the interpolant of $\phi$. Then together with the Galerkin orthogonality (9) we deduce that

[TABLE]

The estimate (20) is obtained by elminating $\|u-u_{h}\|_{L^{2}(\Omega)}$, which completes the proof. ∎

4. Numerical Experiments

In this section, we present some numerical results by solving some benchmark elliptic interface problems. For each case, the source term $f$, the Dirichlet boundary data $g$ and the jump term $a$, $b$ are given according to the solutions. We construct the spaces of order $1\leq m\leq 3$ to solve each problem. For simplicity, we take the $\#S^{i}(K)$ uniformly for all elements and we list the values of $\#S^{i}(K)$ for all $m$ that are used in all experiments in Tab 1. A direct sparse solver is used to solve the resulting sparse linear system. The interface in all numerical experiments is described by a given level set function $\phi(\boldsymbol{x})$.

4.1. 2D Example

Example 1.

We first consider the classical interface problem on the square domain $(-1,1)\times(-1,1)$ with a circular interface $\phi(x,y)=x^{2}+y^{2}-r^{2}$ with radius $r=0.5$ (see Fig 4). The exact solution and coefficient are chosen to be

[TABLE]

With $b=10$, $u$ is continuous over $\Omega$.

By using a series of quasi-uniform triangular meshes, the $L^{2}$ norm and DG energy norm of the error in the approximation to the exact solution with mesh size $h=1/5,1/10,\ldots,1/80$ are reported in Fig 5. For each fixed $m$, we observe that the errors $\|u-u_{h}\|_{L^{2}(\Omega)}$ and $|\!|\!|u-u_{h}|\!|\!|$ converge to zero at the rate $O(h^{m+1})$ and $O(h^{m})$ as the mesh is refined, respectively. Such convergence rates are consistent with the theoretical results.

Example 2.

In this example, we consider the same interface and the same domain as in Example 1. The analytical solution $u(x,y)$ and the coefficient are defined in the same way as in Example 1. But we solve the elliptic interface problem based on a sequence of polygonal meshes as shown in Fig 6, which are generated by PolyMesher [45].

The numerically detected convergence orders are displayed in Fig 7 for both error measurements. It is clear that the orders of convergence in $L^{2}$ norm and DG energy norm are $O(h^{m+1})$ and $O(h^{m})$, respectively, which again are in agreement with the theoretical predicts.

For the Example 3 - 6, the computational domain is $(-1,1)\times(-1,1)$ and we solve the test problems on a sequence of triangular meshes with mesh size $h=1/5,1/10,\cdots,1/80$.

Example 3.

In this case, we consider the problem in [54] which contains the strongly discontinuous coefficient $\beta$ to test the robustness of the proposed method. We consider the elliptic problem with an ellipse interface (see Fig 8),

[TABLE]

The exact solution and the coefficient are given as

[TABLE]

There is a large jump in $\beta$ across the interface $\Gamma$, which may lead to an ill-conditioned linear system. We still use the direct sparse solver to solve the resulting sparse linear system and our method shows the robustness for this case. As can be seen from Fig 9, the computed rates of convergence match with the theoretical analysis.

Example 4.

In this example, we consider solving the elliptic problem with a kidney-shaped interface [25], which is governed by the following level set function

[TABLE]

The boundary data and source term are derived from the exact solution and coefficient

[TABLE]

We present numerical results in Fig 11 and the predicted convergence rates for both norms are verified.

Example 5.

Next, we consider a standard test case with an interface consisting of both concave and convex curve segments [54]. The interface is parametrized with the polar angle $\theta$

[TABLE]

The exact solution is selected to be

[TABLE]

The convergence of the numerical solutions is displayed in Fig 13. Again we observe optimal rates of convergence for both norms as the mesh size is decreased.

Example 6.

In this case, we investigate the performance of our proposed method when dealing with the problem with low regularities. The interface can be found in [22], which is governed by the following level set function

[TABLE]

We note that the interface is only Lipschitz continuous and it has a kink at $(0,0)$, see Fig 14. The analytical solution $u(x,y)$ is given by

[TABLE]

We choose $\beta=1$ over the domain $(-1,1)\times(-1,1)$. The solution $u(x,y)$ is $C^{2}$ continuous but not $C^{3}$ continuous across the line $x+y=1$. The numerical errors in terms of $L^{2}$ norm and DG energy norm are gathered in Tab 2. It is observed that when $m=1,2$ the numerical solutions converge optimally with rate $O(h^{m+1})$ for $L^{2}$ norm and $m$ order for DG energy norm, which matches with the fact that the exact solution $u$ belongs to $H^{3}(\Omega_{0}\cup\Omega_{1})$. When $m=3$ the computed orders of convergence in $\|\cdot\|_{L^{2}(\Omega)}$ and $|\!|\!|\cdot|\!|\!|$ are about $O(h^{3.5})$ and $O(h^{2.5})$, respectively. A possible explanation of the convergence orders can be traced to lack of $H^{4}$-regularity of the exact solution on the domain $\Omega_{1}$.

4.2. 3D Example

Example 7.

Here we consider a three-dimensional elliptic interface problem. The domain $\Omega$ is $(0,1)^{3}$ and the spherical interface is given by

[TABLE]

with radius $r=0.35$. We select $\beta=1$ in the whole domain and the exact solution is taken as

[TABLE]

We adopt a family of tetrahedral meshes with mesh size $h=1/4$, $1/8$, $1/16$, $1/32$ to solve the interface problem (see Fig 15). The numerical solutions on the meshes with $h=1/16$ and $h=1/32$ are depicted in Fig 16 and these two solutions are obtained with the accuracy $m=3$. We display the slices at $y=0.5$ and $z=0.5$ of the numerical approximations on both meshes and both solutions significantly involve a discontinuity across a spherical, which are accordant with the interface. The convergence rates under both norms are shown in Fig 17. Clearly, the numerical results are still consistent with our theoretical predictions.

Example 8.

In this example, we consider a three-dimensional elliptic problem [49] with a smooth interface that is governed by the following level set function (see Fig 18),

[TABLE]

The exact solution is choose to be

[TABLE]

The domain $\Omega$ is taken to be $(0,1)^{3}$ and the coefficient $\beta$ is fixed as $1$. We also use the tetrahedral meshes with mesh size $h=1/4$, $1/8$, $1/16$, $1/32$ for solving the problem (see Fig 15). The slices of the numerical solution with the accuracy $m=3$ on the tetrahedral mesh with $h=1/32$ at $x=0.5$ and at $z=0.5$ are depicted in Fig 19. It is clear that the discontinuity of the numerical solution sketches a curve which matches with the interface given by the level set function (see 18). We also display the convergence history of the numerical approximation under both $L^{2}$ norm and DG energy norm in Fig 20. The convergence rate of $L^{2}$ error may seem less than the predicted value when the accuracy $m=1$. The rate is gradually more close to the theoretical value and we may expect the rate would go back to $O(h^{2})$ as the mesh size tends to zero. For $m=2$ and $m=3$, the computed convergence rates under both error measurements are in agreement with the theoretical results.

4.3. Integrals on Cut Element

In our method, computing the following types of integrals defined on the cut element is an important issue,

[TABLE]

where $K\in\mathcal{T}_{h}^{\Gamma}$ is a cut element and $K^{0}=K\cap\Omega_{0}$, $K^{1}=K\cap\Omega_{1}$ and $\Gamma_{K}=K\cap\Gamma$. Here we list two numerical methods for computing these integrals. The first is we generate highly accurate quadrature points and weights corresponding to the domain $K^{0}$, $K^{1}$ and the interface $\Gamma_{K}$. We refer to [15, 39, 23] for some approaches about finding such quadrature points and weights. The computational cost of the first method is much more expensive than ordinary numerical quadrature methods. The second one is we approximate the interface $\Gamma_{K}$ by planes or lines inside the element $K$, see Fig  21 for an example. In this case, we only need to generate quadrature points and weights for polygons or polyhedrons. The computational cost is much less than the first method but the result is less accurate. We refer to [47] for more details about this method.

Here we make a comparison between two methods. We solve the Example 7 by both two numerical quadrature methods. We call the C subroutines in PHG package [52, 15] to generate highly accurate quadrature points and weights for the cut tetrahedrons. For the second methods and for element $K\in\mathcal{T}_{h}^{\Gamma}$, we let $K^{0}$ be approximated by $\widetilde{K}^{0}$ and let $K^{1}$ be approximated by $\widetilde{K}^{1}$. The actual computational domains $\widetilde{\Omega}_{0}$ and $\widetilde{\Omega}_{1}$ are then given as

[TABLE]

We list the $L^{2}$ errors $\|u-u_{h}^{1}\|_{L^{2}(\Omega_{0}\cup\Omega_{1})}$ and $\|u-u_{h}^{2}\|_{L^{2}(\widetilde{\Omega}_{0}\cup\widetilde{\Omega}_{1})}$ in Tab 3, where $u_{h}^{1}$ and $u_{h}^{2}$ are the numerical solutions obtained by the first and second numerical quadrature methods, respectively. From Tab 3, we observe that the two errors are gradually closer to each other when the mesh size tends to zero. We note that both quadrature methods work in our numerical scheme and the first one is more accurate but much more computational cost is required.

4.4. Efficiency comparison

Hughes et al. [24] point out that the number of unknowns of a discretized problem is a proper indicator for the efficiency of a numerical method. To show the efficiency in DOFs of our method, we make a comparison among the unfitted DG method [37], the unfitted penalty finite element method [50, 51] and our method by solving the two-dimensional elliptic interface problem. The first method adopts the standard discontinuous finite element space, and the second method employs the traditional continuous finite element space. The solution and the partition are taken from Example 1. In Fig 22, we plot the $L^{2}$ norm of the error of three methods against the number of degrees of freedom with $1\leq m\leq 3$.

One see that for the low orders of approximation($m=1$), the penalty FE method is the most efficient method. For $m=2$, our method shows almost the same efficiency as the penalty FE method. For the high order accuracy($m=3$), our method performs better than the other methods.

5. Conclusion

We proposed a new discontinuous Galerkin method for elliptic interface problem. The approximation space is constructed by solving the local least squares problem. We proved optimal convergence orders in both $L^{2}$ norm and DG energy norm. A series of numerical results confirm our theoretical results and exhibit the flexibility, robustness and efficiency of the proposed method.

Acknowledgements

The authors would like to thank the anonymous referees sincerely for their constructive comments that improve the quality of this paper. This research was supported by the Science Challenge Project (No. TZ2016002) and the National Science Foundation in China (No. 11971041).

Appendix A Construction of Element Patch

Here we present the algorithm to the construction of the element patch in Alg  1 and we also give some an example of constructing element patches. We consider a circular interface. Let $\Omega_{1}$ be the domain inside the circle and $\Omega_{0}=\Omega\backslash\Omega_{1}$. For element $K\in\mathcal{T}_{h}^{0}\backslash\mathcal{T}_{h}^{\Gamma}$, the construction of $S^{0}(K)$ is presented in Fig 24.

For element $K\in\mathcal{T}_{h}^{1}\backslash\mathcal{T}_{h}^{\Gamma}$, the construction of $S^{1}(K)$ is presented in Fig 24.

Appendix B 1D Example

Here we present a one-dimensional example to illustrate our method. We consider the interval $\Omega=[-1,1]$ which is divided into two parts $\Omega_{0}=(-1,-0.2)$ and $\Omega_{1}=(-0.2,1)$. We partition $\Omega$ into 8 elements $\left\{K_{1},K_{2},\cdots,K_{8}\right\}$ with uniform spacing.

$\left\{x_{1},x_{2},\cdots,x_{8}\right\}$ are the set of collocations where $x_{i}$ is the midpoint of the element $K_{i}$. Since $\mathcal{T}_{h}^{\Gamma}=\left\{K_{4}\right\}$, we construct element patches for elements in $\mathcal{T}_{h}^{\backslash\Gamma}$. The element patches could be constructed as

[TABLE]

[TABLE]

Then for element $K_{4}$ it is clear that $K_{\circ}^{0}=K_{3}$ and $K_{\circ}^{1}=K_{5}$, and the element patches of $K_{4}$ are

[TABLE]

Then we would solve the least squares problem on every patch. We take $S^{0}(K_{3})$ for an example, for a continuous function $g$ and $m=1$ the least squares problem is written as

[TABLE]

It is easy to get the unique solution

[TABLE]

where

[TABLE]

We note that the matrix $(A^{T}A)^{-1}A^{T}$ has no relationship to the function $g$ and contains all information of all $\lambda_{K}^{i}$ on element $K_{2}$. Hence we store the matrix $(A^{T}A)^{-1}A^{T}$ for every element patch to represent all $\lambda_{K}^{i}$. It is in the same way when we deal with the high dimensional problem.