# EnrichedP1-Conforming Methods for Elliptic Interface Problems with Implicit Jump Conditions

Authors: Gwanghyun Jo; Do Y. Kwak
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

We develop a numerical method for elliptic interface problems with implicit jumps. To handle the discontinuity, we enrich usualP1-conforming finite element space by adding extra degrees of freedom on one side of the interface. Next, we define a new bilinear form, which incorporates the implicit jump conditions. We show that the bilinear form is coercive and bounded if the penalty term is sufficiently large. We prove the optimal error estimates in both energy-like norm and L2-norm. We provide numerical experiments. We observe that our scheme converges with optimal rates, which coincides with our error analysis.


---

## Body

## 1. Introduction

Interface problems arise in various disciplines including mechanical, material, and medical image and petroleum engineering [1–9]. There are several difficulties to solve for the governing equation of such problems.Firstly, partial differential equations may have different coefficients along the interface due to the change of material properties. When the geometry of interface is complex, one needs to generate grids that align with the interfaces. Once a fitted grid is generated, one uses finite element method (FEM) or finite volume method (FVM) based on this grid. Secondly, the problem may have nonhomogeneous jump conditions along the interface. When the jumps along the interface are known explicitly, (say[u]=g1,∂u/∂n=g2, with known g1 and g2), these jumps may be handled effectively by discontinuous Galerkin (DG) [10, 11] by incorporating jumps into the bilinear form with proper penalty terms. For example, an effective DG scheme was developed to describe discontinuous phenomena arising from porous media with discontinuous capillary pressure [12]. The interface problems with known jumps can be solved with immersed interface methods [13–15] or discontinuous bubble-immersed finite element methods [16].However, when the jumps are implicit along the interface problems, numerically solving the governing equations becomes more challenging. Let us consider some problems with interface conditions, where the jumps of primary variables are related to the normal fluxes. Firstly, these problems arise in the medical imaging of cancer cells using MREIT [3, 4] or electrochemotherapy [17], where the jumps of an electric voltage across the cell membrane appear. Next, an elastic body has spring-type jumps that are related to stress [18, 19]. The heat in the material interface may have implicit jump conditions along the interface [20, 21]. Also, a generalized jump condition for Laplace equation or Helmholtz equation has been considered in [22–24].The first attempt to solve the elliptic interface problems having implicit jump conditions seems to be introduced in [25], where the iterative method was used. Recently, some XFEM-based nonfitted methods were developed in [26, 27] for the elliptic problems and elasticity problems, respectively, where the extra degrees of freedom are introduced on elements cut by the interface. On the other hand, an immersed finite element type method was developed in [28].In this work, we introduce a new numerical method to solve elliptic interface problems, where the jumps are related to the normal fluxes and some known functions. A main idea of our work is to include the jump conditions implicitly on the bilinear form so that the numerical solutions for the weak problems satisfy the implicit jump conditions. We enrich the usualP1 FEM space near the interface. We show that our bilinear form is coercive and bounded and prove the optimal error estimates. In numerical section, we provide several numerical examples supporting our analysis.LetΩ be a convex domain in Rn (n=2,3), which is divided into Ω+ and Ω- by a C2 closed interface Γ. The governing equations on Ω are given by(1)-∇·β∇u=f,inΩ≔Ω-∪Ω+,(2)u=g0,on∂Ω,(3)uΓ=α∂u+∂n++g1,onΓ,(4)β∂u∂nΓΓ=g2,onΓ,where f∈L2(Ω) and g0∈H3/2(∂Ω), g1∈H3/2(Γ), g2∈H1/2(Γ), and β is a positive piecewise constant; that is, β=β+ in Ω+ and β=β- in Ω-, where β+ and β- are some positive constants. Here, ns is the outer unit normal vector to Ωs (s=-,+) and ·Γ is the jump along the interface; that is, [u]Γ=u|Ω--u|Ω+. Also, we define nΓ to be an outer normal vector to Ω-. The jump of normal derivatives of u is defined as (5)β∂u∂nΓΓ=β-∂uΩ-∂nΓ-β+∂uΩ+∂nΓ.We assume that α is a positive constant.We introduce some notations. LetO be any domain and let Hm(O), m=1,2, be a usual Sobolev space with norm um,O. We define H01(O) as the set of functions in H1(O) with vanishing trace on ∂O. We define the subspaces of Hm(O), (6)H~mO≔HmΩ-∩O∩HmΩ+∩O,equipped with the (semi)norms:(7)uH~mO≔uHmΩ-∩O+uHmΩ+∩O,uH~mO≔uHmΩ-∩O+uHmΩ+∩O.Finally, we define subspace of H~1(Ω): (8)H~01Ω≔u:u∈H~1O∣u∂Ω=0.We state a theorem regarding the existence and regularity of the problem [29, 30].Theorem 1.
Problem (1)–(4) has a unique solution u∈H~2(Ω) such that, for some constant C>0,(9)uH~2Ω≤CfL2Ω+g0H3/2∂Ω+g1H3/2Γ+g2H1/2Γ.The rest of the paper is organized as follows. In Section2, we derive the variational forms for the problems with implicit jump conditions. We introduce new numerical methods in Section 3 and in Section 4 we prove the error estimates. In Section 5, we give numerical results that support our analysis. The conclusion follows in Section 6.


## 2. Variational Form

In this section, we derive a variational formulation of the model problem. Without loss of generality, we may assume thatg0=0. First, we multiply v∈H~01(Ω) to (1) and apply integration by parts on each subdomain to get (10)∫Ωsβ∇u·∇vdx-∫∂Ωsβs∂us∂nsvds=∫Ωsfvdx,s=+,-.By summation, we have (11)∑s=+,-∫Ωsβ∇u·∇vdx-∑s=+,-∫∂Ωsβs∂us∂nsvds=∑s=+,-∫Ωsfvdx.Using the jump conditions (3) and (4), we see the second terms become(12)-∑s=+,-∫∂Ωsβs∂us∂nsvds=-∫∂Ω-β-∂u-∂nΓv-ds+∫∂Ω+β+∂u+∂nΓv+ds,=-∫Γβ+∂u+∂nΓvΓds-∫Γg2v-ds,by4=-∫Γβ+αg1-uΓvΓds-∫Γg2v-ds,by3=β+α∫ΓuΓvΓds-β+α∫Γg1vΓds-∫Γg2v-ds.We define a bilinear form and a functional on H1~(Ω):(13)au,v≔∫Ω-β∇u·∇vdx+∫Ω+β∇u·∇vdx+β+αuΓ,vΓΓ,(14)F~v≔f,v+β+αg1,vΓΓ+g2,v-Γ,where (·,·) denotes L2(Ω) inner product on Ω and ·,·Γ denotes the L2(Γ) inner product. By (12) and (13), we have the weak problem: find u∈H~2(Ω)∩H~01(Ω) satisfying(15)au,v=F~v,for all v∈H~01(Ω).Now let us show that the weak problem (15) is equivalent to (1)–(4). Suppose that u∈H~2(Ω)∩H~01(Ω) satisfies (15). First, let v be any function H01(Ω-) (or v∈H01(Ω+)). Then, we have (16)∫Ω-β∇u·∇vdx=∫Ω-fvdx,or,∫Ω+β∇u·∇vdx=∫Ω+fvdx.By integration by parts, we see that u satisfies (17)-∇·β∇u=f,onΩ-∪Ω+.Now, assume thatv∈H01(Ω) in (15). By Green’s theorem and the fact that vΓ=0, the left side of (15) becomes (18)∑s=+,-∫∂Ωsβs∂us∂nsvds-∫Ωs∇·β∇uvdx=∫Γβ∂u∂nΓΓvds+∫Ωfvdx.Comparing with the right side of (15), we have (19)β∂u∂nΓΓ=g2.Finally, assume that v∈H~01(Ω) in (15). By integration by parts, the left hand side of (15) becomes (20)∑s=+,-∫∂Ωsβs∂us∂nsvds-∫Ωs∇·β∇uvdx+β+α∫ΓuΓvΓds=∫Γβ+∂u+∂nΓvΓds+∫Γg2v-+β+α∫ΓuΓvΓds+∫Ωfvdx.Comparing with the right side of (15), we see that u satisfies (21)uΓ=α∂u+∂n++g1,onΓ.


## 3. Numerical Methods

In this section, we develop a numerical method for (1)–(4). Our method is obtained by adding extra degrees of freedom to P1-conforming space on one side of the interface. For simplicity, we assume that Ω⊂R2. However, similar constructions are possible for the case of Ω⊂R3 as well.LetTh be a given regular triangulation of Ω fitted with the interface. We let Th+ and Th- be set of elements in Th which belong to Ω+ and Ω-, respectively. We let Sh(Ω) be the usual P1-conforming space; that is, any function in Sh(Ω) is continuous and piecewise linear and is vanishing on the boundary. We use notation Sh(T) for the set of linear functions on T.We letIh be the set of all neighboring elements of interface Γ in Th+; that is, T∈Th belongs to Ih if and only if T∈Th+ and at least one node of T is located on Γ. We let Bh(T) be the space of functions in Sh(T) vanishing on nodes not lying at the interface. For example, suppose that T has three nodes A1, A2, and A3, where A1 and A2 are located on Γ. Then, a function in Bh(T) is linear on T vanishing at A3. In this case, Bh(T) has dimension two. On the other hand, if T have only one node located on Γ, the dimension of Bh(T) is one. A function in Bh(T) is extended to Ω as follows:(22)BhΩ=ϕT=0,onT∈Th-Ih,ϕT∈BhT,onT∈Ih.For example, suppose that there are seven elements aligning with interface (see Figure1). Then function ϕh in Bh(Ω) has a support on grey region. Moreover, ϕh has vanishing values on outside nodes on Ω+. Thus, ϕh has seven degrees of freedom, that is, a,b,c,d,e,f,g,h.Figure 1
The support and the degrees of freedom of a functionϕh in Bh(Ω).We decomposeu as(23)u≔u0+u∗,where u0 belongs to H01(Ω) and u∗ belongs to H~01(Ω). We approximate u0 from Sh(Ω) and approximate u∗ from Bh(Ω). Thus, u=u0+u∗ is approximated in S~h(Ω)≔Sh(Ω)+Bh(Ω).Now, we discretize (1)–(4). Let Eh be edges of elements in Th and let us define subspaces of Eh (see Figure 1).(i)
EhI is the set of edges of Eh whose two endpoints are located on Γ.(ii)
EhN is the set of edges of Eh whose one endpoint is located on Γ and the other is located in the interior Ω-.(iii)
EhO≔Eh-EhI-EhN.We note thatEh=EhI∪EhN∪EhO. For all edges e in Eh, we fix a normal vector ne once and for all. We define jumps and averages across the edges: (24)ue≔uTe1-uTe2,ue≔0.5uTe1+0.5uTe2,where Te1 and Te2 are two neighboring elements of e.We multiplyvh∈S~h(Ω) to (1) and use integration by parts to obtain the following:(25)∑T∈Th∫Tβ∇u·∇vhdx+∑e∈EhDe=∑T∈Thf,vhT,where (26)De≔-∫eβ∇uTe1·nevhTe1-β∇uTe2·nevhTe2ds.Let us classify De into three categories. Firstly, if e∈EhI, then using the similar method used in deriving (12), we have(27)De=β+α∫euevheds-β+α∫eg1vheds-∫eg2vh-ds.Secondly, if e∈EhN, then by using the identity (28)a-b--a+b+=a-+a+b--b+2+a--a+2b-+b+and the fact that (29)β∇u·nee=0,we have (30)De=-∫eβ∇u·neevheds-∫eβ∇u·neevheds,=-∫eβ∇u·neevheds.Finally, if e∈EhO, De vanishes, since both β∇u and vh are continuous across e. Thus, (25) becomes(31)∑T∈Th∫Tβ∇u·∇vhdx-∑e∈EhN∫eβ∇u·neevheds+∑e∈EhIβ+α∫euevheds=∑T∈Thf,vhT+∑e∈EhIβ+α∫eg1vheds+∑e∈EhI∫eg2vh-ds.Now we propose our method based on enrichedP1-conforming space: find uh∈S~h(Ω) satisfying(32)ahϵuh,vh=F~hvh,for all vh∈S~h(Ω), where(33)ahϵu,v≔∑T∈Th∫Tβ∇u·∇vdx-∑e∈EhN∫eβ∇u·neeveds+ϵ∑e∈EhN∫eβ∇v·neeueds+∑e∈EhN∫eσeueveds+∑e∈EhIβ+α∫eueveds,(34)F~hvh=f,vh+∑e∈EhIβ+α∫eg1vheds+∑e∈EhI∫eg2vh-ds.In (33), the parameter σ is positive and the parameter ϵ is 1, 0, or -1, which is motivated by NIPG, IIPG, and SIPG of DG scheme [11].We show that our scheme is consistent.Lemma 2.
Suppose thatu is the solution of (1)–(4). Then, for all vh∈S~h(Ω), the following holds: (35)ahϵu,vh=F~hvh.In other words, we have(36)ahϵu-uh,vh=0,for all vh∈S~h(Ω).Proof.
Since[u]e=0 for all e∈EhN, we have (37)ahϵu,vh≔∑T∈Th∫Tβ∇u·∇vhdx-∑e∈EhN∫eβ∇u·neevheds+∑e∈EhIβ+α∫euevheds,for all vh∈S~h(Ω). By (31), we have the conclusion.


## 4. Error Estimates

We define energy-like norm⦀·⦀h on Hh(Ω)≔H~2(Ω)+S~h(Ω).(38)⦀ϕ⦀2=∑T∈Th∫Tβ∇ϕ·∇ϕdx+∑e∈EhN1e∫eϕe2ds+∑e∈EhI∫eϕe2ds.Let γ0:H1(T)→H1(∂T) be the usual trace operator. Then we have the following theorem [31, 32].Lemma 3 (trace theorem).
There exists a constantC such that(39)γ0vL2∂T≤Ch-1/2vL2T+h1/2∇vL2T,for all v∈H1(T).We define a local interpolation operatorIh:u∈H2(T)→Ihu∈S~h(T) by (40)IhuXi=uXi,i=1,2,3,where Xi’s are nodes of T. The operator Ih is extended to u∈H~2(Ω) by (Ihu)T=Ih(uT) for each element T. We note that if e belongs to EhI, Ihu is discontinuous across e; that is, [Ihu]e≠0. However, Ihu is continuous on each subdomain Ω- and Ω+. Then, by the standard interpolation theory, there exists a constant C such that(41)∑T∈Th+u-Ihu0,T+h∑T∈Th+∇u-Ihu0,T≤Ch2uH~2Ω+,∑T∈Th-u-Ihu0,T+h∑T∈Th-∇u-Ihu0,T≤Ch2uH~2Ω-.From this, we can obtain interpolation estimate in ⦀·⦀ norm.Corollary 4.
There exists a constantCI such that, for all u∈H~2(Ω),(42)⦀u-Ihu⦀≤CIhuH~2Ω.Proof.
Sinceu-uh is continuous across e∈EhN, we have (43)⦀u-Ihu⦀2=∑T∈Thβ1/2∇u-Ihu0,T2+∑e∈EhI∫eu-Ihu2ds≤Ch2uH~2Ω2+∑e∈EhI∫eu-Ihu2ds.It suffices to show that(44)∑e∈EhI∫eu-Ihu2ds≤Ch2uH~2Ω,for some constant C>0. Let ϕ=u-Ihu. By the trace inequality (41), we have (45)∑e∈EhI∫eϕ2ds≤∑e∈EhI∑i=1,2ϕTei0,e2≤C∑e∈EhI∑i=1,2h-1ϕ0,Tei2+h1∇ϕ0,Tei2≤∑e∈EhI∑i=1,2Ch3uH2Tei2≤Ch2∑T∈ThuH2T2.We have the following coercivity property.Theorem 5.
If we chooseσ so that(46)σ>31-ϵ2C122,then the following holds:(47)ahϵvh,vh≥C1⦀vh⦀2,for all vh∈S~h(Ω), where (48)C1=min⁡12,σ,β+α.Proof.
Ifϵ=1, then by definition of ah(·,·) and ⦀·⦀ norm, we have (49)ahϵvh,vh≥min⁡1,σ,β+αvhh2.Now, suppose that ϵ=-1 or 0. We have (50)ahϵvh,vh=∑T∈Thβ1/2∇vh0,T2-1-ϵ∑e∈EhN∫eβ∇vh·neevheds+∑e∈EhNσevhe0,e2+∑e∈EhIβ+αvhe0,e2≔A1+A2+A3+A4.We bound A2 by Cauchy-Schwarz inequality, (39), and Young’s inequality. (51)1-ϵ∑e∈EhN∫eβ∇vh·neevheds≤1-ϵ∑e∈EhNβ1/2vhe0,evhe0,e≤∑e∈EhN1-ϵC1h-1/212β1/2∇vh0,Te1+12β1/2∇vh0,Te2vhe0,e≤1-ϵC1h-1/2∑e∈EhNβ1/2∇vh0,Te121/2∑e∈EhNvhe0,e21/2≤1-ϵC1h-1/23∑T∈Thβ1/2∇vh0,T21/2∑e∈EhNvhe0,e21/2≤ϵ12∑T∈Thβ1/2∇vh0,T2+31-ϵ2C122ϵ1∑e∈EhN1evhe0,e2,for all ϵ1>0. Thus, we have (52)ahϵvh,vh≥1-ϵ12∑T∈Thβ1/2∇vh0,T2+σ-31-ϵ2C122ϵ1∑e∈EhN1evhe0,e2+β+α∑e∈EhIvhe0,e2.Then (47) is obtained by taking ϵ=1.By a similar technique, we can show thatah is bounded.Theorem 6.
There exists a constantCB>0 such that following holds:(53)ahϵuh,vh≤CB⦀uh⦀·⦀vh⦀h,for all uh,vh∈S~h(Ω).Finally, we prove the error estimate in the energy-like norm.Theorem 7.
Assume thatσ>0 satisfies (46). There exists a constant C>0 such that(54)⦀u-uh⦀≤ChfL2Ω+g1H3/2Γ+g2H1/2Γ.Proof.
By (36), we have(55)ahϵu-uh,u-uh=ahϵu-uh,u-ϕh,for all ϕh∈S~h(Ω). By (55), (47), and (53) we have (56)C1⦀u-uh⦀2≤ahϵu-uh,u-ϕh≤CB⦀u-uh⦀·⦀u-ϕh⦀.If we take ϕh=Ihu, then, by (42) and (9), we have (57)⦀u-uh⦀≤CBC1CIhuH~2Ω≤CBC1CIhfL2Ω+g1H3/2Γ+g2H1/2Γ.Next, we proveL2 estimates using duality argument when ϵ=-1.Theorem 8.
Ifσ>0 satisfies (46) and ϵ=-1, there exists a constant C>0 such that (58)u-uhL2Ω≤Ch2fL2Ω+g1H3/2Γ+g2H1/2Γ.Proof.
We define an auxiliary problem. Letψ∈H~2(Ω) be solution of(59)-∇·β∇ψ=eh,onΩ+∪Ω-,(60)ψΓ=-α∂ψ+∂n+,onΓ,(61)β∂ψ∂nΓ=0,onΓ,(62)ψ=0,on∂Ω,where eh≔u-uh∈L2(Ω). We multiply eh to (59) and we use integration by parts to have(63)ehL2Ω2=-∑T∈Th∫T∇·β∇ψehdx=∑T∈Th∫Tβ∇ψ·∇ehdx-∑e∈Eh∫eβ∇ψTe1·neehTe1-β∇ψTe2·neehTe2ds.We use similar techniques in the classification of De of (25) to derive(64)-∑e∈Eh∫eβ∇ψTe1·neehTe1-β∇ψTe2·neehTe2ds=-∑e∈EhN∫eβ∇ψ·neeeheds+∑e∈EhIβ+α∫eψeeheds.Combining with (63) and (64) and the fact that ψ is continuous on each subdomain Ω- and Ω+, we have (65)ehL2Ω2=ahϵψ,eh.By definition of eh and the fact that ahϵ(·,·) is symmetric and by (36), (54), (42), and (9), we have (66)u-uhL2Ω2=ahϵu-uh,ψ=ahϵu-uh,ψ-Ihψ≤ChuH~2ΩChψH~2Ω≤Ch2fL2Ω+g1H3/2Γ+g2H1/2Γu-uhL2Ω.Thus, we have the conclusion.


## 5. Numerical Results

In this section, we provide some numerical experiments of elliptic interface problems with implicit jump conditions. We consider circle- and ellipse-type interface shapes.We let the domainΩ=[-1,1]2 and we let Th be a triangulation of Ω by regular triangles, which aligns with interfaces. We set ϵ=-1 in bilinear form (33). We take σ in (33) as a multiple of β. We report the number of elements, degrees of freedom, L2-errors, and H1-errors for hk=2-k, k=1,2,…, in Tables 1 and 2. For both examples, we observe optimal error convergence, which supports our analysis in Section 4.Table 1
L 2 and H1 errors of Example 1.
Elements DoF u - u h L 2 ( Ω ) Order u - u h 1 , h Order 70 80 3.934 × 1 0 - 1  4.629 × 1 0 - 1 264 221 1.019 × 1 0 - 1 1.949 1.793 × 1 0 - 1 1.368 1062 709 2.568 × 1 0 - 2 1.988 8.009 × 1 0 - 2 1.163 4171 2439 6.432 × 1 0 - 3 1.998 3.940 × 1 0 - 2 1.023 16557 8986 1.610 × 1 0 - 3 1.998 1.967 × 1 0 - 2 1.002 65521 34174 4.023 × 1 0 - 4 2.001 9.832 × 1 0 - 3 1.000 260961 133311 1.006 × 1 0 - 4 2.000 4.911 × 1 0 - 3 1.002Table 2
L 2 and H1 errors of Example 2.
Elements DoF u-uhL2(Ω) Order u-uh1,h Order 110 118 5.964 × 1 0 - 3  1.008 × 1 0 - 1 382 313 2.720 × 1 0 - 3 1.133 7.348 × 1 0 - 2 0.456 1355 920 6.077 × 1 0 - 4 2.162 3.589 × 1 0 - 2 1.034 5005 2983 1.714 × 1 0 - 4 1.826 1.924 × 1 0 - 2 0.900 18993 10461 4.321 × 1 0 - 5 1.988 9.900 × 1 0 - 3 0.959 74225 39037 1.090 × 1 0 - 5 1.987 4.989 × 1 0 - 3 0.989 292404 150041 2.769 × 1 0 - 6 1.976 2.519 × 1 0 - 3 0.986Example 1 (circular interface).
The interface is given byΓ={(x,y):x2+y2=r02} and Ω- and Ω+ are inside and outside of Γ, respectively. The coefficients are (β-,β+)=(10,1), α=1, and r0=0.5. The exact solution is (67)u=r22β--αr0+1β+-1β-r022+1.5-r2,onΩ-,r22β+,onΩ+,where r=x2+y2. We remark that u satisfies jump conditions (3) and (4), where g1 and g2 are given as (68)g1=1.5-r02,g2=-20r0.Table 1 shows the number of elements and L2 and piecewise H1 errors. Figure 2 shows the numerical solution. We observe that our scheme has optimal convergence in L2 and piecewise H1-norms.Figure 2
Numerical solution of Example1.Example 2 (elliptical shape interface).
The interface is given byΓ={(x,y):x2+y2/2=r02} and Ω- and Ω+ are inside and outside of Γ, respectively. The coefficients are (β-,β+)=(1,10) and r0=0.4. In the previous example, α was constant. However, in this example, we set α as a function of (x,y): (69)α=1x2+y2/41/2.The exact solution is(70)u=q2β--12β++1β+-1β-r022+xy,onΩ-,q2β+onΩ+,where q=x2+y2/4. We remark that g1 and g2 in (3) and (4) are given as (71)g1=xy,g2=3xy4x2+y2.Table 2 shows the errors and Figure 3 shows the numerical solution. Again, we observe the optimal convergence.Figure 3
Numerical solution of Example2.


## 6. Conclusion

In this work, we introduce a numerical method for elliptic interface problems, where the jumps are related to the normal fluxes. We enrich usualP1 space by extra degrees of freedom on one side of the interface. We define bilinear form that includes the jump conditions implicitly. We prove that the bilinear form is coercive and bounded. Using Cea’s Lemma, we prove the error estimates in energy-like norm. Next, we prove L2 error estimate using the duality arguments. We provide numerical experiments that support our analysis.



---

*Source: 9891281-2018-05-09.xml*