aa r X i v : . [ m a t h . NA ] J u l A MODIFIED P - IMMERSED FINITE ELEMENT METHOD DO Y. KWAK AND JUHO LEE
Abstract.
In recent years, the immersed finite element methods (IFEM) in-troduced in [20], [21] to solve elliptic problems having an interface in thedomain due to the discontinuity of coefficients are getting more attentions ofresearchers because of their simplicity and efficiency. Unlike the conventionalfinite element methods, the IFEM allows the interface to cut through the inte-rior of the element, yet after the basis functions are altered so that they satisfythe flux jump conditions, it seems to show a reasonable order of convergence.In this paper, we propose an improved version of the P based IFEM byadding the line integral of flux terms on each element. This technique resemblesthe discontinuous Galerkin (DG) method, however, our method has much lessdegrees of freedom than the DG methods since we use the same number ofunknowns as the conventional P finite element method.We prove H and L error estimates which are optimal both in order andregularity. Numerical experiments were carried out for several examples, whichshow the robustness of our scheme. Introduction
In recent years, there have been some developments of immersed finite elementmethods for elliptic problems having an interface. These methods use meshes whichdo not necessarily align with the discontinuities of the coefficients [20], [21], thusviolate a basic principle of triangulations in the conventional finite element methods[4], [11]. However, when the basis functions are modified so that they satisfy theinterface conditions, they seem to work well [10], [20], [21]. These methods wereextended to the case of Crouzeix-Raviart P nonconforming finite element method[12] by Kwak et al. [18], and to the problems with nonzero jumps in [7]. Somerelated works on interface problems can be found in [5], [16], [17], [19], [22], [23],[26].On the other hand, the discontinuous Galerkin methods (DG) where one usescompletely discontinuous basis functions were developed and have been studiedextensively, see [1], [2], [13], [24] and references therein. The DG methods workquite well for problems with discontinuous coefficient in the sense that they capturethe sharp changes of the solutions well, yet they require large number of unknownsand the meshes have to be aligned with the discontinuity.The purpose of this paper is to combine the advantages of the two methods. Weuse a DG type idea of adding the consistency terms to the IFEM, thus proposing amodified version of IFEM based on the P - Lagrange basis functions on triangulargrids. In spirit, it resembles [15] in the sense that the standard linear basis func-tions are used for noninterface elements and line integrals are added, but in ourmethod the line integrals along the edges, not along the interface, are added. Fur-thermore, our method incorporate the flux jump conditions to the basis functionshence requires no extra unknowns along the interface as in [15]. Mathematics Subject Classification. primary 65N30, secondary 74S05, 76S05 .
Key words and phrases. modified P -immersed finite element, flux jump, discontinuousGalerkin, NIPG, SIPG.This author was supported by the National Research Foundation of Korea,No.2014R1A2A1A11053889. − Ω + Γ Figure 1.
A domain Ω with interfaceWe prove error estimates in the mesh dependent H - norm and L - norm whichare optimal both in the order and the regularity. We carry out various numericaltests to confirm our theory and compare the performance with the unmodifiedscheme. 2. Preliminaries
Let Ω be a connected, convex polygonal domain in R which is divided into twosubdomains Ω + and Ω − by a C interface Γ = ∂ Ω + ∩ ∂ Ω − , see Figure 1. Weassume that β ( x ) is a positive function bounded below and above by two positiveconstants. Although our theory applies to the case of nonconstant β ( x ), we assume β ( x ) is piecewise constant for the simplicity of presentation: there are two positiveconstants β + , β − such that β ( x ) = β + on Ω + and β ( x ) = β − on Ω − . Consider thefollowing elliptic interface problem − ∇ · ( β ( x ) ∇ u ) = f in Ω s ( s = + , − )(2.1) u = 0 on ∂ Ω(2.2)with the jump conditions along the interface[ u ] Γ = 0 , (cid:20) β ( x ) ∂u∂n (cid:21) Γ = 0 , (2.3)where f ∈ L (Ω) and u ∈ H (Ω) and the bracket [ · ] Γ means the jump across theinterface: [ u ] Γ := u | Ω + − u | Ω − . Let p ≥ m ≥ D , we let W mp ( D ) be theusual Sobolev space with (semi)-norms and denoted by | · | m,p,D and k · k m,p,D .For m ≥
1, let f W mp ( D ) := { u ∈ W m − p ( D ) : u | D ∩ Ω s ∈ W mp ( D ∩ Ω s ) , s = + , − } , with norms; | u | p f W mp ( D ) := | u | pm,p,T ∩ Ω + + | u | pm,p,D ∩ Ω − , k u k p f W mp ( D ) := k u k pm − ,p,D + | u | p f W mp ( D ) . When p = 2, we write e H m ( D ) and denote the (semi)-norms by | u | e H m ( D ) and k u k e H m ( D ) . H (Ω) is the subspace of H (Ω) with zero trace on the boundary. Also,when some finite element triangulation {T h } is involved, the norms are understoodas piecewise norms ( P T ∈T h | u | p f W mp ( T ) ) /p and ( P T ∈T h k u k p f W mp ( T ) ) /p , etc. If p = 2, MODIFIED P - IMMERSED FINITE ELEMENT METHOD 3 we denote them by | u | m,h and k u k m,h . We also need some subspaces of e H ( T ) and e H (Ω) satisfying the jump conditions: e H ( T ) := { u ∈ H ( T ) : u | T ∩ Ω s ∈ H ( T ∩ Ω s ) , s = + , − , (cid:20) β ∂u∂n (cid:21) Γ = 0 on Γ ∩ T } e H (Ω) := { u ∈ H (Ω) : u | T ∈ e H ( T ) , ∀ T ∈ T h } . Throughout the paper, the constants
C, C , C , etc., are generic constants inde-pendent of the mesh size h and functions u, v but may depend on the problem data β, f and Ω, and are not necessarily the same on each occurrence.The usual weak formulation for the problem (2.1) - (2.3) is: Find u ∈ H (Ω)such that(2.4) Z Ω β ( x ) ∇ u · ∇ vdx = Z Ω f vdx, ∀ v ∈ H (Ω) . We have the following existence and regularity theorem for this problem; see [5],[8], [25].
Theorem 2.1.
Assume that f ∈ L (Ω) . Then the variational problem (2.4) has aunique solution u ∈ e H (Ω) which satisfies k u k e H (Ω) ≤ C k f k L (Ω) . (2.5) 3. P -immersed finite element methods We briefly review the immersed finite element space based on the P - Lagrangebasis functions ([20], [21]). Let {T h } be the usual quasi-uniform triangulations ofthe domain Ω by the triangles of maximum diameter h which may not be alignedwith the interface Γ. We call an element T ∈ T h an interface element if the interfaceΓ passes through the interior of T , otherwise we call it a noninterface element. Let T Ih be the collection of all interface elements. We assume that the interface meetsthe edges of an interface element at no more than two points.We construct the local basis functions on each element T of the partition T h . Fora noninterface element T ∈ T h , we simply use the standard linear shape functionson T whose degrees of freedom are functional values on the vertices of T , and use S h ( T ) to denote the linear spaces spanned by the three nodal basis functions on T : S h ( T ) = span { φ i : φ i is the standard linear shape function } We let S h (Ω) denote the space of usual continuous, piecewise linear polynomialswith vanishing boundary values.Now we consider a typical interface element T ∈ T Ih whose geometric configu-ration is given as in Fig. 2. Here the curve between the two points D and E is apart of the interface and DE is the line segment connecting the intersections of theinterface and the edges. A A A e e − e +3 e e E T − T + D Γ Figure 2.
A typical interface triangle
D. Y. KWAK AND JUHO LEE
We construct piecewise linear basis functions ˆ φ i , i = 1 , , φ i ( X ) = (cid:26) a + + b + x + c + y, X = ( x, y ) ∈ T + , a − + b − x + c − y, X = ( x, y ) ∈ T − ,(3.1)satisfying ˆ φ i ( A j ) = δ ij , j = 1 , , , (3.2) [ ˆ φ i ( D )] = [ ˆ φ i ( E )] = 0 , (3.3) " β ∂ ˆ φ i ∂ n DE = 0 . (3.4)These are continuous, piecewise linear functions on T satisfying the flux jumpcondition along DE , whose uniqueness and existence are known [10], [20]. Remark . Since ˆ φ i is continuous, piecewise linear, it is clear that the tangentialderivative along DE is continuous, i.e., ∂ ˆ φ + i ∂ t DE = ∂ ˆ φ − i ∂ t DE , where t DE is the tangential vector to DE .We denote by b S h ( T ) the space of functions generated by ˆ φ i , i = 1 , , immersed finite element space b S h (Ω) tobe the set of all functions φ ∈ L (Ω) such that φ ∈ b S h ( T ) if T ∈ T Ih , and φ ∈ S h ( T ) if T
6∈ T Ih , having continuity at all vertices of the triangulationand vanishes on the boundary vertices. We note that a function in b S h (Ω), in general, is not continuous across an edgecommon to two interface elements. Let H h (Ω) := H (Ω) + b S h (Ω) and equip it withthe piecewise norms | u | ,h := | u | e H (Ω) , k u k ,h := k u k e H (Ω) . Next, we define theinterpolation operator. For any u ∈ e H ( T ), we let ˆ I h u ∈ b S h ( T ) be such thatˆ I h u ( A i ) = u ( A i ) , i = 1 , , , where A i , i = 1 , , T and we call ˆ I h u the local interpolant of u in b S h ( T ). We naturally extend it to e H (Ω) by ( ˆ I h u ) | T = ˆ I h ( u | T ) for each T . Thenwe have the following approximation property [18], [21]. Proposition 3.2.
There exists a constant
C > such that X T ∈T h ( k u − ˆ I h u k ,T + h | u − ˆ I h u | ,T ) ≤ Ch k u k e H (Ω) (3.5) for all u ∈ e H (Ω) . With P - Lagrange basis function introduced in [20], [21], the IFEM reads: ( P -IFEM) Find u h ∈ b S h (Ω) such that a h ( u h , v h ) = ( f, v h ) , ∀ v h ∈ b S h (Ω) , (3.6)where a h ( u, v ) = X T ∈T h Z T β ∇ u · ∇ v dx, ∀ u, v ∈ H h (Ω) . The error estimate for this scheme is shown in [10].
MODIFIED P - IMMERSED FINITE ELEMENT METHOD 5 Modified P -IFEM In this section, we modify the P -IFEM above by adding the line integrals forjumps of fluxes and functional values. The method resembles the discontinuousGalerkin methods (see [2], [14], [24] and references therein) which use completelydiscontinuous basis functions, but the degrees of freedom in our method are muchsmaller than the DG methods since our method has the same number of basisfunctions as the conventional P -FEM.In order to describe the new method, we need some additional notations. Letthe collection of all the edges of T ∈ T h be denoted by E h and we split E h into twodisjoint sets; E h = E oh ∪ E bh , where E oh is the set of edges lying in the interior of Ω,and E bh is the set of edges on the boundary of Ω. In particular, we denote the setof edges cut by the interface Γ by E Ih . For every e ∈ E oh , there are two element T and T sharing e as a common edge. Let n T i , i = 1 , T i , but for the edge e , we choose a direction of the normalvector, say n e = n T and fix it once and for all. For functions v defined on T ∪ T ,we let [ · ] e and {·} e denote the jump and average across e respectively, i.e.[ v ] e = v − v , { v } e = 12 ( v + v ) . We also need the mesh dependent norm ||| · ||| on the space H h (Ω), ||| v ||| := X T ∈T h (cid:18)(cid:13)(cid:13)(cid:13)p βv (cid:13)(cid:13)(cid:13) ,T + (cid:13)(cid:13)(cid:13)p β ∇ v (cid:13)(cid:13)(cid:13) ,T (cid:19) + X e ∈E oh (cid:18) h (cid:13)(cid:13)(cid:13) { p β ∇ v · n e } e (cid:13)(cid:13)(cid:13) ,e + h − (cid:13)(cid:13)(cid:13) [ p βv ] e (cid:13)(cid:13)(cid:13) ,e (cid:19) . Multiplying both sides of the equation (2.1) by v ∈ H ( T ), applying Green’sformula and adding, we get X T ∈T h (cid:18)Z T β ∇ u · ∇ vdx − Z ∂T β ∇ u · n T vds (cid:19) = Z Ω f vdx. By using the preassigned normal vectors n e and adding the unharmful term ǫ R e { β ∇ v · n e } e [ u ] e for any ǫ , we see the above equation becomes(4.1) X T ∈T h Z T β ∇ u · ∇ vdx − X e ∈E oh Z e { β ∇ u · n e } e [ v ] e ds + ǫ X e ∈E oh Z e { β ∇ v · n e } e [ u ] e ds = Z Ω f vdx which is valid for v ∈ L (Ω) such that v ∈ H ( T ) for all T ∈ T h . We define thefollowing bilinear forms b ǫ ( u, v ) := − X e ∈E oh Z e { β ∇ u · n e } e [ v ] e ds + ǫ X e ∈E oh Z e { β ∇ v · n e } e [ u ] e ds,j σ ( u, v ) := X e ∈E oh Z e σh [ u ] e [ v ] e ds, for some σ > a ǫ ( u, v ) := a h ( u, v ) + b ǫ ( u, v ) + j σ ( u, v ) . Now, for each ǫ = 0, ǫ = − ǫ = 1, we define the modified P -IFEM for theproblem (2.1)-(2.3): (Modified P -IFEM) Find u mh ∈ b S h (Ω) such that(4.2) a ǫ ( u mh , v h ) = ( f, v h ) , ∀ v h ∈ b S h (Ω) . D. Y. KWAK AND JUHO LEE
This is similar to a class of DG methods, corresponding to IP, SIPG, NIPG andOBB ([1], [14], [13], [3]), if ǫ = 0, ǫ = − ǫ = 1, and ǫ = 0 , σ = 0, respectively. Remark . For the line integrals in b ǫ ( u, v ), it suffices to consider the integrals onthe edges of the interface elements only since both [ u h ] , [ v h ] vanish for e ∈ E oh \ E Ih .5. Error analysis
In this section, we prove an optimal order of error estimates in H and L -normsof our schemes. For simplicity, we present the case with ǫ = − Lemma 5.1.
There exist positive constants C , C independent of the function v h such that for all v h ∈ P k ( T ) ∪ ˆ S h ( T ) , k v h k ,T ≤ C h − k v h k ,T , k v h k ,∂T ≤ C h − k v h k ,T . (5.1) There exists a positive constant C independent of the function v such that for all v ∈ H ( T ) k v k ,e ≤ C ( h − k v k ,T + h | v | ,T ) . (5.2)Now we show the following interpolation error estimate for the mesh dependentnorm ||| · ||| . Proposition 5.2.
There exist positive constants C , C I independent of the function u such that for all u ∈ H (Ω) ∩ e H (Ω) , (5.3) X e ∈E oh h (cid:13)(cid:13)(cid:13) {∇ ( u − ˆ I h u ) · n e } e (cid:13)(cid:13)(cid:13) ,e + X e ∈E oh h − (cid:13)(cid:13)(cid:13) [ u − ˆ I h u ] e (cid:13)(cid:13)(cid:13) ,e ≤ Ch k u k e H (Ω) . Consequently, we have (5.4) ||| u − ˆ I h u ||| ≤ C I h k u k e H (Ω) . Proof.
We first consider ∇ ( u − ˆ I h u ). Since ∇ ( u − ˆ I h u ) is not in H ( T ), we cannotapply (5.2) directly. Instead, we decompose it as ∇ ( u − ˆ I h u ) = ( ∇ ( u − ˆ I h u ) · n Γ ) n Γ + ( ∇ ( u − ˆ I h u ) · t Γ ) t Γ := w + z , where n Γ and t Γ are the unit normal and tangent vector to the interface Γ, respec-tively. We have (cid:13)(cid:13)(cid:13) ∇ ( u − ˆ I h u ) · n e (cid:13)(cid:13)(cid:13) ,e ≤ k w · n e k ,e + k z · n e k ,e ≤ β min k β w · n e k ,e + k z · n e k ,e . We can easily check that β w is in H ( T ). For the smoothness of z we proceedas follows: Since u ∈ H ( T ), we have ( ∇ u · t Γ ) | T + ∩ Γ = ( ∇ u · t Γ ) | T − ∩ Γ . Hence( ∇ u · t Γ ) | T has well defined trace on Γ, which implies ∇ u · t Γ is in ˜ H ( T ) also. MODIFIED P - IMMERSED FINITE ELEMENT METHOD 7 Therefore, we can apply (5.2) to β w · n e and z · n e . Hence h (cid:13)(cid:13)(cid:13) ∇ ( u − ˆ I h u ) · n e (cid:13)(cid:13)(cid:13) ,e ≤ C β min (cid:0) k β w · n e k ,T + h | β w · n e | ,T (cid:1) + C (cid:0) k z · n e k ,T + h | z · n e | ,T (cid:1) ≤ C (cid:18) β max β min (cid:16) k w k ,T + h | w | H ( T ) (cid:17) + k z k ,T + h | z | H ( T ) (cid:19) ≤ C α (cid:18)(cid:13)(cid:13)(cid:13) ∇ ( u − ˆ I h u ) (cid:13)(cid:13)(cid:13) ,T + h (cid:12)(cid:12)(cid:12) ∇ ( u − ˆ I h u ) (cid:12)(cid:12)(cid:12) H ( T ) (cid:19) ≤ Ch k u k H ( T ) , where β min = min( β + , β − ) , β max = max( β + , β − ) and we have set α = β max β min . HereProposition 3.2 was used to derive the last estimate.The estimate of the second term follows easily from (5.2) and Proposition 3.2: h − (cid:13)(cid:13)(cid:13) u − ˆ I h u (cid:13)(cid:13)(cid:13) ,∂T ≤ C (cid:18) h − (cid:13)(cid:13)(cid:13) u − ˆ I h u (cid:13)(cid:13)(cid:13) ,T + (cid:12)(cid:12)(cid:12) u − ˆ I h u (cid:12)(cid:12)(cid:12) H ( T ) (cid:19) ≤ Ch k u k H ( T ) . Thus, the estimate (5.4) follows. (cid:3)
The following discrete Poincar´e inequality holds for ˆ S h (Ω), (see [10]). Lemma 5.3.
There exists a constant C p > such that (5.5) C p k v h k , Ω ≤ | v h | ,h , ∀ v h ∈ ˆ S h (Ω) . Now we show some basic properties of a ǫ ( · , · ). Clearly, a ǫ ( · , · ) is bounded on H h (Ω) with respect to ||| · ||| : | a ǫ ( u, v ) | ≤ C b ||| u |||||| v ||| , ∀ u, v ∈ H h (Ω) . Next, we prove the coercivity of the form a ǫ ( · , · ) on the space ˆ S h (Ω). We need alemma. Lemma 5.4.
For all v ∈ ˆ S h (Ω) , there exists a positive constant C independent ofh such that (5.6) X e ∈E oh h (cid:13)(cid:13)(cid:13)np β ∇ v · n o(cid:13)(cid:13)(cid:13) ,e ≤ C α X T ∈T h (cid:13)(cid:13)(cid:13)p β ∇ v (cid:13)(cid:13)(cid:13) ,T where α = β max β min is the same as before.Proof. We decompose ∇ v h as(5.7) ∇ v h = ( ∇ v h · n Γ ) n Γ + ( ∇ v h · t Γ ) t Γ := w + z . The rest of the proof is almost the same as that of (5.3). (cid:3)
Proposition 5.5.
There exists a positive constant C c independent of v h such thatfor all v h ∈ ˆ S h (Ω) the following holds: a ǫ ( v h , v h ) ≥ C c ||| v h ||| . D. Y. KWAK AND JUHO LEE
Proof.
First of all, we consider b ǫ ( v h , v h ), the second part of a ǫ ( v h , v h ). By Lemma5.4, Cauchy-Schwarz and arithmetic-geometric inequality, we have X e ∈E oh Z e { β ∇ v h · n } [ v h ] ds ≤ X e ∈E oh h k{ β ∇ v h · n }k ,e / X e ∈E oh h − k [ v h ] k ,e / ≤ C α X T ∈T h (cid:13)(cid:13)(cid:13)p β ∇ v h (cid:13)(cid:13)(cid:13) ,T ! / X e ∈E oh h − (cid:13)(cid:13)(cid:13) [ p βv h ] (cid:13)(cid:13)(cid:13) ,e / ≤ γ X T ∈T h (cid:13)(cid:13)(cid:13)p β ∇ v h (cid:13)(cid:13)(cid:13) ,T ! + C α γ X e ∈E oh h − (cid:13)(cid:13)(cid:13) [ p βv h ] (cid:13)(cid:13)(cid:13) ,e for every γ >
0. Hence by Lemma 5.3, we have a ǫ ( v h , v h ) = a h ( v h , v h ) + b ǫ ( v h , v h ) + j σ ( v h , v h )= X T ∈T h Z T β ∇ v h · ∇ v h dx − X e ∈E oh Z e { β ∇ v h · n } [ v h ] ds + X e ∈E oh Z e σh [ v h ] ds ≥ C p (cid:13)(cid:13)(cid:13)p βv h (cid:13)(cid:13)(cid:13) , Ω + ( 12 − γ ) (cid:12)(cid:12)(cid:12)p βv h (cid:12)(cid:12)(cid:12) ,h + (cid:18) σ − C αγ (cid:19) X e ∈E oh h (cid:13)(cid:13)(cid:13) [ p βv h ] (cid:13)(cid:13)(cid:13) ,e ≥ C p (cid:13)(cid:13)(cid:13)p βv h (cid:13)(cid:13)(cid:13) , Ω + ( 14 − γ ) (cid:12)(cid:12)(cid:12)p βv h (cid:12)(cid:12)(cid:12) ,h + 14 C α X e ∈E oh h (cid:13)(cid:13)(cid:13)np β ∇ v h · n o(cid:13)(cid:13)(cid:13) ,e + (cid:18) σ − C αγ (cid:19) X e ∈E oh h (cid:13)(cid:13)(cid:13) [ p βv h ] (cid:13)(cid:13)(cid:13) ,e , where we have set σ = σ/β . If we choose γ = and σ large enough so that( σ − C α ) ≥ . Then with C c := min (cid:16) C p , , C α (cid:17) , we have a ǫ ( v, v ) ≥ C c ||| v ||| . (cid:3) Remark . We can take any positive σ when ǫ = 1, because b ǫ ( v, v ) becomes zero.If ǫ = 0 or −
1, it seems that σ > σ or even σ = 0 works for all the cases we have tested.This is in contrast to the usual DG schemes, where sufficiently large σ is necessary.The reason seems to be that, unlike the usual DG, the term b ǫ ( v, v ) is small enoughto be dominated by a h ( v, v ), since the jump [ v ] vanishes at the vertices of each T ∈ T h . In fact, using the techniques in [9] and the proof of Proposition 5.2 we canshow | b ǫ ( v, v ) | ≤ C ( h | log h | ) / k v k ,h , but the details are complicated. This willbe shown in the subsequent paper.5.1. H -error analysis. First we check that the modified P -IFEM is consistent. Lemma 5.7.
Let u be the solution of (2.1)-(2.3) and let u mh be the solution of(4.2). For any v h ∈ b S h (Ω) , we have (5.8) a ǫ ( u, v h ) = ( f, v h ) . In other words, a ǫ ( u − u mh , v h ) = 0 . MODIFIED P - IMMERSED FINITE ELEMENT METHOD 9 Proof.
By (4.1), the definition of the a ǫ form and the homogeneous jump conditionof u , we have a ǫ ( u, v h ) − a ǫ ( u mh , v h ) = a ǫ ( u, v h ) − ( f, v h ) = X e ∈E oh Z e σh [ u ] e [ v h ] e ds = 0 . (cid:3) Now we can prove the H -error estimate which is optimal both in order and theregularity. Theorem 5.1.
Let u be the solution of (2.1)-(2.3) and let u mh be the solution of(4.2). Then there exists a positive constant C independent of u and h such that ||| u − u mh ||| ≤ Ch k u k e H (Ω) . Proof.
By Proposition 5.5, (5.8) and boundedness of a ǫ ( · , · ) with respect to ||| · ||| ,we have ||| u mh − ˆ I h u ||| ≤ C − c a ǫ ( u mh − ˆ I h u, u mh − ˆ I h u )= C − c a ǫ ( u − ˆ I h u, u mh − ˆ I h u ) ≤ C − c C b ||| u − ˆ I h u |||||| u mh − ˆ I h u ||| . By the triangle inequality and Proposition 5.2, we get ||| u − u mh ||| ≤ ||| u − ˆ I h u ||| + ||| u mh − ˆ I h u |||≤ ( C − c C b + 1) C I h k u k e H (Ω) . (cid:3) L -error analysis.Theorem 5.2. For the solution u mh of (4.2), there exists a positive constant C independent of u and h such that k u − u mh k L (Ω) ≤ Ch k u k e H (Ω) . Proof.
Consider the dual equation: −∇ ( β ∇ Ψ) = w in Ω s ( s = + , − )[Ψ] Γ = 0 , (cid:20) β ( x ) ∂ Ψ ∂n (cid:21) Γ = 0 , Ψ = 0 on ∂ Ω . Then by Theorem 2.1 the solution satisfies(5.9) k Ψ k ˜ H (Ω) ≤ C k w k L (Ω) . Let Ψ h be the modified IFEM solution of this problem. Then with e h := u − u mh ,we have by Lemma 5.7( e h , w ) = a ǫ ( e h , Ψ) = a ǫ ( e h , Ψ − Ψ h ) . Then by boundedness of a ǫ , Theorem 5.1 and (5.9) | ( e h , w ) | ≤ C ||| e h |||||| Ψ − Ψ h ||| ≤ Ch k Ψ k ˜ H (Ω) ||| e h ||| ≤ Ch k w k L (Ω) ||| e h ||| . Taking w = e h , we obtain k u − u mh k L (Ω) ≤ Ch ||| u − u mh ||| ≤ Ch k u k e H (Ω) . (cid:3) /h x k u − u h k order k u − u h k ,h order k u − u h k ∞ order8 1.344e-2 3.315e-1 2.761e-216 3.453e-3 1.961 1.709e-1 0.955 8.715e-3 1.663 P -IFEM 32 8.9002-4 1.956 8.727e-3 0.970 3.069e-3 1.50664 2.161e-4 2.043 4.507e-2 0.953 1.295e-3 1.245128 5.541e-5 1.963 2.347e-2 0.941 5.786e-4 1.162256 1.851e-5 1.582 1.288e-2 0.865 3.598e-4 0.686512 8.193e-6 1.176 7.297e-3 0.820 1.776e-4 1.0181 /h x k u − u mh k order k u − u mh k ,h order k u − u mh k ∞ order8 1.233e-2 3.306e-1 2.345e-216 3.260e-3 1.919 1.694e-1 0.965 6.765e-3 1.793Modified 32 8.269e-4 1.979 8.554e-2 0.986 1.775e-3 1.931 P -IFEM 64 2.094e-4 1.982 4.300e-2 0.992 4.621e-4 1.941128 5.286e-5 1.986 2.156e-2 0.996 1.185e-4 1.964256 1.328e-5 1.993 1.078e-2 0.999 2.991e-5 1.986512 3.308e-6 2.005 5.399e-3 0.998 7.557e-6 1.985 Table 1.
Example 6.1 (Cubic curve): β − = 1 , β + = 106. Numerical Experiments
For numerical tests, we solve the problem (2.1)-(2.3) on the rectangular domainΩ = [ − , × [ − ,
1] partitioned into unform right triangles with h x = h y = 1 / n − for n = 4 , · · · ,
10. Three types of interface problems are considered with variousvalues of parameter β . We measured k u − u h k and k u − u h k ,h which are veryclose to the theoretical orders of convergence, 2 and 1 respectively. Although notreported, we also measured P e k u − u h k ,e and P e k ∂ ( u − u h ) /∂n k ,e , the ordersof which agree with the theoretical value 1 . . L ∞ norm also. Ω − Ω + y = 3 x ( x − . x − .
8) + 0 .
34 Ω + Ω − θ y = (( x − .
6) tan θ ) ( x + 0 .
4) Ω − Ω + x / . + y / . = 1 Figure 3.
Interfaces of Examples 1,2 and 3
Example 6.1 (Cubic curve) . The 0-set of function L ( x, y ) = y − x ( x − . x − . − . is used in this example as the interface. The exact solution is u = L ( x, y ) /β , where β = β ± on Ω ± . We test the cases when β + /β − = 10 and . The comparison with error surfaces in Figure ?? . shows that modified methodgives much more accurate results than the original P -IFEM when β + /β − = 10and 1 /h x = 128. The smaller the mesh, the more accurate results the modifiedmethod shows.Table 1 shows the comparison of errors between the two methods when β + /β − =10. We can see the original P -IFEM has suboptimal convergence as the grids arerefined (1 /h x = 256 and 512). However, the modified method shows a robust orderof convergence for all grids. MODIFIED P - IMMERSED FINITE ELEMENT METHOD 11 /h x k u − u h k order k u − u h k ,h order k u − u h k ∞ order8 1.923e-2 3.530e-1 5.617e-216 4.002e-3 2.264 1.716e-1 1.040 1.470e-2 1.934 P -IFEM 32 9.196e-4 2.122 8.453e-2 1.022 3.854e-3 1.93264 2.291e-4 2.005 4.221e-2 1.002 1.288e-3 1.582128 5.408e-5 2.083 2.105e-2 1.004 2.836e-4 2.183256 1.337e-5 2.016 1.056e-2 0.995 1.159e-4 1.291512 3.336e-6 2.002 5.304e-3 0.994 5.258e-5 1.1411 /h x k u − u mh k order k u − u mh k ,h order k u − u mh k ∞ order8 1.266e-2 3.216e-1 2.470e-216 3.205e-3 1.982 1.643e-1 0.969 6.836e-3 1.854Modified 32 8.163e-4 1.973 8.293e-2 0.986 1.784e-3 1.938 P -IFEM 64 2.068e-4 1.981 4.172e-2 0.991 4.642e-4 1.943128 5.199e-5 1.992 2.093e-2 0.996 1.185e-4 1.970256 1.302e-5 1.998 1.048e-2 0.998 3.009e-5 1.977512 3.259e-6 1.998 5.243e-3 0.999 7.564e-6 1.992 Table 2.
Example 6.1 (Cubic curve) : β − = 1 , β + = 1000 /h x k u − u h k order k u − u h k ,h order k u − u h k ∞ order8 3.359e-3 7.958e-2 1.036e-216 9.014e-4 1.898 4.185e-3 0.927 4.118e-3 1.332 P -IFEM 32 2.219e-4 2.022 2.161e-3 0.954 1.958e-3 1.07364 5.686e-5 1.965 1.197e-3 0.852 9.568e-4 1.033128 1.463e-5 1.958 6.573e-3 0.865 5.063e-4 0.918256 6.070e-6 1.269 3.967e-3 0.728 2.462e-4 1.040512 2.942e-6 1.045 2.439e-3 0.702 1.241e-4 0.9881 /h x k u − u mh k order k u − u mh k ,h order k u − u mh k ∞ order8 3.056e-3 7.817e-2 9.005e-316 7.441e-4 2.038 3.956e-2 0.983 2.316e-3 1.959Modified 32 1.930e-4 1.947 1.990e-2 0.991 6.221e-4 1.896 P -IFEM 64 4.716e-5 2.033 1.000e-2 0.993 1.608e-4 1.952128 1.216e-5 1.956 5.015e-3 0.996 4.090e-5 1.975256 3.010e-6 2.014 2.510e-3 0.999 1.031e-5 1.989512 7.621e-7 1.982 1.256e-3 0.999 2.633e-6 1.968 Table 3.
Example 6.2 (Sharp corner) : θ = 45 ◦ , β − = 1 , β + = 10On the other hands, Table 2 shows both methods has an optimal convergence in L and H norms when β + /β − = 1000. Remark . Comparing Tables 1 and 2, we see the original P -IFEM behavesbetter when β + /β − = 1000 than β + /β − = 10. This is a common phenomenon forall the examples we tested. This seems to contradict the usual behavior of standardFEM. We guess the reason is that the large ratio between the coefficients masks thediscontinuity of basis functions. Figure ?? shows the behavior of P -IFEM basisbetween β + /β − = 10 and β + /β − = 1000. When β + /β − = 10, the gap betweenadjacent elements is conspicuous. However, when β + /β − = 1000, the gap is almostinvisible. Example 6.2 (Sharp corner) . In this example, we consider an interface with asharp corner having interior angle θ . Let Γ be the zero set of L ( x, y ) = − y +(( x − .
6) tan θ ) ( x +0 . for x ≤ . . We test the case with θ = 45 ◦ and β + /β − = 10 and / . The exact solution is u = L ( x, y ) /β . /h x k u − u h k order k u − u h k ,h order k u − u h k ∞ order8 1.238e-2 3.013e-1 1.613e-216 3.159e-3 1.971 1.513e-1 0.994 4.327e-3 1.899 P -IFEM 32 7.949e-4 1.991 7.572e-2 0.998 1.174e-3 1.88264 2.030e-4 1.969 3.821e-2 0.987 7.475e-4 0.651128 5.366e-5 1.920 1.933e-2 0.983 4.704e-4 0.668256 1.528e-5 1.812 9.919e-3 0.963 2.452e-4 0.940512 4.898e-6 1.642 5.155e-3 0.944 1.199e-4 1.0331 /h x k u − u mh k order k u − u mh k ,h order k u − u mh k ∞ order8 1.238e-2 3.010e-1 1.610e-216 3.094e-3 2.000 1.507e-1 0.998 4.107e-3 1.971Modified 32 7.787e-4 1.990 7.543e-2 0.999 1.037e-3 1.986 P -IFEM 64 1.947e-4 2.000 3.773e-2 0.999 2.605e-4 1.993128 4.876e-5 1.998 1.887e-2 1.000 6.528e-5 1.997256 1.219e-5 2.000 9.435e-3 1.000 1.634e-5 1.998512 3.051e-6 1.998 4.718e-3 1.000 4.087e-6 1.999 Table 4.
Example 6.2 (Sharp corner) : θ = 45 ◦ , β − = 10 , β + = 1 /h x k u − u h k order k u − u h k ,h order k u − u h k ∞ order8 8.550e-2 1.585e-0 2.415e-116 2.931e-2 1.544 9.840e-1 0.688 1.025e-1 1.237 P -IFEM 32 7.954e-3 1.882 5.538e-1 0.829 4.174e-2 1.29564 2.002e-3 1.990 3.033e-1 0.869 1.568e-2 1.413128 4.825e-4 2.053 1.665e-1 0.865 8.471e-3 0.888256 1.206e-4 2.000 8.948e-2 0.896 4.393e-3 0.947512 3.461e-5 1.802 5.063e-2 0.822 2.132e-3 1.0431 /h x k u − u mh k order k u − u mh k ,h order k u − u mh k ∞ order8 8.652e-2 1.572e-0 2.150e-116 2.867e-2 1.593 9.704e-1 0.696 9.448e-2 1.187Modified 32 8.049e-3 1.833 5.368e-1 0.854 3.656e-2 1.370 P -IFEM 64 2.195e-3 1.874 2.889e-1 0.894 1.097e-2 1.736128 5.585e-4 1.975 1.485e-1 0.959 3.055e-3 1.845256 1.437e-4 1.958 7.550e-2 0.976 8.386e-4 1.865512 3.649e-5 1.978 3.809e-2 0.987 2.209e-4 1.925 Table 5.
Example 6.3 (Variable coefficient)This example is not covered by analysis of this work because the problem haslow regularity at the interface corner. However, we see that the modified methodworks better; See the Table 3.
Remark . We have also computed other cases such as β + /β − = 1000 , .
001 withvarious angles. The results of our scheme are always optimal while the unmodified P immersed method deteriorates for some cases. Example 6.3 (Variable coefficient) . Finally, we consider the case with variablecoefficient. The 0-set of function L ( x, y ) = x / (0 . + y / (0 . − . is used inthis example as the interface. The exact solution is u = L ( x, y ) /β ( x, y ) where β ( x, y ) = (cid:26) ( x + y − on Ω − , on Ω + . In this case, both methods show an optimal order of convergence in H -norm. Butthe modified method performs much better in L ∞ -norm; See Table 5. MODIFIED P - IMMERSED FINITE ELEMENT METHOD 13 Conclusion
We introduced a modified IFEM for solving elliptic interface problems. Byadding the line integral terms similar to the DG methods, we overcome the sub-optimal behavior of the original IFEM proposed in [20], [21]. (The computationalresult there seemed to show optimal order. However, more numerical experimentsshow the original P -IFEM is not optimal for some problems. The proof in [10]seems incorrect. However, our modified scheme is always robust for all problemstested including unreported ones). The optimal convergence rates in H and L norms are shown by a similar technique as in DG methods; the modified IFEM isconsistent, the coercivity and boundedness of the bilinear form hold. Several nu-merical tests show the errors are O ( h ), O ( h ) order in respective norms. Althoughno proof is given, we also obtain O ( h ) order in L ∞ norm.Some of the limitation of our scheme might be these: for problems with sharpinterface, one has to arrange the grids so that the cusp point is located at a vertexof an element, for problems with highly oscillating interface, further refinement arenecessary to apply the IFEM.We now comment on the computational aspects: The matrix structure are ex-actly the same as usual P -FEM, i.e., 5-point stencil; the number of unknowns arealso the same. When ǫ = −
1, the scheme becomes symmetric. The assembly ofstiffness matrix requires slightly more time than the unmodified P -IFEM, but thetime for iterative solver such as conjugate gradient is almost the same.An obvious advantage of (both) IFEM is that we can use fast solver such asmultigrid methods since we can use uniform meshes.Future works related to this topic are:(1) Local refinement near singularity.(2) Problems with nonhomogenous jumps, tensor coefficients, etc.(3) Q -IFEM for rectangular elements.(4) 3-dimensional problems.(5) Problems with a moving interface.(6) Two phase Stokes/Navier-Stokes problems.(7) Development of fast solver such as multigrid methods. References [1] D. N. Arnold,
An interior penalty finite element method with discontinuous elements , SIAMJ. Numer. Anal., 19 (1982), pp. 742–760.[2] D. Arnold, F. Brezzi, B. Cockburn, and D. Marini,
Discontinuous Galerkin methods forelliptic problems , in Discontinuous Galerkin Methods. Theory, Computation and Applica-tions, B. Cockburn, G. E. Karniadakis, and C.-W. Shu, eds., Lecture Notes in Comput. Sci.Engrg. 11, Springer-Verlag, NewYork, 2000, pp. 89–101.[3] C. E. Baumann and J. T. Oden,
A discontinuous hp finite element method for convection-diffusion problems , Comput. Meth. Appl. Mech. Engrg., 175 (1999), pp. 311–341.[4] D. Braess,
Finite elements: Theory, fast solvers, and applications in solid mechanics , Sec-ond edition. Cambridge University Press, Cambridge, 2001.[5] J. H. Bramble and J. T. King,
A finite element method for interface problems in domainswith smooth boundary and interfaces , Adv. Comp. Math. (1996), pp. 109-138.[6] S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods , Springer-Verlag, New York, (1994).[7] K. S. Chang and Do Y. Kwak,
Discontinuous Bubble scheme for elliptic problems withjumps in the solution , Comp. Meth. Appl. Mech. Engrg., 200 (2011) pp. 494–508.[8] Z. Chen and J. Zou,
Finite element methods and their convergence for elliptic and parabolicinterface problems , Numer. Math. (1998), pp. 175–202. fitted grid method.[9] S. H. Chou, D. Y. Kwak and K. T. Wee, Some error estimates for an immersed interfacefinite element method , unpublished manuscript (submitted to Computer Methods in AppliedMechanics and Engineering), December 6, 2006. [10] S. H. Chou, D. Y. Kwak and K. T. Wee,
Optimal convergence analysis of an immersedinterface finite element method , Adv Comput Math, V. 33 (2010), pp. 149–168.[11] P. G. Ciarlet,
The finite element method for elliptic problems , North Holland, 1978.[12] M. Crouzeix and P. A. Raviart,
Conforming and nonconforming finite element methods forsolving the stationary Stokes equations , RAIRO Anal. Num´er. (1973), pp. 33-75.[13] C. Dawson, S. Sun, and M. F. Wheeler,
Compatible algorithms for coupled flow and trans-port , Comput. Meth. Appl. Mech. Eng., 194 (2004), pp. 2565- 2580.[14] J. Douglas, Jr. and T. Dupont,
Interior Penalty Procedures for Elliptic and ParabolicGalerkin Methods , Lecture Notes in Phys. 58, Springer-Verlag, Berlin, 1976.[15] A. Hansbo and P. Hansbo,
An unfitted finite element method, based on Nitsche’s method,for elliptic interface problems , Comput. Methods Appl. Mech. Engrg. 191, (2002), pp. 5537–5552.[16] S. Hou and X. Liu,
A numerical method for solving variable coefficient elliptic equation withinterfaces , J. Comput. Phys. (2005), no. 2, pp. 411-445.[17] S. Hou, L. Wang and W. Wang,
Numerical method for solving matrix coefficient ellipticequation with sharp-edged interfaces , J. Comput. Phys. (2010), no. 19, pp. 7162-7179.[18] Do Y. Kwak, K. T. Wee and K. S. Chang,
An analysis of a broken P -nonconforming finiteelement method for interface problems , SIAM J. Numer. Anal. (2010), pp. 2117–2134[19] M. Lai, Z. Li and X. Lin, Fast solvers for 3D Poisson equations involving interfaces in afinite or the infinite domain , J. Comput. Appl. Math. (2006), no. 1, pp. 106-125.[20] Z. Li, T. Lin and X. Wu,
New Cartesian grid methods for interface problems using the finiteelement formulation , Numer. Math. (2003), pp. 61-98.[21] Z. Li, T. Lin, Y. Lin and R. C. Rogers, An immersed finite element space and its approxi-mation capability , Numer. Methods. Partial Differential Equations (2004), pp. 338-367.[22] Z. Li and K. Ito, The immnersed interface method: Numerical solutions of PDEs involvinginterfaces and irregular domains , Frontiers in Applied Mathematics , SIAM, 2006.[23] M. Oevermann, C. Scharfenberg, R. Klein, A sharp interface finite volume method forelliptic equations on Cartesian grids , J. Comput. Phys. (2009) No. 14, pp. 5184-5206.[24] B. Rivi´ere, M. F. Wheeler, and V. Girault,
Improved energy estimates for interior penalty,constrained and discontinuous Galerkin methods for elliptic problems I , Comput. Geosci.,3 (1999), pp. 337–360.[25] J. A. Roitberg and Z. G. Seftel,
A theorem on homeomorphisms for elliptic systems and itsapplications , Math. USSR-Sb. 7, (1969), pp. 439-465.[26] S. Yu, Y. Zhou, and G.W. Wei
Matched interface and boundary (MIB) method for ellipticproblems with sharp-edged interfaces , J. Comp. Physics, (2007), pp. 729–756.
Department of Mathematical Sciences, KAIST,, 291 Daehak-ro, Yuseong-gu, Dae-jeon, Korea 305-701.
E-mail address : [email protected] Department of Mathematical Sciences, KAIST,, 291 Daehak-ro, Yuseong-gu, Dae-jeon, Korea 305-701.
E-mail address ::