High order symmetric direct discontinuous Galerkin method for elliptic interface problems with fitted mesh

Journal Pre-proof High order symmetric direct discontinuous Galerkin method for elliptic interface problems with fitted mesh

Hongying Huang, Jin Li, Jue Yan

PII:

S0021-9991(20)30075-9

DOI:

https://doi.org/10.1016/j.jcp.2020.109301

Reference:

YJCPH 109301

To appear in:

Journal of Computational Physics

Received date:

20 May 2019

Revised date:

19 December 2019

Accepted date:

29 January 2020

Please cite this article as: H. Huang et al., High order symmetric direct discontinuous Galerkin method for elliptic interface problems with fitted mesh, J. Comput. Phys. (2020), 109301, doi: https://doi.org/10.1016/j.jcp.2020.109301.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier.

High order symmetric direct discontinuous Galerkin method for elliptic interface problems with fitted mesh Hongying Huang∗, Jin Li†, Jue Yan‡ February 10, 2020

Abstract In this article, we aim to develop a high order symmetric direct discontinuous Galerkin method solving elliptic interface problems with zero or non-zero solution jump and flux jump interface conditions. We focus on the case, in which the mesh is partitioned along the curved interface. The two interface jump conditions are naturally and simultaneously built into the numerical flux definition on the curved triangular elements’ edges that overlap with the interface. We obtain a stable and high order method, regardless of the combination of the two interface jump conditions. The two interface jump conditions are both weakly enforced in the scheme formulation. Optimal (k +1)th order L2 norm error estimate is proved for polygonal interfaces. A sequence of numerical examples are carried out to verify the optimal convergence of the symmetric direct discontinuous Galerkin method with high order P2 , P3 and P4 approximations. Uniform optimal convergence order, which is independent of the diffusion coefficient ratio inside and outside of the interface, is obtained. The symmetric direct discontinuous Galerkin method is shown to be capable of handling interface problems with complicated geometries.

1

Introduction

In this article, we aim to develop a symmetric direct discontinuous Galerkin (DDG) method [16] to solve the following elliptic interface problems, −∇ · (α∇u± (x)) = f ± (x), x ∈ Ω± , u = u0 , on ∂Ω, JuK = φ, on Γ, Jαun K = ψ, on Γ.

(1.1) (1.2) (1.3) (1.4)

For simplicity of presentation, we focus on two-dimensional problems. We have Ω ⊂ R2 denoting an open bounded domain. Let Γ be an interface that divides Ω into disjoint open subdomains, Ω− and Ω+ . Hence we have Ω = Ω+ ∪Ω− (see Figure 1). Assume that the domain boundary ∂Ω and the interface Γ are both Lipschitz continuous. The diffusion coefficient α(x) is assumed to be uniformly positive and smooth in each subdomain, but may be discontinuous across the interface, and there exist positive constants αl and αu satisfying, 0 < αl ≤ α(x) ≤ αu ,

x ∈ Ω.

(1.5)

Without the loss of generality, we consider two subdomains with diffusion coefficients piecewise defined as,   α− (x), x ∈ Ω− , α= (1.6)  α+ (x), x ∈ Ω+ . ∗ School of Mathematics, Physics and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang and Key Laboratory of Oceanographic Big Data Mining & Application of Zhejiang Province, Zhoushan, Zhejiang, China. e-mail: [email protected] † College of Science, North China University of Science and Technology, 21 Bohai Road, Caofeidian Xincheng, Tangshan 063210, Hebei, China. e-mail: [email protected] ‡ Department of Mathematics, Iowa State University, Ames, IA 50011. e-mail: [email protected] The research of this author is supported by NSF grant DMS-1620335.

1

Figure 1: Sample domain Ω = Ω+ ∪ Ω− ∪ Γ. On the interface Γ, the solution jump and flux jump conditions are given as JuK = u− − u+ = φ and Jαun K = α− ∇u− · n− + α+ ∇u+ · n+ . Here n+ and n− denote the outward-pointing unit normal vectors of Ω+ and Ω− , respectively (see Figure 1). The two non-zero jump conditions of (1.3) and (1.4) impose a challenging task to obtain a high order method solving the interface problem. Based on mesh partition over the underlined domain Ω as fitted or unfitted mesh, numerical methods are classified into two categories in literature. With mesh cut by the interface, we have the Immersed Boundary Method [27] [28], Immersed Interface Methods [20] [19], Immersed Finite Element methods or Immersed-Interface Finite element methods [21] [12] [22], and many others such as [9] [4] [14] [23] [1] [13] [35] [18], etc. In this article, we consider the case of body-fitted mesh (Figure 2) and high order Pk (k ≥ 2) polynomial approximations. Up to the authors’ knowledge, there is little result for numerical methods of higher than 2nd order to solve interface problem (1.1)-(1.4). Since the numerical solution of discontinuous Galerkin (DG) method is designed to be discontinuous across element edges, it is natural to consider DG method to incorporate the non-zero solution jump interface condition (1.3). Direct discontinuous Galerkin method [24] involves a concept of numerical flux u cn , which approximates the flux or the normal derivative un = ∇u · n on element edges. All these features imply that DDG method is a good choice to solve elliptic interface problems coupled with zero or non-zero solution jump (1.3) and flux jump (1.4) interface conditions. The focus of this article is to study the definition of numerical flux u cn to incorporate the two zero or non-zero jump conditions and obtain a uniform high order method. Different to [3] [31], purpose of the article is not toward mesh generation. Mesh is generated through MATLAB, even for challenging numerical tests involving complicated geometries. We should mention the recent work of Weak DG method [26], which is a second order method solving elliptic interface problems. In [24], we have developed the DDG method to discretize the diffusion term for time dependent convection diffusion equations. The key idea of DDG method is the concept of numerical flux u cn that approximates solution’s normal derivative un = ∇u · n over the discontinuous element edges. The numerical flux formula of u cn involves the solution jump JuK, the normal derivative average {{un }}, and second order normal derivative jump Junn K terms of the two piecewise polynomials. Thus, we can build the given two jump conditions of (1.3) and (1.4) into the definition of u cn , which is defined on element edges that overlap with the interface, and thereby we obtain a uniform high order method. The DDG method is closely related to the classical IPDG method [2] [32]. With low order piecewise constant and linear approximations, DDG method degenerates to the IPDG method. On the other hand, DDG method shows quite a few advantages over IPDG method with high order approximations (Pk with k ≥ 2). In [25], numerically we observe that a small fixed-penalty coefficient, which is independent of the approximation polynomial degrees, can be applied to stabilize the scheme for any high order approximations. Under the topic of maximum principle, DDG polynomial solutions are proved to satisfy the strict maximum principle on unstructured triangular meshes [7] with at least third order of accuracy. In [33], DDG methods are proved to be super convergent to the solution’s spatial derivative. On the other hand, we have not observed

2

Figure 2: A sample body-fitted triangular mesh partition of the domain Ω = Ω+ ∪ Ω− ∪ Γ. this super convergence phenomena with IPDG method. The study of DDG method has been recently extended to 2nd order elliptic problems in [16]. When comparing with IPDG method, which is the leading DG type diffusion solver, the symmetric DDG method [16] saves at least one digit on the mass matrix condition number and roughly saves 7-10% on CPU time for high order approximations. In [16], we carry out a numerical test (Example 4.3) to an elliptic interface problem with a straight-line interface and zero solution jump (φ = 0 in (1.3)) and zero flux jump (ψ = 0 in (1.4)) conditions. For this example, we do not need to modify the symmetric DDG scheme formulation, because the zero flux jump condition (ψ = 0) implies the continuity of the flux. The successful numerical experiment (Example 4.3 in [16]) does motivate us to further study symmetric DDG method for elliptic interface problems under general interface conditions of (1.3) and (1.4). Due to the discontinuous feature of DG method and the concept of αu dn , the numerical flux on element edges that overlap with the interface can be carefully designed. We consider body-fitted mesh aligned with the curved interface. Thus curved triangular elements with one curved edge are applied in the symmetric DDG scheme implementation. We pay special attention to the curved triangular elements, along the curved edge that the two zero or non-zero jump conditions of (1.3) and (1.4) are simultaneously built into the definition of the numerical flux. Thus, we can obtain a high order method solving the elliptic interface problem (1.1)-(1.4). The two interface conditions are both enforced in the weak sense. Different from classical finite element, in which the non-zero solution jump condition is manually removed by pasting the inside and outside solutions together or a homogenization step is needed, we do not differentiate the solution jump as being zero or nonzero. A uniform high order method is obtained, which is capable of handling any combination of zero or non-zero solution jump (1.3) and flux jump (1.4) interface conditions. Around the interface, we have curved triangular elements applied in the symmetric DDG method. An accurate enough volume integration approximation over the curved triangular element is applied in the symmetric DDG scheme formulation. A Gaussian-type formula is applied to approximate the line integration on the curved edge. Accuracy loss is observed if we simplify the volume integration approximation over straight triangles instead of the curved triangles. In this article, we focus on high order method development and skip the discussion on efficiency issues. A combination of symmetric DDG method around the interface and continuous finite element method away from the interface is a better option to save the total number of degrees of freedom and is more practical toward real applications. A simple symmetric DDG method is obtained, and the scheme is consistent with the original symmetric DDG method for smooth elliptic problems with no interface. Optimal (k + 1)th order L2 norm error estimate is obtained for polygonal interfaces. A sequence of numerical examples are carried out to verify the optimal (k + 1)th order convergence with high order P2 , P3 and P4 approximations. Numerical tests also demonstrate that uniform optimal convergence order, which is independent of the diffusion coefficient ratio inside and outside of the interface, is obtained. The symmetric DDG method is shown to be capable of handling interface problems with complicated geometries. Both mesh generation and symmetric DDG method implementation are carried out with MATLAB code. No special technique or handling are applied to the mesh around

3

the curved interface, and the symmetric DDG method can be easily extended to three-dimensional interface problems. For simplicity of presentation, we focus on two-dimensional problems in this paper. The rest of the article is organized as follows. In §2, we lay out the symmetric DDG scheme formulation to problem (1.1) - (1.4). Definition of the numerical flux on the curved triangular element edge that falls on the interface, which incorporates the two jump conditions, is presented. In §3 and §4, we present the stability and a priori error estimate with a standard energy norm and L2 norm. Finally numerical examples are shown in §5. Throughout this paper, we let | · |H s (E) and k · kH s (E) denote the semi-norm and the norm of H s (E), s ≥ 0. We have the solution u = (u+ , u− ) piecewise defined on Ω, such that u|Ω+ = u+ and u|Ω− = u− , and u+ ∈ H s (Ω+ ) and u− ∈ H s (Ω− ) for s ≥ 0.

2

Scheme formulation

We start with mesh partition, notations and the introduction of direct discontinuous Galerkin (DDG) finite element method. In this article, we consider the mesh partitioned along with the interface Γ. See Figure 2 for a sample body-fitted mesh partition. Let Th be a body-fitted and shape regular triangular mesh with Ω = Ω− ∪ Ω+ = ∪E∈Th E. Notice we have each triangular element E ∈ Th sitting either inside Ω+ or Ω− , and each triangle has at most two vertices on the interface Γ. For each element E ∈ Th , we have hE denote its diameter and h = maxE∈Th hE . When Γ is a curve and the triangle E has two vertices on Γ, E is a curved triangle with one curved edge. The triangulation Th can be decomposed into three categories: Th+ := {E ∈ Th |E ⊂ Ω+ , E has at most one vertex on Γ}, Th− := {E ∈ Th |E ⊂ Ω− , E has at most one vertex on Γ}, ThΓ := {E ∈ Th |E has two vertices on Γ}.

(2.1)

The edges of Th are defined as follows. A boundary edge of Th is the edge of ∂E ∩ ∂Ω, where E ∈ Th has two vertices on ∂Ω. An interface edge of Th is the edge of ∂E ∩ Γ, where E ∈ ThΓ has two vertices on Γ. The remaining edges of Th are called the interior edges. We denote Eh as the collection of all edges of triangular element E ∈ Th . We define EhI as the set of interior edges in Ω+ and Ω− , EhD as the set of boundary edges, and EhΓ as the set of interface edges. Thus we have Eh = EhI ∪ EhD ∪ EhΓ . For each edge e, we use he to represent its length. Regarding the regularity of the mesh, there exists a positive constant Ch such that for all E ∈ Th , and all e ∈ Eh , we have hE ≤ Ch he . We have Pk (E) representing the space of polynomials of degree at most k on element E. The DG solution space and test functions space is defined as, Vkh := {v ∈ L2 (Ω) : v|E ∈ Pk (E), ∀E ∈ Th }. To define the numerical flux of a DG method, we need to further introduce following notations. Suppose E and E 0 are two adjacent elements and share one common edge e. There are two traces of v along the edge e, where we add or subtract those values to obtain the average and the jump. We denote by n = (n1 , n2 )T the outward unit normal vector pointing from E into its neighbor element E 0 . Now the average and the jump of v on edge e, evaluated from current element E, are defined and denoted as follows, {{v}} =

1 (v|E + v|E 0 ) , 2

JvK = v|E 0 − v|E ,

∀e = ∂E ∩ ∂E 0 ,

1 (v|E + v|E 0 ) · n, Jv · nK = (v|E 0 − v|E ) · n, ∀e = ∂E ∩ ∂E 0 . 2 With no ambiguity, for the rest of this article we use the same letter u instead of the standard notation of uh to represent the piecewise DG polynomial solution. {{v · n}} =

2.1

Symmetric direct discontinuous Galerkin method

We use the following Poisson equation over one simple domain Ω with no interface to review the symmetric DDG method for 2nd order elliptic equations. −∆u = f,

x ∈ Ω,

4

u|∂Ω = u0 .

Multiply the Poisson equation with test function v ∈ Vkh , integrate over any element E ∈ Th , and we perform integration by part. After adding the test function numerical flux terms to symmetrize the scheme, we obtain the following symmetric DDG method. We seek the piecewise polynomial solution u ∈ Vkh , such that, Z Z Z Z ∇u · ∇vdx − u cn vds + vf f vdx, for all v ∈ Vkh , (2.2) n JuKds = E

∂E

∂E

E

with, ( \ u cn = ∇u · n = β0u JuK he + {{un }} + β1 he Junn K, JvK ^ vf = ∇v · n = β + {{vn }} + β1 he Jvnn K, n 0v

∀e = ∂E ∩ ∂E 0 .

(2.3)

he

The numerical flux u cn is designed to approximate the normal derivative un = ∇u · n on the discontinuous element edge shared by E and its neighbor element E 0 . The numerical flux formula of u cn involves the solution jump JuK = u|E 0 − u|E , the normal derivative average {{un }} = 12 (∇u|E + ∇u|E 0 ) · n, and the 2nd order normal derivative jump term Junn K. Here we have unn = ∇un · n and Junn K = unn |E 0 − unn |E . Same format numerical flux formula is applied to the test function vf n , which is added into the original DDG formulation, and we have the symmetric DDG method [16]. Since the test function v ∈ Vkh is taken to be non-zero only from the current element E, the term vf n is essentially calculated as the following, when implemented over one element (E ∈ Th ) in (2.2),   (−v) 1 (2.4) + vn + β1 he (−vnn ) |∂E . vf n = β0v he 2 To obtain the primal or global formulation of the symmetric DDG method, we add (2.2) over all elements, correspondingly the test function numerical flux vf n will take the format of (2.3). In the primal formulation, test function values from current element E and its neighbor element E 0 both contribute to the calculation of vf n as in (2.3). We take β0u = β0v , because these two involve the same jump term. We also denote β0 = β0u + β0v . The coefficients, namely (β0 , β1 ) in the numerical flux formula, need to be suitably chosen to stabilize the symmetric DDG scheme. In [30], we show that any (β0 , β1 ) coefficients pair that satisfies the following quadratic form inequality (k is the degree of approximation polynomial space),   2 2 2 k 2 (k 2 − 1) k 2 2 k (k − 1) β0 > 4 (β1 ) − β1 + , 3 2 4 leads to an admissible numerical flux and stabilizes the symmetric DDG method. Summing up the scheme formulation (2.2) over all elements, the primal or global formulation of the symmetric DDG method for Poisson equation can be shown bilinear and symmetric. Correspondingly we can prove the coercivity and bounded-ness of the bilinear form, and we have the stability and energy norm error estimate. We refer to [16] for details. To apply Dirichlet type boundary condition of (1.2), i.e. on the edge e ∈ EhD ⊂ ∂Ω, we have, u cn = β0

JuK + un |e with JuK = u0 − u|e , and vf n = vn |e . he

(2.5)

For piecewise constant (k = 0) and linear (k = 1) approximations, the second derivative jump term Junn K degenerates to zero and has no contribution to the calculation of the numerical flux (2.3). Thus, the symmetric DDG method degenerates to the IPDG method with low order approximations. For higher order approximations (k ≥ 2) and with the second order derivative jump term Junn K included in (2.3), the DDG methods are shown to be more flexible and have several advantages over the IPDG method and the LDG method [10]. The DDG polynomial solutions have been proven to satisfy the strict maximum principle with at least third order of accuracy in [7]. We can only obtain up to second order maximum principle satisfying DG solution for the IPDG and LDG methods. In [33], the DDG solutions are shown to be super convergent to the solution’s spatial derivative. On the other hand, we observe no such super convergence result for the IPDG method. With the concept of u cn that approximates the solution’s normal derivative un = ∇u · n, we can incorporate the interface jump conditions (1.3)-(1.4) into its definition of (2.3), specifically on element edges that overlap with the interface. In the following section, we present the symmetric DDG method for the elliptic interface problem to uniformly accommodate the zero or non-zero solution jump and flux jump interface conditions.

5

2.2

Symmetric DDG method for the interface problem

Now we study the scheme formulation of symmetric DDG method for elliptic interface problem (1.1)-(1.4) with body-fitted mesh (Figure 2). Again, we multiply equation (1.1) with test function v ∈ Vhk , integrate over any element E ∈ Th , and we perform integration by parts. After adding the test function numerical flux terms, we obtain the following symmetric DDG method for (1.1)-(1.4). We seek piecewise DG polynomial solution u ∈ Vkh , such that, Z Z Z Z φ α∇u · ∇vdx − αu dn vds + αv gn JuK ds = f vdx, for all v ∈ Vkh . (2.6) E

∂E

∂E

E

The above scheme formulation is almost identical to the symmetric DDG method (2.2) for the Poisson equation, except the new jump term notion of JuKφ . In this article, we consider body-fitted mesh, thus a triangular element E sits either inside Ω− or Ω+ . In a word, we have E ∈ Th− or E ∈ Th+ . For volume integration terms in (2.6), it is straight forward to take α = α± and f = f ± correspondingly, based on E ∈ Th+ or E ∈ Th− . Now we only need to discuss numerical flux definitions of αu dn on the edges. Recall that we classify all triangular element edges into three categories of Eh = EhI ∪ EhD ∪ EhΓ . For an interior edge (∂E ∈ EhI ), the numerical flux αu dn will be just as (2.3) for Poisson equation. For a boundary edge (∂E ∈ EhD ), the numerical flux will be similar to (2.5). In this section, we mainly focus on the definition of αu dn on the curved edge (∂E ∈ EhΓ ) that overlaps with the curved interface. Since body-fitted mesh is considered, around the interface, curved triangular element (E ∈ ThΓ ) is applied in the symmetric DDG scheme formulation (2.6). Over an interior edge e ∈ EhI , which neither falls on the interface Γ nor on the domain boundary ∂Ω, the numerical fluxes of αu dn and αv gn are taken just as those in (2.3) and we have,    JuK   + { {u } } + β h Ju K , αu d = α β  n 1 e nn n 0      he 1 ∀e ∈ EhI . (2.7) vn − β1 he vnn |e , αv gn = α   2     JuKφ = JuK, Around an interior edge e ∈ EhI , the diffusion coefficient α = α(x) is smooth and continuous, thus we have αu dn = αc un and αv gn = αf vn . Here we only write down the non-zero part of αv gn as in (2.4), since the symmetric DDG method is considered over a single element. Again, the numerical flux u cn quantity, which is evaluated on the edge e = E ∩ E 0 with E 0 as the neighbor element, involves the solution jump JuK = u|E 0 − u|E , the normal derivative average {{un }} = (un |E + un |E 0 )/2, and the second order normal derivative jump Junn K = unn |E 0 − unn |E . Here n is the outward-pointing unit normal vector from E to E 0 . The notation of JuKφ is introduced to those elements that interact with the interface, with which we incorporate the non-zero solution jump condition of (1.3). For an interior edge e ∈ EhI , we simply have JuKφ = JuK. Next we discuss how to define αu dn and αv gn on element edge that overlaps with the interface Γ. Suppose we have E ∈ ThΓ . If the edge ∂E does not fall on the interface, we follow (2.7) on the numerical flux definition. Now we consider the edge ∂E ∩ Γ 6= ∅. To better understand how each jump interface condition is weakly enforced, we discuss two cases. For each case, the problem is associated with non-zero flux jump condition ψ 6= 0 of (1.4). But one case is coupled with zero solution jump φ = 0 of (1.3), and the other coupled with non-zero solution jump φ 6= 0 of (1.3). Since the exact solution satisfies non-zero flux jump condition Jαun K = α− ∇u− · n− + α+ ∇u+ · n+ = ψ, it does not make sense to define one unique value of αu dn on e = E − ∩ E + with E − ⊂ Ω− and E + ⊂ Ω+ . In a word, we assign two separate values of αu dn on the edge e ⊂ Γ. This is different to the regular DG method, in which a unique numerical flux value is defined on the element edge. Thus we maintain the conservative property of a DG method. Now we choose to take αu dn = α \ un− , when being evaluated from the side of E − ⊂ Ω− , and αu dn = αu [ n+ , when being evaluated from the side of E + ⊂ Ω+ . With n+ = −n− , the difference of these two values exactly matches the given jump condition Jαun K = ψ. The flux jump interface condition is thus enforced in the weak sense and is consistent to the weak formulation of the interface problem. Since the flux is already discontinuous over the edge or the interface, we will not include second order normal derivative jump term Junn K into the definition of αu dn . Here we take α = α− or α = α+ correspondingly, given the value being evaluated from E − + or E respectively.

6

Similar idea applies to the non-zero solution jump condition φ 6= 0, where we assign the information of φ into the definition of the solution jump JuKφ whenever being applied across the interface. Thus the solution jump interface condition is also weakly enforced. In summary, we have the following definition for αu dn over an interface edge, ∂E ∈ EhΓ . 1. non-zero flux jump ψ 6= 0 but zero solution jump φ = 0 On the interface edge e = ∂E − ∩ Γ with E − ⊂ Ω− , we define,  ψ JuK   α + {{αun− }} + , \ un− = β0 max{α− , α+ } he 2 1   αv ] = αv − , n− 2 n

∀e ∈ EhΓ .

(2.8)

Here we have αu dn = α \ un− , which is being applied from the side of E − ⊂ Ω− . Notice e = ∂E − ∩ ∂E + , + + with E ∈ Ω as the neighbor element of E − . In (2.8), we have JuK = u|E + − u|E − and {{αun− }} = 1 − + − − + 2 (α ∇u|E + α ∇u|E ) · n . On the interface edge e = ∂E + ∩ Γ with E + ⊂ Ω+ , we define,  JuK ψ   αu [ = β0 max{α− , α+ } + {{αun+ }} + , n+ he 2 1   αv + + = , ] αv n 2 n

∀e ∈ EhΓ .

(2.9)

+ ⊂ Ω+ . Again we have e = Here we have αu dn = αu [ n+ , which is being applied from the side of E − + − − + ∂E ∩ ∂E , and E ⊂ Ω is the neighbor element of E . Thus the jump term is evaluated as JuK = u|E − − u|E + and the average term is evaluated as {{αun+ }} = 12 (α− ∇u|E − + α+ ∇u|E + ) · n+ in (2.9). With the outward unit vector n+ = −n− , the difference of (2.8) and (2.9) in the global formulation exactly matches the given flux jump condition of (1.4). For this case with zero solution jump φ = 0, we have JuKφ = JuK being applied in the symmetric DDG scheme formulation of (2.6).

2. non-zero flux jump ψ 6= 0 and non-zero solution jump φ 6= 0 On the interface edge e = ∂E − ∩ Γ with E − ⊂ Ω− , we define,  φ ψ  − + JuK  − α \ u = β max{α , α } + {{αun− }} + ,  0 n   he 2 1 = αv − , αv ] n−   2 n    JuKφ = JuK + φ,

∀e ∈ EhΓ .

(2.10)

We see the only difference between (2.10) and (2.8) is the term JuKφ = JuK + φ with JuK = u|E + − u|E − . Notation JuKφ is introduced to weakly enforce the solution jump condition (1.3). On the interface edge e = ∂E + ∩ Γ with E + ⊂ Ω+ , we define,  φ ψ  − + JuK  + αu [ = β max{α , α } + {{αun+ }} + ,  0 n   he 2 1 + + αv ] = αv , n   2 n    JuKφ = JuK − φ,

∀e ∈ EhΓ .

(2.11)

When being evaluated on e = ∂E − ∩ E + , from the side of E + , the jump term is JuK = u|E − − u|E + . Thus we define JuKφ = JuK − φ to incorporate the jump condition (1.3). In the symmetric DDG scheme formulation (2.6), we follow (2.10) and (2.11) on the definition of JuKφ . Remark 2.1. To implement the symmetric DDG method (2.6), we simply apply numerical fluxes (2.10) and (2.11) for an interface edge, ∂E ∈ EhΓ . These numerical fluxes degenerate or are consistent to (2.8) and (2.9), when we have zero solution jump (φ = 0) interface condition. 7

D

A

C

B

Figure 3: Left: volume integration over curved triangular element. Right: sample curved region To apply Dirichlet type boundary condition of (1.2), i.e. on the edge e ∈ EhD ⊂ ∂Ω, we have,   JuK αu dn = α β0 + un |e with JuK = u0 − u|e , and αv gn = αvn |e . he

(2.12)

To obtain the primal or global formulation, we sum up (2.6) over all elements E ∈ Th . After applying numerical flux (2.7) to interior edges EhI , numerical flux (2.10) and (2.11) to interface edges EhΓ , and numerical flux (2.12) to domain boundary edges EhD , we obtain the primal formulation of the symmetric DDG method for (1.1)-(1.4) as the following. Seek numerical solution u ∈ Vkh , such that, Bh (u, v) = F (v),

∀v ∈ Vkh .

The bilinear form Bh (w, v) is listed below as,  X Z  β0 X Z X Z β0 wv − vn w − wn v ds Bh (w, v) = α∇w · ∇vdx + α α JwKJvKds + he e e he D I E∈Th E e∈Eh e∈Eh XZ + α {({{wn }} + β1 he Jwnn K) JvK + ({{vn }} + β1 he Jvnn K) JwK} ds I e∈Eh

+

e

(2.13)

(2.14)

 X Z  β0 − + max{α , α }JwKJvK + {{αwn }}JvK + {{αvn }}JwK ds, he e Γ

e∈Eh

with the right hand side F (v) given as,  X Z  β0 − + F (v) = max{α , α }φJvK + {{αvn }}φ + ψ{{v}} ds he e Γ e∈Eh  Z X Z  β0 + α u0 v − vn u0 ds + f vdx. he e Ω D

(2.15)

e∈Eh

Remark 2.2. The interface solution jump (1.3) and flux jump (1.4) conditions are both weakly enforced. The symmetric DDG scheme formulation (2.13) is consistent to the symmetric DDG method for smooth problems with zero solution jump φ = 0 and zero flux jump ψ = 0, or smooth problems with no interface. Remark 2.3. Similar ideas for numerical fluxes (2.10)-(2.11) showed up in Local DG methods solving mixed boundary problems [5] and boundary element methods [6] and [17]. 8

P3

P2 P1

Figure 4: Sample curved mesh/element generation. Left: interface description; Middle: initial curved elements; Right: refined curved elements. Next, we discuss issues related to the implementation of symmetric DDG method. Under a general setting, the interface Γ is a curve and curved triangular elements are applied in the symmetric DDG method around the interface. Now we discuss in details how the volume integration over a curved element (E ∈ ThΓ ), the left figure in Figure 3, is approximated. We also discuss how we generate the curved mesh through Matlab. Suppose E is a curved element with P1 , P2 , P3 as its three vertices, and P2 P3 is the curved edge (left figure in Figure 3). For the illustrated example with a convex curved edge, the curved element E is divided into two subdomains with E1 denoting the straight triangle part and E2 as the other curved part. In a word, we have Z Z Z g(x, y)dxdy = g(x, y)dxdy + g(x, y)dxdy. E

E1

E2

The volume integration on E1 is over a straight triangle. We apply the 13 points Gaussian quadrature rule [11], which has 7th order of accuracy, to approximate the volume integration over E1 . To approximate the volume integration over the curved region, we use the right figure in Figure 3 to illustrate the quadrature rule. Suppose the curved region is enclosed by a straight line segment AB and a curved segment ADB. Denote the center of line segment AB as point C. We have point D taken on the curved segment ADB such that line CD is perpendicular to line AB. Now we set up a new coordinate system with point C as the origin, line AB as the new x-axis and line CD as the new y-axis. If the curved segment ADB is not convex, line CD points to the negative y-axis. Suppose we have (η, ξ) denoting the new coordinates, and the curved line equation of ADB is given with ξ = ψ(η). Letter J denotes the Jacobian matrix. We have, Z Z Z 1 Z ψ(η) g(x, y)dxdy = g(η, ξ)|J|dηdξ = dη g(η, ξ)|J|dξ. −1

ACBDA

0

Then we apply the dimension by dimension Gaussian quadrature rule to approximate the above volume integration. Suppose we use ηi denoting the ith Gaussian quadrature point over [−1, 1] and wi as its weight, we have, Z X Z ψ(ηi ) g(x, y)dxdy ≈ wi g(ηi , ξ)|J|dξ. ACBDA

i

0

R ψ(η ) Finally, a similar one-dimensional Gaussian quadrature rule is applied to approximate 0 i g(ηi , ξ)|J|dξ. Notice that we need to pay attention to the concavity property of the interface curve Γ around the Gaussian point (ηi , ξi ). Here ηi is the ith Gaussian point over [−1, 1]. If the curve segment around (ηi , ξi ) is concave, R0 R ψ(η ) the integration should be ψ(ηi ) g(ηi , ξ)|J|dξ, instead of the convex case of 0 i g(ηi , ξ)|J|dξ. Regarding the curved mesh/element generation, we apply the ”pdetool” tool box in Matlab for mesh partition. For example, we perform the following three steps to generate the mesh of Figure 4. 9

Step 1. Apply the ”pdegeom” function to define the geometry of the domain. (See Figure 4). The domain is represented by parameterized edge segments. Notice the edge segments cannot overlap. Each pair of the edge segments must intersect with the endpoint. For example, the domain in Figure 4 include three parameterized edge segments, which are straight lines P1 P2 , P3 P1 , and the curved line segment P2 P3 ; Step 2. Call ”initmesh” function to draw the initial mesh; Step 3. Call ”refinemesh” function in the next step to obtain the next level refined mesh partition.

3

Continuity, coercivity and well-posedness

In this section and the following, we assume the interface Γ being in polygonal shape and we have straight triangular mesh partition Th of the domain Ω. Now we study the well-posedness of the variational problem of the symmetric DDG method. We first discuss the boundedness and stability of the bilinear form Bh (·, ·) of (2.14). We define the energy norm of v ∈ Vkh as:  |||v|||h := 

1/2 X Z v2 X Z JvK2 ds + ds . ∇v · ∇vdx + E e he e he D I

X Z

E∈Th

(3.1)

e∈Eh

e∈Eh

Below we list the trace inequality and inverse inequality of the solution space Vkh ([29]). Lemma 3.1. Consider any triangular element E ∈ Th and v ∈ H s (E) with s ≥ 2. On the edge e of the triangle E, there exists a positive constants Cg independent of hE such that, −1/2

kvkL2 (e) ≤ Cg hE

(kvkL2 (E) + hE |v|H 1 (E) ), −1/2

kvn kL2 (e) = k∇v · nkL2 (e) ≤ Cg hE

(|v|H 1 (E) + hE |v|H 2 (E) ).

Lemma 3.2. For any element E ∈ Th and v ∈ Pk (E), there exist positive constants Ct and Ci which are independent of hE such that, −1/2 kvkL2 (e) ≤ Ct hE kvkL2 (E) , ∀e ∈ ∂E, −1/2

kvn kL2 (e) = k∇v · nkL2 (e) ≤ Ct hE

|v|H j (E) ≤ Ci h−j+1 |v|H 1 (E) , E

|v|H 1 (E) ,

∀e ∈ ∂E,

0 ≤ j ≤ k.

To study the continuity and coercivity of the bilinear form (2.14) and error estimates of the symmetric DDG method solution, we first prove the following inequalities. Lemma 3.3. Let e = ∂E ∩ ∂E 0 where E and E 0 are two adjacent elements of Th . Assume that w ∈ H s (E ∪ E 0 ), s ≥ 3 and v ∈ Vkh , we have the following estimates: Z Z 1/2 JwK2 {{α∇v · n}}JwKds ≤ 1 Ct αu |v|H 1 (E∪E 0 ) ds , 2 e he e Z Z 1/2 JwK2 β1 he Jαvnn KJwKds ≤ 2β1 Ct Ci αu |v|H 1 (E∪E 0 ) ds , e e he Z Z 1/2
JvK2 {{α∇w · n}}JvKds ≤ 1 Cg αu |w|H 1 (E∪E 0 ) + h|w|H 2 (E∪E 0 ) ds , 2 e e he Z Z 1/2
JvK2 β1 he Jαwnn KJvKds ≤ 3β1 Cg αu h |w|H 2 (E∪E 0 ) + h|w|H 3 (E∪E 0 ) ds . e e he We obtain similar results on the boundary edges e ∈ EhD .

10

Proof. Given w ∈ H s (E ∪ E 0 ) and v ∈ Vkh , let’s use w, v and w0 , v 0 to denote the functions defined in element E and element E 0 , respectively. Notice that e = ∂E ∩ ∂E 0 and the fact he ≤ min{hE , hE 0 }, with Lemma 3.2 we have, Z {{α∇v · n}}JwKds

1/2 Z
1 p JwK2 0 0 αu he k∇v · n kL2 (e) + k∇v · nkL2 (e) ds 2 e he 1/2 Z
1 JwK2 0 . Ct αu |v|H 1 (E) + |v |H 1 (E 0 ) ds 2 e he

≤

e

≤

Recall αu is the upper bound of the diffusion coefficient α of (1.5). By applying Lemma 3.1 and the notation h = maxE∈Th hE , we have the following estimate on the third inequality, Z {{α∇w · n}}JvKds

1/2 JvK2 ds e he 1/2 Z
JvK2 ds . |w|H 1 (E) + h|w|H 2 (E) + |w0 |H 1 (E 0 ) + h|w0 |H 2 (E 0 ) e he

≤


1 p αu he k∇w0 · n0 kL2 (e) + k∇w · nkL2 (e) 2

≤

1 Cg αu 2

e

Z

For the forth inequality, we first have the fact of n = (n1 , n2 ) and |n1 | ≤ 1, |n2 | ≤ 1. After applying Lemma 3.2 we have, Z e



2 Z X 2 2 X 2 p 2 X X ∂2w

∂ w

he |wnn JvK|ds ≤ ∂xi ∂xj |JvK|ds ≤ ∂xi ∂xj e i=1 j=1

i=1 j=1

L2 (e)

i=1 j=1

2

2 2 X X

∂ w

≤ Cg

∂xi ∂xj

L2 (E)

≤ 3Cg |w|H 2 (E) + hE |w|H 3 (E)

2 ∂ w + hE ∂xi ∂xj


Z e

Z

! Z H 1 (E)

e

e

1/2 JvK2 ds he

1/2 JvK2 ds he

1/2 JvK2 ds . he

Finally we obtain, Z Z β1 he Jαwnn KJvKds ≤ β1 he αu (|wnn | + |wn0 0 n0 |)|JvK|ds e

e

≤ 3β1 Cg αu h |w|H 2 (E) + hE |w|H 3 (E) + |w0 |H 2 (E 0 ) + hE 0 |w0 |H 3 (E 0 )




Z e

1/2 JvK2 ds . he

Similar discussion applies to the second inequality in this Lemma. Now we are ready to establish the continuity and coercivity of the bilinear form of (2.14). Theorem 3.1. With admissible parameter pair (β0 , β1 ) taken in the numerical flux, there exist positive constants Cs and Cb for the symmetric DDG method (2.13) such that for any w, v ∈ Vkh , we have |Bh (w, v)| ≤ Cb |||w|||h |||v|||h ,

(3.2)

Bh (v, v) ≥ Cs |||v|||2h .

(3.3)

Proof. The proof is similar to Theorem 3.1 in [16]. Theorem 3.2. With admissible parameter pair (β0 , β1 ) taken in the numerical flux, problem (2.13) exists a unique solution u ∈ Vkh . Proof. Since (2.13) is a linear problem in finite dimensional space Vkh , existence is equivalent to uniqueness. Suppose that there are two numerical solutions u1 and u2 , then we have the difference w = u1 − u2 satisfying Bh (w, w) = 0. By the coercivity result of (3.3), we have |||w|||h = 0 which directly implies that u1 = u2 . 11

Theorem 3.3. Let’s denote uex ∈ H s (Ω+ ∪ Ω− ), s ≥ 3 as the exact solution of the elliptic interface problem (1.1)-(1.4). Assume that admissible parameter pair (β0 , β1 ) is chosen in the numerical flux (2.7) and (2.10)(2.10). Denoting u ∈ Vkh as the symmetric DDG method (2.13) solution, we have the following orthogonality equation Bh (uex − u, v) = 0, ∀v ∈ Vkh . (3.4) Proof. Since the symmetric DDG solution u ∈ Vkh satisfies the bilinear form (2.13), we only need to prove the exact solution uex satisfying Bh (uex , v) = F (v), ∀v ∈ Vkh . For any v being a piecewise polynomial in Vkh , we multiply (1.1) with v, integrate over an element E ∈ Th , after integration by parts we have Z Z Z α∇uex · ∇vdx − α(uex )n vds = f vdx. E

∂E

E

Here (uex )n = ∇uex · n is the outward normal derivative along the edge ∂E. We sum the above equation over all elements E ∈ Th and formally we have, X Z X Z X Z α∇uex · ∇vdx − α(uex )n vds = f vdx. E

E∈Th

E∈Th

∂E

E∈Th

E

Notice the edges are classified into three groups with Eh = EhI ∪ EhD ∪ EhΓ , thus we have, X Z XZ XZ X Z + + + − − − α(uex )n vds = − α(uex )n JvKds + α (uex )n+ v + α (uex )n− v ds + α(uex )n vds.

E∈Th

∂E

I e∈Eh

e

Γ e∈Eh

e

D e∈Eh

e

(3.5) Here we use the fact that the diffusion coefficient α is continuous across an interior edge e ∈ EhI . The exact solution’s normal derivative is also continuous across interior edges e ∈ EhI , thus we have {{α(uex )n }} = α(uex )n ,

Juex K = 0, Jα(uex )n K = 0, β1 he Jα(uex )nn K = 0,

∀e ∈ EhI .

(3.6)

The above arguments hold true since uex ∈ H 3 (Ω+ ∪ Ω− ) and α+ , α− are continuous function in Ω+ , Ω− , respectively. Over the edges on the interface e ∈ EhΓ we have, + − − − (u+ ex )n+ v + α (uex )n− v

1 + + 1 + + − − + − + + − (α ∇uex + α− ∇u− ex ) · n (v − v ) + (α (uex )n+ + α (uex )n− )(v + v ) 2 2 = −{{α(uex )n }}JvK + ψ{{v}}. (3.7) =

Here we need the interface condition (1.4), the notion of n+ = −n− and the jump and average definition of v ∈ Vkh over an edge. Collecting the details of (3.5), (3.6) and (3.7), we have the exact solution satisfying the following weak formulation, Z X Z X Z X Z XZ α∇uex · ∇vdx + {{α(uex )n }}JvKds − α(uex )n vds = f vdx + ψ{{v}}. E∈Th

E

I ∪E Γ e∈Eh h

e

D e∈Eh

e

Ω

Γ e∈Eh

e

Now combine the facts of (3.6) over interior edges, the solution jump interface condition (1.3) and the Dirichlet boundary condition (1.2), we have the exact solution satisfies the bilinear form Bh (·, ·) of (2.14) as, Bh (uex , v) = F (v),

4

∀v ∈ Vkh .

Error estimate

In this section, we carry out energy norm and L2 norm error estimates for the symmetric DDG method (2.13). We first list the approximation properties of the solution space Vkh , then establish the energy norm error estimate. At the end of this section, we obtain the optimal error estimate of the symmetric DDG solution under L2 norm. 12

Lemma 4.1. (see Theorem 2.6 in[29]) Consider any triangular element E ∈ Th and v ∈ H s (E) with s ≥ 1. Denote Ih v as the kth degree polynomial L2 projection of v in Pk (E). There exists a constant CI independent of v and hE such that min(k+1,s)−q kv − Ih vkH q (E) ≤ CI hE |v|H s (E) , 0 ≤ q ≤ s. (4.1) Theorem 4.1. Assume that admissible parameter pair (β0 , β1 ) is chosen in the numerical flux (2.7) and (2.10)-(2.10). We have uex ∈ H s (Ω+ ∪ Ω− ), s ≥ 3, as the exact solution of the interface problem (1.1)-(1.4) and u ∈ Vkh as the symmetric DDG solution of (2.13), then we have |||uex − u|||h ≤ Chmin(k+1,s)−1 |uex |H s (Ω− ∪Ω+ ) ,

(4.2)

where the constant C depends on (β0 , β1 ) and diffusion coefficient α, is independent of mesh size h and the exact solution uex . Proof. Denote Ih uex ∈ Vkh as the L2 projection of the exact solution uex . With the coercivity of the bilinear form (Theorem 3.1) and the error equation (Theorem 3.3), we have Bh (uex − Ih uex , u − Ih uex ) = Bh (u − Ih uex , u − Ih uex ) ≥ Cs |||u − Ih uex |||2h . Now introduce Φ = uex − Ih uex and ξ = u − Ih uex . We have ξ ∈ Vkh and consider a general edge e = ∂E ∩ ∂E 0 with elements E ∈ Th and E 0 ∈ Th . Again we have αu denoting the upper bound of the diffusion coefficient α over the whole domain Ω. With Cauchy-Schwartz inequality and applying the results of Lemma 3.3 and approximation results of Lemma 4.1 we have, Z α∇Φ · ∇ξdx ≤ αu |Φ|H 1 (E) |ξ|H 1 (E) ≤ αu CI hmin(k+1,s)−1 |uex |H s (E) |ξ|H 1 (E) , E

Z Z 1/2 JξK2 {{αΦn }}JξKds ≤ 1 Cg αu CI hmin(k+1,s)−1 |uex |H s (E∪E 0 ) ds , 2 e he e Z Z 1/2 JξK2 β1 he JαΦnn KJξKds ≤ 3β1 αu Cg CI hmin(k+1,s)−1 |uex |H s (E∪E 0 ) ds . he e

e

Similarly applying Lemma 3.3 and Cauchy-Schwartz inequality we obtain the following estimates across any edge e = ∂E ∩ ∂E 0 , Z Z 1/2 JΦK2 β1 he Jαξnn KJΦKds ≤ 2β1 Ct Ci αu |ξ|H 1 (E∪E 0 ) ds , e he e Z Z 1/2 Z 1/2 JξK2 JΦK2 α β0 JΦKJξKds ≤ αu β0 ds ds , e he e he e he Z Z 1/2 JΦK2 {{α∇ξ · n}}JΦKds ≤ 1 Ct αu |ξ|H 1 (E∪E 0 ) ds . 2 e he e Since we have a regular mesh partition Th , there exists a positive constant Ch such that for any E ∈ Th we have hE ≤ Ch he with Ch independent of E. With the trace inequality of Lemma 3.1 and approximation result of Lemma 4.1, we have Z Z 2 Φ + Φ02 JΦK2 ds ≤ 2 ds he e he e
2
2 2Cg2 2Cg2 ≤ kΦkL2 (E) + hE |Φ|H 1 (E) + kΦ0 kL2 (E 0 ) + hE 0 |Φ0 |H 1 (E 0 ) he hE he hE 0 ≤ 8Ch Cg2 CI2 h2(min(k+1,s)−1) |uex |2H s (E∪E 0 ) . We obtain similar estimates on the boundary edge e ∈ EhD . Now we plug in above inequalities into the bilinear form Bh (uex − Ih uex , u − Ih uex ) of (2.14) and sum up over all elements, we have Bh (uex − Ih uex , u − Ih uex ) ≤ C1 CI hmin(k+1,s)−1 |uex |H s (Ω+ ∪Ω− ) |||u − Ih uex |||h ,

13

√ √ √ √ √ where C1 = αu max{1 + (12 2β1 Ci + 3 2)Ct Cg Ch , Cg (1.5 + 9β1 + 6 2β0 Ch )}. With the coercivity of the bilinear form we directly obtain, |||u − Ih uex |||h ≤

C1 CI min(k+1,s)−1 h |uex |H s (Ω+ ∪Ω− ) . Cs

By the definition of the energy norm ||| · |||h and the L2 projection Ih uex ∈ Vkh , we have the following projection error estimate, q |||uex − Ih uex |||h = |||Φ|||h ≤ CI 1 + 24Ch Cg2 hmin(k+1,s)−1 |uex |H s (Ω− ∪Ω+ ) . Finally, a simple triangle inequality leads to the optimal estimate of C1 q |||uex − u|||h ≤ CI ( + 1 + 24Ch Cg2 )hmin(k+1,s)−1 |uex |H s (Ω− ∪Ω+ ) . Cs Next, we carry out the optimal L2 norm error estimate for the symmetric DDG solution (2.13) with the standard duality argument. We define auxiliary function Ψ to be the solution of the following adjoint problem,  ± ± ±   in Ω± ,  −∇ · (α ∇Ψ ) = (uex − u) := χ ,  (4.3) Ψ = 0, on ∂Ω,     Ψ− − Ψ+ = 0, α− ∇Ψ− · n− + α+ ∇Ψ+ · n+ = 0, on Γ. With Ψ as the solution of the adjoint problem and the regularity result [8], we have Ψ ∈ H 1 (Ω) ∩ H 2 (Ω+ ∪ Ω− ) such that kΨkH 2 (Ω+ ∪Ω− ) ≤ Cr kuex − ukL2 (Ω− ∪Ω+ ) . (4.4) Theorem 4.2. Let uex ∈ H s (Ω+ ∪ Ω− ), s ≥ 3 be the exact solution of the interface problem (1.1)-(1.4) and u ∈ Vkh denote the symmetric DDG method solution of (2.13). Let Ψ be the solution of adjoint problem (4.3) with regularity Ψ ∈ H q (E), q > 5/2 over every element E ∈ Th . There exists constant C exists, such that kuex − ukL2 (Ω) ≤ Chmin(k+1,s) |uex |H s (Ω− ∪Ω+ ) ,

(4.5)

where the constant C depends on (β0 , β1 ) and diffusion coefficient α but is independent of mesh size h and the exact solution uex . Proof. According to the regularity result (4.4) of the solution Ψ and the fact that the diffusion coefficient α is smooth inside Ω+ and Ω− , over the interior edges e ∈ EhI we have J∇Ψ · nK = 0,

JΨK = 0,

{{α∇Ψ · n}} = α∇Ψ · n on e ∈ EhI .

(4.6)

Recall that the adjoint problem (4.3) is associated with zero interface conditions of JΨK = 0 and Jα∇Ψ · nK = 0 over Γ and zero boundary condition Ψ|∂Ω = 0. Following the bilinear form (2.14) we formally have, X Z XZ Bh (χ, Ψ) = α∇χ · ∇Ψdx + ({{αΨn }} + β1 he JαΨnn K) JχKds E∈Th

+

E

X Z E∈Th

e

X Z

{{αΨn }}JχKds −

e

Γ e∈Eh

=

I e∈Eh

XZ

α∇χ · ∇Ψdx −

E

αΨn χds

e

D e∈Eh

X Z E∈Th

αΨn χds +

∂E

XZ I e∈Eh

e

β1 he JαΨnn KJχKds.

Now multiple the adjoint problem (4.3) with χ± and integrate over Ω± , have integration by parts over each element E ∈ Th and we obtain, Z 2 kχkL2 (Ω) = − ∇ · (α∇Ψ)χdx Ω+ ∪Ω− Z X X Z = α∇Ψ · ∇χdx − αΨn χds E∈Th

E

= Bh (χ, Ψ) −

E∈Th

XZ I e∈Eh

14

e

∂E

β1 he JαΨnn KJχKds.

We choose Ih Ψ to denote the continuous linear polynomial solution of the adjoint problem (4.3), then we have (see [29]), kΨ − Ih ΨkH 1 (Ω− ∪Ω+ ) ≤ CI h|Ψ|H 2 (Ω− ∪Ω+ ) . (4.7) By the orthogonality equation (3.4) we have Bh (uex − u, Ih Ψ) = 0 since the linear polynomial solution Ih Ψ ∈ Vkh . Now combining the above argument we obtain, XZ β1 he JαΨnn KJχKds. kuex − uk2L2 (Ω) = kχk2L2 (Ω) = Bh (χ, Ψ − Ih Ψ) − e

I e∈Eh

Since Ih Ψ ∈ Vkh is a continuous linear polynomial over Th , we have JIh ΨKe = 0,

and

(Ih Ψ)nn |e = 0,

∀e ∈ EhI ∪ EhΓ .

Furthermore we have the adjoint problem solution Ψ ∈ H 2 (Ω+ ∪ Ω− ) and Ψ is continuous in Ω± and the fact of JΨK|Γ = 0. Also the boundary condition is Ψ|∂Ω = 0, and we have Ih Ψ|∂Ω = 0. Plugging in all these facts, the bilinear form is simplified as X Z X Z Bh (χ, Ψ − Ih Ψ) = α∇(Ψ − Ih Ψ) · ∇χdx + {{α∇(Ψ − Ih Ψ) · n}}JχKds E∈Th

E

−

α∇(Ψ − Ih Ψ) · nχds +

e

D e∈Eh

e

I ∪E Γ e∈Eh h

X Z

XZ I e∈Eh

e

β1 he JαΨnn KJχKds.

Using Cauchy-Schwartz inequality, the approximation result of (4.7) and Lemma 3.3, we obtain the following estimate, X kχk2L2 (Ω) ≤ αu |Ψ − Ih Ψ|H 1 (E) |χ|H 1 (E) E∈Th

1 + Cg αu 2 +Cg αu

X

|Ψ − Ih Ψ|H 1 (E∪E 0 ) + h|Ψ|H 2 (E∪E 0 )




Z e

I ∪E Γ e∈Eh h

X

|Ψ − Ih Ψ|H 1 (E) + h|Ψ|H 2 (E)




Z e

D e∈Eh

1/2 JχK2 ds he

1/2 χ2 ds he

q ≤ αu h CI2 + 4Cg2 (CI + 1)2 |Φ|H 2 (Ω+ ∪Ω− ) |||χ|||h . Applying the regularity (4.4) of the solution Φ and the fact of χ = uex − u, we have q kuex − ukL2 (Ω) ≤ αu CI2 + 4Cg2 (CI + 1)2 Cr h|||uex − u|||. Finally combining with Theorem 4.1, we obtain the optimal L2 norm error estimate.

5

Example

In this section, we present a sequence of numerical examples to validate the high order accuracy of the proposed symmetric DDG method (2.6)-(2.7) with interface numerical fluxes (2.10)-(2.11). All examples have exact solutions available, thus we mainly describe the geometric of the interface Γ and the discontinuous diffusion coefficient α of (1.6). We skip the right hand side terms f + and f − presentations in the elliptic equations (1.1), since they can be easily computed based on the given exact solutons in the examples. The Dirichlet boundary data (1.2) and the two interface jump conditions (1.3)-(1.4) can also be derived from the given exact solutions. The symmetric DDG method degenerates to the IPDG method with low order piecewise constant and linear polynomial approximations (Pk with k ≤ 1), thus we only present high order implementations (Pk with k ≥ 2) in all examples. The mesh generation and the symmetric DDG method implementation are

15

Figure 5: Two domains of Example 5.1. all conducted in the MATLAB environment. In this section, we use notation x = (x, y) ∈ Ω as the spatial variables of the solution. In general, the interface Γ is a curve (except Example 5.1) and curved elements are applied to implement the symmetric DDG method (2.6). Volume integration over a curved element E ∈ ThΓ is required (refer to §2.2). For a regular straight triangular element E ∈ Th \ ThΓ , both the volume integration and the line integration are calculated out either exactly or approximated with a Gaussian quadrature rule (13 points, up to 7th order of accuracy for volume integration approximation). For the rest of this section, we use uex and uh to denote the exact solutions and the symmetric DDG solutions of the interface problem (1.1)-(1.4). In each example, errors with following notations are measured: keh kL∞ :=

max

x∈Ω− ∪Ω+

|uh − uex |,

keh kL2 := kuh − uex kL2 (Ω+ ∪Ω− ) ,

|eh |H 1 := |uh − uex |H 1 (Ω+ ∪Ω− ) .

We should mention the L∞ norm errors are computed over the Gaussian points (13 points) in each triangular element. Example 5.1. Piecewise constant diffusion coefficient and straight line interface The computational domain is a square with Ω = (0, 1) × (0, 1). The interface is a straight line of Γ : x = b = 0.5. We define Ω− = {(x, y) | x < b, (x, y) ∈ Ω} and Ω+ = {(x, y) | x > b, (x, y) ∈ Ω}, see Figure 5. The diffusion coefficient is two piecewise constants with α = α− in Ω− and α = α+ in Ω+ . We carry out three accuracy tests in this example with zero or non-zero interface jump conditions (1.3)-(1.4). The interface is a straight line, thus the mesh is partitioned long with the interface. Regular straight triangular elements are used in the symmetric DDG method implementations and we have no geometric error introduced in this example. Case I: zero solution jump φ = 0 and zero flux jump ψ = 0 The exact solution is available with, uex (x, y) =

 πx  1 sin (x − b)(y − 0.5)(1 + x2 + y 2 ). α 2

The piecewise diffusion coefficient is taken as (α− , α+ ) = (1, 100). We obtain optimal (k + 1)th order convergence on the L∞ and L2 errors and kth order on the energy norm errors. In Table 1 we list the errors and orders obtained for this case. Case II: zero solution jump φ = 0 and non-zero flux jump ψ 6= 0

16

Table 1: Case I of Example 5.1 with zero solution jump and zero flux jump. k, β0 , β1

h

keh kL∞

2

0.3536

2.0500e-3

β0 = 4.5

0.1768

3.2015e-4

2.68

4.3465e-5

3.07

4.1391e-3

2.05

β1 = 1/40

0.0884

4.4374e-5

2.85

5.2365e-6

3.05

1.0089e-3

2.04

0.0442

5.8262e-6

2.93

6.4176e-7

3.03

2.4861e-4

2.02

0.0221

7.4591e-7

2.97

7.9447e-8

3.01

6.1687e-5

2.01

3

0.3536

6.5600e-5

β0 = 10

0.1768

4.3086e-6

3.92

7.3838e-7

3.96

1.1177e-4

3.05

β1 = 1/40

0.0884

2.7705e-7

3.96

4.7057e-8

3.97

1.3665e-5

3.03

0.0442

1.8958e-8

3.87

2.9758e-9

3.98

1.6872e-6

3.02

0.0221

1.2370e-9

3.94

1.872e-10

3.99

2.0955e-7

3.01

order

keh kL2

order

3.6622e-4

|eh |H 1

order

1.7110e-2

1.1470e-5

9.2320e-4

Table 2: Case II of Example 5.1 with zero solution jump but non-zero flux jump k, β0 , β1

h

keh kL∞

2

0.3536

1.4079e-2

β0 = 4.5

0.1768

2.0168e-3

2.80

1.1660e-4

3.36

8.3194e-3

2.32

β1 = 1/40

0.0884

2.7304e-4

2.88

1.1449e-5

3.35

1.7276e-3

2.27

0.0442

3.5609e-5

2.94

1.1687e-6

3.29

3.7869e-4

2.19

0.0221

4.5491e-6

2.97

1.2578e-7

3.22

8.7168e-5

2.12

3

0.3536

4.6874e-4

β0 = 4.5

0.1768

2.9175e-5

4.01

2.9725e-6

4.68

3.1529e-4

3.62

β1 = 1/40

0.0884

1.8681e-6

3.97

1.2988e-7

4.52

2.9335e-5

3.43

0.0442

1.1847e-7

3.98

6.0956e-9

4.41

2.9921e-6

3.29

0.0221

7.4592e-9

3.99

3.068e-10

4.31

3.2898e-7

3.19

order

keh kL2

order

1.2006e-3

|eh |H 1

order

4.1661e-2

7.6113e-5

3.8714e-3

The exact solution is given with uex (x, y) = sin

 πx  2

(x − b)(y − 0.5)(1 + x2 + y 2 ).

The piecewise diffusion coefficient is taken as (α− , α+ ) = (1, 100). Symmetric DDG method with interface numerical fluxes (2.8)-(2.9) is implemented and optimal (k + 1)th order convergence is obtained on the L∞ and L2 errors, see Table 2. Case III: non-zero solution jump φ 6= 0 and non-zero flux jump ψ 6= 0 The exact solution is available with uex (x, y) = sin

 πx  2

(x − b)(y − 0.5)(1 + x2 + y 2 ) + α.

The diffusion coefficient is taken with (α− , α+ ) = (1, 1000). Numerical fluxes (2.10)-(2.11) are applied on elements with edges falling on the interface. We carry out accuracy test with P2 and P3 polynomial approximations and obtain optimal 3rd and 4th orders of convergence on the L∞ and L2 errors, see Table 3. Example 5.2. Curved interface 17

Table 3: Case III of Example 5.1 with non-zero solution jump and non-zero flux jump k, β0 , β1

h

keh kL∞

2

0.3536

7.1907e-2

β0 = 4.5

0.1768

9.7213e-3

2.89

5.5198e-4

3.66

4.1657e-2

2.64

β1 = 1/40

0.0884

1.2632e-3

2.94

4.5080e-5

3.61

6.8755e-3

2.60

0.0442

1.6100e-4

2.97

3.8241e-6

3.56

1.1804e-3

2.54

0.0221

2.0320e-5

2.99

3.3600e-7

3.51

2.1087e-4

2.48

3

0.3536

3.6816e-3

β0 = 4.5

0.1768

2.2553e-4

4.03

1.2983e-5

4.85

1.3118e-3

3.81

β1 = 1/40

0.0884

1.4019e-5

4.01

4.8829e-7

4.73

1.0171e-4

3.69

0.0442

8.7414e-7

4.00

1.9597e-8

4.64

8.4238e-6

3.59

0.0221

5.4566e-8

4.00

8.2997e-10

4.56

7.3884e-7

3.51

order

keh kL2

order

6.9754e-3

|eh |H 1

order

2.6034e-1

3.7466e-4

1.8456e-2

In this example, we consider a classical circular interface problem with non-zero solution jump φ 6= 0 and non-zero flux jump ψ 6= 0 conditions across the interface. The domain is a square pΩ = (−1, 1) × (−1, 1) with a circular interface Γ = {(x, y) | x2 + y 2 = 14 }. We introduce the notation of r = x2 + y 2 for any (x, y) ∈ Ω. The diffusion coefficient is piecewise defined as,   α− = 2, r ≤ 1 , 2 α=  α+ = 10, r > 1 . 2

The inhomogeneous right side is given with,   −4α− , r ≤ 12 , f (x, y) =  −8r2 − 4, r > 1 . 2 The exact solution is available with,    r2 − 1,  4  uex (x, y) = 1 1 1 r 2  − + + r /α+ ,  − 4 32α+ 4α+ 2

r ≤ 12 , r > 21 .

In this example, we carry out three tests with the proposed symmetric DDG method (2.6)-(2.7). The first two tests are based on straight triangular mesh, which geometrically approximates the circular interface, and the volume integrations are calculated on straight triangles. The third test is applied on curved triangular elements and the volume integration is calculated on the curved element. In the very first test, the interface jump conditions are still applied on the curved edge of the triangle with numerical fluxes (2.10)-(2.11). And the element edge line integral is calculated on the circular curve. On the other hand, the volume integration is calculated over the approximated straight triangle. Different from [34], we observe one order loss with P2 and P3 polynomials approximations, see Table 4. In the second test, the interface can be considered to be first approximated with a polygon with edges in the length scale of the mesh size h. In a word, the circular interface is approximated with the collection of straight triangular edges. The polygon interface sits totally inside the circle, thus we keep the exact solution expression within the circle as the inside solution for the polygon interface. We then extend the outside of the circle exact solution to the polygon interface. With the two piecewise defined modified exact solution, we can calculate the solution jump and flux jump interface conditions across the polygon, with which the numerical fluxes can be defined suitably. We carry out accuracy test with P2 and P3 polynomials and obtain optimal (k + 1)th order convergence on the L∞ and L2 errors, see Table 5. Generally the jump conditions on the polygon should be approximated or extended from the given jump conditions on the curved interface. 18

Figure 6: Example 5.2 with mesh h = 0.0658 (left) and symmetric DDG solution(right) of P2 approximation.

Table 4: First test of Example 5.2 with order loss k, β0 , β1

h

keh kL∞

2

0.2150

1.0088e-2

β0 = 18

0.1075

2.3178e-3

2.12

1.3337e-3

2.19

2.2570e-2

1.56

β1 = 1/40

0.0538

5.5709e-4

2.06

3.0933e-4

2.11

7.7702e-3

1.54

0.0269

1.3646e-4

2.03

7.4308e-5

2.06

2.7062e-3

1.52

3

0.2150

1.5458e-2

β0 = 20

0.1075

2.9945e-3

2.37

1.5606e-3

2.39

3.5531e-2

1.59

β1 = 1/40

0.0538

6.4940e-4

2.21

3.3617e-4

2.21

1.2123e-2

1.55

order

keh kL2

order

6.1004e-3

|eh |H 1

order

6.6651e-2

8.1531e-3

1.0696e-1

Table 5: Second test of Example 5.2 on the approximated polygonal interface k, β0 , β1

h

keh kL∞

2

0.2150

7.6002e-5

β0 = 18

0.1075

1.0069e-5

2.92

4.1754e-6

2.95

4.3491e-4

2.01

β1 = 1/40

0.0538

1.2783e-6

2.98

5.3226e-7

2.97

1.0847e-4

2.00

3

0.2150

1.6054e-6

β0 = 20

0.1075

9.2994e-8

4.11

2.7265e-8

4.15

4.9740e-6

3.31

β1 = 1/40

0.0538

5.7761e-9

4.01

1.6819e-9

4.02

5.7630e-7

3.11

order

keh kL2

order

3.2253e-5

order

1.7557e-3

4.8440e-7

19

|eh |H 1

4.9366e-5

Table 6: Third test of Example 5.2 on curved triangular elements k, β0 , β1

h

keh kL∞

2

0.2420

7.6011e-5

β0 = 18

0.1282

1.0069e-5

2.92

4.1737e-6

2.95

4.3495e-4

2.01

β1 = 1/40

0.0658

1.2783e-6

2.98

5.3223e-7

2.97

1.0847e-4

2.00

3

0.2420

1.6054e-6

β0 = 20

0.1282

9.2994e-8

4.11

2.7264e-8

4.15

4.9750e-6

3.31

β1 = 1/40

0.0658

5.7760e-9

4.01

1.6819e-9

4.02

5.7633e-7

3.11

4

0.2420

4.9350e-7

β0 = 25

0.1282

1.6655e-8

4.89

2.3298e-9

5.36

6.0851e-7

4.43

β1 = 1/40

0.0658

5.4526e-10

4.93

5.4407e-11

5.42

2.7672e-8

4.46

order

keh kL2

order

3.2150e-5

|eh |H 1

order

1.7572e-3

4.8408e-7

4.9419e-5

9.5604e-8

1.3149e-5

In the third test, we implement the symmetric DDG method directly over curved elements that intersect with the interface. The two jump conditions (1.3)-(1.4) are applied on the curved edge of the curved triangular element and there is no geometric error introduced for the circular interface. The volume integration is calculated over the curved triangle as discussed in section §2. We obtain optimal convergence orders on the L∞ and L2 errors with P2 , P3 and P4 polynomials approximations, see Table 6. We present the P2 numerical solution and a sample mesh with h = 0.0658 in Figure 6. We also carry out tests with the outside and inside diffusion coefficient ratio, α+ /α− , varying from 10−4 to 103 . Similar convergence orders, as in Table 6, are observed. We highlight that, for all rest examples in this section, the symmetric DDG method is implemented on curved triangular elements, which fall around the interface. Example 5.3. Variable diffusion coefficients and curved interface In this example, we consider a problem with piecewise defined variable diffusion coefficient α in (1.1) associated with non-zero solution jump φ 6= 0 and non-zero flux jump ψ 6= 0 across the interface. Geometries of the computational domain Ω and the circular interface Γ are the same as Example 5.2. The diffusion coefficient α is piecewise defined as,   α− = (xy + 2)/5, r ≤ 21 α(x, y) =  α+ = (x2 − y 2 + 3)/7, r > 1 . 2

The exact solution is available with,   2, r ≤ 12 , uex (x, y) =  sin(x + y), r > 1 . 2 We implement symmetric DDG method on this variable coefficient problem with P2 and P3 polynomials. Numerical solution errors and convergence orders are presented in Table 7. We obtain optimal (k + 1)th orders of convergence on the L2 and L∞ errors and kth order on the energy norm errors. Example 5.4. Three domains with two interfaces In this example, we study a benchmark problem from [26] or [15] with two interfaces and mixed jump conditions. The computational domain is set up as Ω = (−1, 1) × (−1, 1), which is divided into three sub domains of Ω− , Ω2 and Ω3 , see Figure 7. We name Ω+ = Ω2 ∪ Ω3 , since the diffusion coefficient α+ is continuously defined over Ω+ . The first interface Γ1 is defined and connected by two pieces of a straight line, y = 2x for x + y > 0, and a quadratic curve, y = 2x + x2 for x + y ≤ 0. Interface Γ1 divides the domain Ω into Ω− and Ω+ to be correspondingly the left and the right part of Γ1 . The second interface Γ2 is a straight line, x + y = 0 with x ∈ (0, 1), which further divides Ω+ into two sub domains of Ω2 and Ω3 . Let

20

Table 7: Example 5.3 with piecewise defined diffusion coefficients k, β0 , β1

h

keh kL∞

2

0.1249

2.2116e-5

β0 = 18

0.0638

2.9039e-6

2.93

5.8294e-7

3.00

1.2864e-4

1.98

β1 = 1/40

0.0322

3.7635e-7

2.95

7.2747e-8

3.00

3.2325e-5

2.00

3

0.2397

1.7451e-6

β0 = 9

0.1249

1.4149e-7

3.62

1.8467e-8

3.93

2.4232e-6

2.86

β1 = 1/40

0.0638

9.6208e-9

3.88

1.1547e-9

4.00

3.1211e-7

2.96

order

keh kL2

order

4.6639e-6

|eh |H 1

order

5.0726e-4

2.8090e-7

1.7579e-5

Figure 7: Illustration of three subdomains Ω = Ω− ∪ Ω2 ∪ Ω3 and two interfaces Γ1 and Γ2 of Example 5.4. Ω2 = {(x, y)| x + y < 0, (x, y) ∈ Ω+ } and Ω3 = {(x, y)| x + y > 0, (x, y) ∈ Ω+ } (see Figure 7). The diffusion coefficient α is two piecewise defined as,   α− = (xy + 2)/5, (x, y) ∈ Ω− , α(x, y) =  α+ = (x2 − y 2 + 3)/7, (x, y) ∈ Ω+ . The exact solution uex (x, y) is three piecewise defined with,    2, (x, y) ∈ Ω− ,   u(x, y) = sin(x + y), (x, y) ∈ Ω2 ,     x + y, (x, y) ∈ Ω3 . This problem is equipped with two interfaces and mixed interface conditions. Across interface Γ1 we have non-zero solution jump φ 6= 0 and non-zero flux jump ψ 6= 0 conditions. Over Γ2 , we have zero solution jump φ = 0 and zero flux jump ψ = 0. The problem is solved with symmetric DDG method and we obtain optimal convergence orders on the L∞ and L2 errors, see Table 8. We also output the numerical solution with P2 approximations and a sample mesh of h = 0.0384 in Figure 8. Example 5.5. High-contrast diffusion coefficient 21

Figure 8: Example 5.4 with mesh h = 0.0384 (left) and P2 symmetric DDG solution (right).

Table 8: Example 5.4 with three domains and two interfaces. k, β0 , β1

h

keh kL∞

2

0.0768

2.2811e-5

β0 = 32

0.0384

3.2577e-6

2.81

5.5411e-7

3.00

6.7251e-5

1.96

β1 = 1/40

0.0192

3.4826e-7

3.23

6.9406e-8

3.00

1.6581e-5

2.02

3

0.1536

1.9777e-6

β0 = 60

0.0768

1.4330e-7

3.79

1.8859e-8

3.99

1.7486e-6

3.08

β1 = 1/40

0.0384

1.0211e-8

3.81

1.1864e-9

4.00

2.2017e-7

2.99

order

keh kL2

order

4.4172e-6

order

2.6094e-4

2.9731e-7

22

|eh |H 1

1.4738e-5

Table 9: L2 norm errors with (α+ , α− ) = (1, α ˆ ), β0 = 18, β1 = 1/40, k = 2. h

α ˆ = 10

α ˆ = 100

α ˆ = 1000

α ˆ = 10000

α ˆ = 100000

0.1282

2.7416e-5

2.7428e-5

2.7566e-5

2.7847e-5

2.7920e-5

0.0658

3.4532e-6

3.4504e-6

3.4513e-6

3.4625e-6

3.4721e-6

0.0333

4.3385e-7

4.3335e-7

4.3335e-7

4.3353e-7

4.3421e-7

rate

2.99

2.99

2.99

3.00

3.00

Table 10: L2 norm errors with (α+ , α− ) = (α ˆ , 1), β0 = 20, β1 = 1/40, k = 2. h

α ˆ = 10

α ˆ = 100

α ˆ = 1000

α ˆ = 10000

α ˆ = 100000

0.2420

1.0274e-4

1.0127e-4

1.0418e-4

1.0925e-4

1.1009e-4

0.1282

1.2923e-5

1.2896e-5

1.2747e-5

1.3471e-5

1.3844e-5

0.0658

1.5710e-6

1.5633e-6

1.5574e-6

1.5773e-6

1.6539e-6

0.0333

1.9407e-7

1.9207e-7

1.9207e-7

1.9180e-7

1.9607e-7

rate

3.02

3.02

3.02

3.04

3.08

This example is from [9] that involves a circular interface problem with high contrast diffusion coefficient across the interface. The computational domain is a square Ω =p(−1, 1) × (−1, 1). The interface Γ is a circle centered at the origin with radius R = π/6.28. Denote r = x2 + y 2 for any (x, y) ∈ Ω. The diffusion coefficient α is two piecewise constants with,   α− , r < R, α(r, θ) =  α+ , r ≥ R. The exact solution in polar coordinates is available with,   r3 , α− uex (r, θ) =  r3 + 1 − α+

α−

r < R, 1 α+




R3 ,

(5.1)

r ≥ R.

We have the right hand side function f = −9r. This problem is associated with zero solution jump φ = 0 but non-zero flux jump ψ 6= 0 interface conditions. We carry out a sequence of numerical tests with diffusion − −5 coefficient ratio α to 105 . We demonstrate the capability of symmetric DDG α+ taken from the range of 10 method that obtains uniform optimal convergence orders, which is independent of the diffusion coefficient setup. The tests are classified into two groups. The first group refers to the inside to outside diffusion coefficient − ˆ gradually increasing from 10 to 105 . Errors in L2 norm with P2 approximations are listed in ratio α α+ = α − 1 Table 9. We observe uniform 3rd order of convergence. The second group of tests involves the ratio α α+ = α ˆ −1 −5 2 decreasing from 10 to 10 . Errors in L norm with P2 approximations are listed in Table 10. Again we obtain uniform 3rd order of accuracy. We also list symmetric DDG method solution errors and orders with P2 and P3 approximations in Table 11 and Table 12, corresponding to the two extreme cases of (α+ , α− ) = (1, 105 ) and (α+ , α− ) = (105 , 1), respectively. Uniform optimal convergence orders are observed. Finally, we present symmetric DDG solutions with P2 approximations in Figure 9, which match well with those in literature. Example 5.6. Complex non convex interface In this example, we consider an elliptic interface problem [26] with a flower pedal shape interface that consists both concave and convex curved segments. The computational domain is Ω = (−1, 1) × (−1, 1). The interface Γ is parameterized with polar coordinates (r, θ) as, r(θ) =

1 sin(5θ) + , 2 7 23

θ ∈ [0, 2π].

Table 11: Errors and orders of Example 5.5 with (α+ , α− ) = (1, 105 ). k, β0 , β1

h

keh kL∞

2

0.2420

4.2419e-4

β0 = 18

0.1282

5.1424e-5

3.04

2.2016e-5

3.05

2.8200e-3

2.03

β1 = 1/40

0.0658

6.5144e-6

2.98

2.7211e-6

3.02

6.9996e-4

2.01

0.0333

8.2405e-7

2.98

3.4078e-7

3.00

1.7454e-4

2.00

3

0.2420

1.1449e-5

β0 = 20

0.1282

7.8297e-7

3.87

4.3745e-7

2.78

2.7233e-5

3.03

β1 = 1/20

0.0658

7.1091e-8

3.46

1.9349e-8

4.50

3.3480e-6

3.02

order

keh kL2

order

1.8217e-4

|eh |H 1

order

1.1548e-2

2.9999e-6

2.2201e-4

Table 12: Errors and orders of Example 5.5 with (α+ , α− ) = (105 , 1). k, β0 , β1

h

keh kL∞

2

0.2420

3.6107e-4

β0 = 18

0.1282

6.7466e-5

2.42

1.3465e-5

3.00

1.5212e-3

1.97

β1 = 1/40

0.0658

9.8349e-6

2.78

1.5987e-6

3.07

3.7583e-4

2.01

0.0333

1.0319e-6

3.25

1.8870e-7

3.08

9.2599e-5

2.02

3

0.2420

1.7903e-5

β0 = 20

0.1282

3.5908e-6

2.32

3.6934e-7

3.27

4.3243e-5

2.80

β1 = 1/20

0.0658

4.8378e-7

2.89

2.0474e-8

4.17

5.9176e-6

2.87

order

keh kL2

order

1.0773e-4

|eh |H 1

order

5.9615e-3

3.5506e-6

3.0210e-4

Figure 9: Example 5.5 symmetric DDG solutions for the case of (α+ , α− ) = (1, 105 ) (left) and (α+ , α− ) = (105 , 1) (right), P2 approximations with mesh h = 0.0658.

24

Table 13: Errors and orders of Example 5.6 with non convex interface k, β0 , β1

h

keh kL∞

2

0.1525

9.1300e-5

β0 = 12

0.0630

1.4280e-5

2.68

2.5402e-6

2.95

4.6804e-4

1.94

β1 = 1/25

0.0324

1.9295e-6

2.88

3.1657e-7

3.00

1.1832e-4

1.98

3

0.2470

3.7443e-5

β0 = 18

0.1235

5.5438e-6

2.76

4.4100e-7

3.86

7.8116e-5

2.95

β1 = 1/50

0.0630

5.7268e-7

3.28

2.5833e-8

4.09

8.9921e-6

3.12

order

keh kL2

order

1.9663e-5

|eh |H 1

order

1.7940e-3

6.3834e-6

6.0333e-4

Figure 10: Example 5.6 with mesh h = 0.0630 (left) and P2 symmetric DDG solution(right). The diffusion coefficient α is chosen to be two piecewise constants with α− = 1 and α+ = 10 for the inside and outside of the interface Γ. The analytical solution is given as,   er 2 , inside Γ, uex (r, θ) =  0.1r4 − 0.01 ln(2r), outside Γ. Across the interface Γ, we have non-zero solution jump φ 6= 0 and non-zero flux jump ψ 6= 0 conditions. We perform P2 and P3 symmetric DDG methods tests and obtain optimal convergence orders on L∞ and L2 errors, see Table 13. We output P2 symmetric DDG approximation with mesh h = 0.0324 and present it in Figure 10. To further test the capability of symmetric DDG method, we carry out accuracy tests with the outside and inside diffusion coefficient ratio, α+ /α− , varying from 10−3 to 103 . Similar convergence orders, as in Table 13, are observed.

6

Appendix

We have also considered another group of numerical fluxes in the symmetric DDG scheme formulation. We test out the numerical Example 5.4 with three domains and two interfaces with following numerical fluxes and obtain optimal accuracy orders on the L2 norm, L∞ norm and energy norm errors as Table 8. Now we just list the set of numerical fluxes we tested. Over an interior edge e = ∂E ∩ ∂E 0 ∈ EhI , the numerical flux is 25

defined as,    αu dn    αv gn      JuKφ

JuK + {{αun }} + β1 he J∇ · (α∇u)K, := β0 max{α, α0 } he 1 := αvn − β1 he ∇ · (α∇u) 2 := JuK.

On a boundary edge e ∈ EhD , the numerical flux is defined as,  u0 − u  αu dn := β0 α + αun + β1 he (f − ∇ · (α∇u)) ,   he  αv gn := αvn − β1 he ∇ · (α∇v),     JuKφ := u0 − u.

On an interface edge e = Γ ∩ ∂E + with E + ⊂ Ω+ , the numerical flux is defined as,  φ ψ  + − JuK  αu d := β max{α , α } + {{αun+ }} + + β1 he (J∇ · (α∇u)K + Jf K) ,  n 0   he 2 1 αvn − β1 he ∇ · (α∇v), αv gn :=   2    JuKφ := JuK − φ. On an interface edge e = Γ ∩ ∂E − with E − ⊂ Ω− , the numerical flux is defined as,  JuKφ ψ   αu dn := β0 max{α+ , α− } + {{αun− }} + + β1 he (J∇ · (α∇u)K + Jf K) ,    he 2 1 αv gn := αvn − β1 he ∇ · (α∇v),   2    JuKφ := JuK + φ.

(6.1)

(6.2)

(6.3)

(6.4)

Acknowledgements. Huang’s work is supported by Natural Science Foundation of Zhejiang Province grant No.LY14A010002, and is subsidized by the National Natural Science Foundation of China under grant No. 11771398, No.11471195 and No. 61673352.

References [1] S. Adjerid, M. Ben-Romdhane, and T. Lin. Higher degree immersed finite element methods for secondorder elliptic interface problems. International Journal of Numerical Analysis and Modeling, 11(3):541– 566, 2014. [2] D. N. Arnold. An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal., 19(4):742–760, 1982. [3] C. B¨orgers. A triangulation algorithm for fast elliptic solvers based on domain imbedding. SIAM J. Numer. Anal., 27(5):1187–1196, 1990. [4] E. Burman and P. Hansbo. Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method. Comput. Methods Appl. Mech. Engrg., 199(41-44):2680–2686, 2010. [5] R. Bustinza and G. N. Gatica. A local discontinuous Galerkin method for nonlinear diffusion problems with mixed boundary conditions. SIAM J. SCI. COMPUT., 26:152–177, 2004. [6] R. Bustinza, G. N. Gatica, and F. J. Sayas. On the coupling of local discontinuous Galerkin and boundary element methods for non-linear exterior transmission problems. IMA Journal of Numerical Analysis, 28(2):225–244, 4 2008. [7] Z. Chen, H. Huang, and J. Yan. Third order maximum-principle-satisfying direct discontinuous Galerkin methods for time dependent convection diffusion equations on unstructured triangular meshes. J. Comput. Phys., 308:198–217, 2016. 26

[8] Z. M. Chen and J. Zou. Finite element methods and their convergence for elliptic and parabolic interface problems. Numer. Math., 79:175–202, 1998. [9] C. C. Chu, I. G. Graham, and T. Y. Hou. A new multiscale finite element method for high-contrast elliptic interface problems. Mathematics of Computation, 79(272):1915–1955, 10 2010. [10] B. Cockburn and C.-W. Shu. The local discontinuous Galerkin method for time-dependent convectiondiffusion systems. SIAM J. Numer. Anal., 35(6):2440–2463 (electronic), 1998. [11] D. A. Dunavant-1985. Gaussian integration points and corresponding weights for isosceles right triangular quadrature domain. International Journal for Numerical Methods for Engineering, 21:1129–1148, 1985. [12] Y. Gong, B. Li, and Z. L. Li. Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions. SIAM J. Numer. Anal., 46(1):472–495, 2008. [13] J. Guzm´an, M. A. S´ anchez, and M. Sarkis. Higher-order finite element methods for elliptic problems with interfaces. ESAIM: Mathematical Modelling and Numerical Analysis, 50:1561–1583, 2016. [14] X. M. He, T. Lin, and Y. P. Lin. Immersed finite element methods for elliptic interface problems with nonhomogeneous jump conditions. International Journal of Numerical Analysis and Modeling, 8(2):284–301, 2011. [15] S. M. Hou, W. Wang, and L. Q. Wang. Numerical method for solving matrix coefficient elliptic equation with sharp-edged interfaces. Journal of Computational Physics, 229:7162–7179, 2010. [16] H. Huang, Z. Chen, J. Li, and J. Yan. Direct discontinuous Galerkin method and its variations for second order elliptic equations. J. Sci. Comput., 70(2):744–765, 2017. [17] H. Y. Huang. The direct coupling of local discontinuous Galerkin and natural boundary element method for nonlinear interface problem in R3 . Applied Mathematics and Computation, 262:232–248, 2015. [18] P. Q. Huang, H. J. Wu, and Xiao Y. M. An unfitted interface penalty finite element method for elliptic interface problems. Comput. Methods Appl. Mech. Engrg., 323:439–460, 2017. [19] L. Lee and R. J. Leveque. An immersed interface method for incompressible Navier-Stokes equations. SIAM J. Sci. Comput., 25(3):832–856, 2003. [20] R. J. LeVeQue and Z. Li. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J. Numer. Anal., 31:1019–1044, 1994. [21] Z. Li. The immersed interface method using a finite element formulation. Appl. Numer. Math., 27:253–267, 1998. [22] Z. Li, T. Lin, Y. Lin, and R. C. Rogers. An immersed finite element space and its approximation capability. Numer. Methods Partial Differential Equations, 20(3):338–367, 2004. [23] T. Lin, Y. P. Lin, and X. Zhang. Partially penalized immersed finite element methods for elliptic interface problems. SIAM J. NUMER. ANAL., 53(2):1121–1144, 2015. [24] H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) methods for diffusion problems. SIAM J. Numer. Anal., 47(1):475–698, 2009. [25] H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) method for diffusion with interface corrections. Commun. Comput. Phys., 8(3):541–564, 2010. [26] L. Mu, J. P. Wang, G. W. Wei, X. Ye, and S. Zhao. Weak Galerkin methods for second order elliptic interface problems. Journal of Computational Physics, 250:106–125, 2013. [27] C. S. Peskin. Numerical analysis of blood flow in the heart. J. Comput. Phys., 25:220–252, 1977. [28] C. S. Peskin. The immersed boundary method. Acta Numer., 11:479–517, 2002. [29] B. Rivi`ere. Discontinuous Galerkin methods for solving elliptic and parabolic equations, volume 35 of Frontiers in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008. Theory and implementation. 27

[30] C. Vidden and J. Yan. A new direct discontinuous Galerkin method with symmetric structure for nonlinear diffusion equations. J. Comput. Math., 31(6):638–662, 2013. [31] H. Y. Wei, L. Chen, Huang Y. Q., and B. Zheng. Adaptive mesh refinement and superconvergence for two-dimensional interface problems. SIAM J. SCI. COMPUT., 36(4):A1478–A1499, 2014. [32] M. F. Wheeler. An elliptic collocation-finite element method with interior penalties. SIAM J. Numer. Anal., 15:152–161, 1978. [33] M. Zhang and J. Yan. Fourier type super convergence study on DDGIC and symmetric DDG methods. J. Sci. Comput., 73(2-3):1276–1289, 2017. [34] X. Zhang. A curved boundary treatment for discontinuous Galerkin schemes solving time dependent problems. J. Comput. Phys., 308:153–170, 2016. [35] Y. C. Zhou and G. W. Wei. On the fictitious-domain and interpolation formulations of the matched interface and boundary (MIB) method. Journal of Computational Physics, 219:228–246, 2006.

28

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

Jue Yan