An Unfitted Discontinuous Galerkin Method for Elliptic Interface Problems

An unfitted discontinuous Galerkin method is proposed for the elliptic interface problems. Based on a variant of the local discontinuous Galerkin method, we obtain the optimal convergence for the exact solution u in the energy norm and its flux p in the L norm. These results are the same as those in the case of elliptic problems without interface. Finally, some numerical experiments are presented to verify our theoretical results.


Introduction
Elliptic interface problems are often encountered in many multiphysics and multiphase applications in science computing and engineering.For example, second order elliptic equations with discontinuous coefficients are often used to model problems in material sciences and fluid dynamics when two or more distinct materials or fluids with different conductivities, densities, or permeability are involved.It is well known that, when the interface is smooth enough, the solution of elliptic interface problems has higher regularity in individual material or fluid region than in the entire physical domain.
To numerically solve such interface problems, first we need to generate a mesh.One approach is to use a body fitted mesh.However, for those problems where the interface moves with time, repeated remeshing of the domain to obtain a fitted mesh is very costly.Another one is to use an unfitted grid independent of the location of the interface.This technique is particularly preferred to simulate time-dependent problems with moving interfaces.The major advantage for using an unfitted mesh is that it avoids repeatedly remeshing the domain for fitting the moving interfaces.
As for fitted mesh method for elliptic problems with interface, Chen and Zou in [1] considered the finite element method for solving elliptic and parabolic interface problems, and almost-optimal error estimates in the  2 norm and energy norm were obtained.In [2], the authors studied a class of discontinuous Galerkin method for elliptic interface problems, which was shown to be optimally convergent in  2 norm.Recently, a high-order HDG method was presented to solve elliptic interface problems by Huynh et al. in [3], which was extended to solve Stokes interface flow in [4].
Various unfitted grid methods for interface problems have been proposed in the literature.Finite difference methods are very popular unfitted grid methods due to their simplicity, for example, the immersed interface method [5,6], the immersed boundary method [7], the boundary condition capturing method [8], and many others.And there exist many works for finite element methods on unfitted grid as well.Babuška in [9] studied the elliptic interface problem on unfitted mesh and derived suboptimal convergence behavior.Li et al. proposed an immersed interface finite element method in [10], which modified the basis functions near interface to satisfy the homogeneous jump conditions.Later, this method was applied to elliptic and elasticity interface problems with nonhomogeneous jump conditions in [11,12].Recently, an unfitted finite element method based on Nitsche's method was presented by A. Hansbo and P. Hansbo in [13] and optimal order of convergence was proved without restrictions on the location of the interface relative to the mesh, which was used to solve incompressible elasticity with discontinuous modulus in [14].More recently, Massjung in [15] considered an ℎ unfitted discontinuous Galerkin method which was viewed as a generalization of Hansbo's method in [13].An optimal convergence rate with respect to ℎ and a suboptimal convergence rate with respect to  in energy norm were proved.Later, Wu and Xiao also presented an unfitted ℎ interface penalty finite element method, which was extended to the three dimensional case in [16].
The local discontinuous Galerkin (LDG) method was proposed by Cockburn and Shu in [17] to solve general timedependent convection-diffusion problems.Later, the method was carried to elliptic problems for mixed discontinuous Galerkin formulation by Castillo et al. in [18].The purpose of this paper is to extend the LDG method to a class of elliptic problems with a smooth interface.However, employing an unfitted mesh method, the interface can divide regular grid cells into degenerated subcells.If this situation happens, the standard inverse estimates can no longer be valid.In this paper, we use the weighted average instead of the arithmetic average in the classic LDG method to retrieve the inverse estimates (see Lemmas 8 and 10).Thus, we propose an unfitted discontinuous Galerkin method, based on a variant of LDG method.We prove the optimal convergence rate of the method for the exact solution  in the energy norm and its flux p in the  2 norm, respectively.
The rest of this paper is organized as follows.In Section 2 we propose our DG method and present some necessary preliminaries.We prove optimal order error estimates for our DG method in Section 3. In Section 4, some numerical experiments are presented to justify our theoretical results.Finally, conclusions are given in Section 5.
Let us now end this section with some notation to be used in this paper.We will use the standard notations for Sobolev spaces and norms in this paper (see [19,20]).In particular, for a bounded open set Ω in R 2 , if Ω = ⋃  =1 Ω  and Ω  ∩ Ω  = 0 ( ̸ = ), we denote by   (⋃  =1 Ω  ) the Sobolev space of functions  such that | Ω  ∈   (Ω  ), where   (Ω  ) denotes the standard Sobolev space with norm ‖ ⋅ ‖ ,Ω  .As usual we define the broken norm: Throughout the paper, the generic constant  is always independent of the mesh parameter ℎ.

Discontinuous Galerkin Method and Preliminaries
Let Ω be a bounded domain in R 2 with convex polygonal boundary Ω and Ω + ⊂ Ω an open domain with  2 boundary Γ = Ω + ⊂ Ω.Let Ω − = Ω \ Ω + (see Figure 1).We consider the following elliptic interface problem: where n is the outward pointing unit normal to Ω + and [] :=  + | Γ −  − | Γ is the jump of  across the interface Γ, where  ± is the restrictions of  on Ω ± .For the sake of simplicity, we assume that the coefficient  is a positive and piecewise constant; that is, | Ω ± =  ± > 0.
Regarding the regularity for the solution of the interface problem (1), we state without proof the following theorem.
Theorem 1 (cf.[1]).Assume that  ∈  2 (Ω) and  ∈  1/2 (Γ).Then problem (1) has a unique solution  ∈  2 (Ω + ∪ Ω − ), and the following a priori estimate holds: By introducing the flux p = ∇, the interface problem (1) can be rewritten into a first order system as Let T ℎ be a shape regular and locally quasi-uniform simplicial triangulation of Ω, generated independently of the location of the interface Γ.For the definition of shape regular and locally quasi-uniform, we refer to [20,21].Suppose T ℎ to be made of straight triangles  with diameter ℎ  .As usual, let ℎ := max ∈T ℎ ℎ  .The set of edges of the triangulation T ℎ is denoted by E ℎ , and E  ℎ is the set of interior edges of T ℎ .For any element  ∈ T ℎ , denote the part of  in Ω ± by  ± ; that is,  ± =  ∩ Ω ± .For any edge  ∈ E ℎ , let  ± =  ∩ Ω ± .We call the elements whose interiors are cut through by Γ "interface elements", and denote the set of the interface elements by T  ℎ .For an interface element  ∈ T  ℎ , assume that Γ  is the part of interface Γ intersecting .For the geometrical features of the interface Γ, we give the following plausible assumptions (cf.[13,15]).
Assumption 2. We assume that Γ intersects the boundary  of an element  ∈ T  ℎ exactly twice and each (open) edge at most once.Assumption 3. Let Γ ,ℎ be the straight line segment connecting the points of intersection between Γ and .We assume that Γ  is a function of length on Γ ,ℎ , in local coordinates: To formulate our numerical scheme, first we define two usual discontinuous finite element spaces as where   () denotes the space of polynomials of degree less than or equal to  ( ≥ 1) on each element .
We define our discontinuous finite element spaces as Following the notation of [22], let  be an interior edge shared by two triangles  1 and  2 in T ℎ .For a scalar valued function V, piecewise smooth on T ℎ with V  = V|   , we define the jump and the weighted average of V as Similarly, for a vector valued function w, piecewise smooth on T ℎ with w  = w|   , we set where n  is the unit normal of  pointing towards the outside of   and If  is an edge on the boundary of Ω, we define on  ⟦V⟧ = Vn, {{w}}  = w, (10) where n denotes the unit outer normal of  pointing towards the outside of Ω.
For the weight average across interface Γ of any piecewise smooth function V discontinuous on Γ, we set where  + +  − = 1 whose specific definitions will be given in Lemma 10.
We define the bilinear and linear forms And integral by parts yields that Hence, our DG approximation can be written as the following mixed variational problem: find (q ℎ ,  ℎ ) ∈ W ℎ × ℎ , such that For the exact  of the interface problem (1) and p = ∇, using Theorem 1, we have ⟦⟧ = 0, ⟦p⟧ = 0 on  ∈ E  ℎ , and [] = 0, [p ⋅ n] =  across Γ.Then the following consistency property holds: Let the mesh-dependent norm |‖ ⋅ ‖| ℎ be defined by Theorem 5. Suppose that the stabilization parameter  is positive; then the DG method (16) defines a unique approximate solution (p ℎ ,  ℎ ) ∈ W ℎ ×  ℎ .
The following lemma comes from the famous Stein's extension theorem.
Next, we state a standard approximation lemma.
Lemma 10.Let  ∈ T  ℎ ; then there exists a positive constant  such that where Proof.Under Assumptions 2 and 3, the trace inequality (29) follows from Lemma 3 in [13] and a scaling argument.

Error Estimate of Our DG Method
Now, we define interpolation operator Π by Π|  := Π  and To obtain the convergence result, we need to show the following approximate error estimates.
Lemma 12. Suppose that  ∈  +1 (Ω + ∪ Ω − ) and p ∈ ( +1 (Ω + ∪ Ω − )) 2 ; then the following approximate estimates hold: Proof.We only need to show inequality (30); inequality (31) can be shown similarly.For the first term on the left-hand side of (30), by the inequality ( 22) from Lemma 7 we obtain Summing over all triangles, it follows by Lemma 6 that Next, we estimate the second term on the left-hand side of (30) as follows.Suppose that  =  1 ∩  2 with  1 ,  2 ∈ T ℎ ; using the inequality (23) from Lemma 7 yields Summing over all edges, using Lemma 6 implies that Similarly, due to the inequality (24) from Lemma 7 and Lemmas 6 and 11, we find Thus the inequality (30) follows combining ( 33)-( 36), which completes the proof.
Next we present a priori error estimate of the exact solution  in the energy norm |‖ ⋅ ‖| ℎ and its flux p in the  2norm.
Theorem 13.Let (p, ) be the solution of (3) and (p ℎ ,  ℎ ) the solution of (16), respectively.Then, for (p, ) ∈ (  (Ω + ∪ Ω − )) Proof.By using the consistency property (17), we have  ( p , q ℎ ) −  (  , q ℎ ) =  ( p , q ℎ ) −  (  , q ℎ ) , Set Using (31) in Lemma 12 and -Cauchy-Schwartz inequality, we estimate the first term on the right-hand side of (39) as For the second term on the right-hand side of (39), by using (30) in Lemma 12 and trace inequalities from Lemmas 8 and 10, we obtain Similarly, by using the approximate estimation (31), the third term on (39) can be estimated as Then, by virtue of inequality (30) in Lemma 12, we bound the fourth term on (39) as which completes the proof.
Using the standard duality argument, we can obtain the following error estimate in the  2 norm.(50) Proof.Consider the following so-called adjoint problem Using Theorem 1, we get By introducing an auxiliary variable q = ∇, we obtain Since [] = 0 on Γ and ⟦⟧ = 0 on E ℎ , using integration by parts and the consistency property (17), we deduce Due to q = ∇, as in the proof of Theorem

Conclusions
In this paper, we have discussed an unfitted discontinuous Galerkin method for elliptic problems with a smooth interface.Based on a variant of local discontinuous Galerkin method, we have obtained the optimal order error estimates  in the energy norm in  and its flux p.And by using the standard duality argument the optimal convergence rate in  2 norm for  has also been derived.These presented results are the same as that of elliptic problems without interface.Finally, numerical experiments are given to confirm our theoretical results.We note that the convergence behavior in most existing works concerning the elliptic interface problems depends on the jump in the discontinuous coefficients.
It will be one of our future subjects to design an efficient numerical scheme that is robust with respect to the jump in the discontinuous coefficients.

Figure 3 :
Figure 3: The convergence rates of  2 norm in the exact  and its flux p for  1 :  2 = 1 : 10.

Figure 4 :
Figure 4: The convergence rates of  2 norm in the exact  and its flux p for  1 :  2 = 1 : 1000.

Figure 5 :Figure 6 :
Figure 5: The convergence rates of  2 norm in the exact  and its flux p for  1 :  2 = 10 : 1.