An Analysis of broken P_{1}-Nonconforming Finite Element Method For Interface Problems

# An Analysis of broken P1-Nonconforming Finite Element Method For Interface Problems

Do Y. Kwak Korea Advanced Institute of Science and Technology, Daejeon, Korea 305-701. email:kdy@kaist.ac.kr. This work was supported by a grant from Korea Science and Engineering Foundation.(R01-2007-000-10062-0)    K. T. Wee Korea Advanced Institute of Science and Technology, Daejeon, Korea 305-701 email:ktwee@kaist.ac.kr .
###### Abstract

We study some numerical methods for solving second order elliptic problem with interface. We introduce an immersed interface finite element method based on the ‘broken’ -nonconforming piecewise linear polynomials on interface triangular elements having edge averages as degrees of freedom. This linear polynomials are broken to match the homogeneous jump condition along the interface which is allowed to cut through the element. We prove optimal orders of convergence in and -norm. Next we propose a mixed finite volume method in the context introduced in  using the Raviart-Thomas mixed finite element and this ‘broken’ -nonconforming element. The advantage of this mixed finite volume method is that once we solve the symmetric positive definite pressure equation(without Lagrangian multiplier), the velocity can be computed locally by a simple formula. This procedure avoids solving the saddle point problem. Furthermore, we show optimal error estimates of velocity and pressure in our mixed finite volume method. Numerical results show optimal orders of error in -norm and broken -norm for the pressure, and in -norm for the velocity.

Key words. Immersed interface, -nonconforming finite element method, uniform grid, mixed finite volume method, average degrees of freedom

AMS(MOS) subject classifications. 65N15, 65N30, 35J60.

## 1 Introduction

There are many physical problems where the underlying partial differential equations have an interface. For example, second order elliptic equations with discontinuous coefficients are often used to model problems in material sciences and porous media when two or more distinct materials or media with different conductivities, densities or permeability are involved. The solution of these interface problems must satisfy interface jump conditions due to conservation laws. If the interface is smooth enough, then the solution of the interface problem is also smooth in individual regions where the coefficient is smooth, but due to the jump of the coefficient across the interface, the global regularity is usually low and the solution usually belongs to for some . Because of the low global regularity, achieving accuracy is difficult with standard finite element methods, unless the elements fit with the interface of general shape.

Immersed interface method using uniform grids has many advantages over usual fitted grid method. Using uniform grid, one does not need to generate a grid. This is quite convenient in several aspects. First of all, the structure of stiffness matrix is the same as that of the standard finite element method, where many known efficient solvers can be exploited. Second, when a moving interface problem is involved one does not need to generate a new grid as time evolves. This saves considerable amount of time and storage.

The first attempt to avoid fitted grid for interface problem was made by LeVeque and Z. Li , where they proposed an immersed interface method for finite difference method where the jump condition was properly incorporated in the scheme. Cartesian grids is most natural in this case. They subsequently applied the same idea to other interface problems such as the Stokes flow problem, one-dimensional moving interface problem and Hele-Shaw flow, etc. [26, 31, 32, 33] The resulting linear systems from these methods are non-symmetric and indefinite even when the original problem is self-adjoint and uniformly elliptic. Although these methods were demonstrated to be very effective, convergence analysis of related finite difference methods are extremely difficult and are still open.

For finite element methods, T. Lin et. al. [35, 36, 37, 38] recently studied an immersed interface finite element method using uniform grids and they proved the approximation property of the finite element space of their scheme. Their numerical examples demonstrated optimal orders of the error. Other related works in this direction can be found in [10, 25, 29, 34, 39, 43] and references therein.

On the other hand, -nonconforming finite element method introduced in  for solving Stokes equation is being widely used in solving elliptic equations and shown to be quite effective [15, 16, 20, 27]. Especially, it is extremely useful in solving mixed finite element method by hybridization [1, 2] or finite volume formulation [15, 16, 17, 19].

The mixed finite element method based on the dual formulation is well-known [5, 6, 7, 8, 11, 21, 23, 24, 42]. The motivation of mixed method is to obtain an accurate approximation of the flow variable and has been widely used in the study of flow in porous media such as petroleum engineering, underground water flow, and electrodynamics, etc. (e.g., [12, 22, 41]) But this scheme leads to a saddle point problem for which many well known fast iterative methods fails. To overcome this difficulty, mixed hybrid methods have been introduced [2, 8, 40], where the problem reduces to a symmetric positive definite system in Lagrange multiplier only. The flow and pressure variables are obtained via some post -processing.

Recently, there has been some development of the mixed finite element method in another direction: A mixed finite volume method was proposed in  and extended in [15, 16]. In this method, one use Raviart-Thomas space and -nonconforming space as trial spaces for velocity and pressure, and integrates the mixed system of equation on each volume. Then one can eliminate velocity variable and obtain the equation of pressure variable only (in terms of -nonconforming FEM) directly from the formulation without using Lagrange multiplier. The resulting linear system is again symmetric positive definite, and velocity can be recovered from pressure locally in a simple manner.

The purpose of this paper is two-folded. First, we propose a finite element method on a uniform triangular grid using ‘broken’ -functions having degrees of freedom on edges. This is a Galerkin type -nonconforming finite element method with the basis functions having the average on edges as degrees of freedom, broken along the interface to match the flux condition. Then we show optimal error estimates in -norms. Here, we emphasize that the meaning of ‘nonconforming’ is different from the context of Li et al. [37, 35] where the basis function has degrees of freedom at vertices, discontinuous along edges of interface elements. Meanwhile, the basis function here are Crouzeix-Raviart type. Hence it is discontinuous along all edges intrinsically. Furthermore, since we use the average of linear function(possibly broken) along edges as degrees of freedom, the overhead of dealing with nonconformity in the proof of error estimate is significantly reduced.(See section 3)

Next, we propose a mixed finite volume method using Raviart-Thomas space and the immersed interface finite element introduced above. This is similar to the scheme studied in , but the usual nonconforming basis function is replaced by a ‘broken’ one on the interface element. We provide an optimal error analysis of pressure and velocity.

The rest of the paper is organized as follows. In the next section, we will describe the model problem and some preliminaries. We construct an immersed interface -nonconforming space with average degrees of freedom which preserves flux continuity weakly along the interface, and prove an interpolation error estimate. In Sections 3 and 4, we propose an immersed interface finite element scheme and prove and -error estimates. In Section 5, we propose a mixed finite volume method using Raviart-Thomas mixed finite element and our -nonconforming immersed interface finite element method, where the problem reduces to symmetric positive definite system in pressure variables. The velocity can be computed locally after pressure computation. Finally, in Section 6, some numerical results are presented which indicate optimal orders convergence of our methods.

## 2 Preliminaries

Let be a convex polygonal domain in which is separated into two sub-domains and by a -interface with as in Fig. 1. We consider the following elliptic interface problem

 {−div(β∇p)=f % in Ω∖Γ,p=0 on ∂Ω (2.1)

with the jump conditions on the interface

 [p]=0,   [β∂p∂n]=0   across Γ, (2.2)

where and . We assume that the coefficient is positive and piecewise constant, that is,

We take as usual the weak formulation of the interface problem: Find such that

 ∫Ωβ∇p⋅∇qdx=∫Ωfqdx,   ∀q∈H10(Ω). (2.3)

Now we introduce the space

 ˜H2(Ω):={p∈H1(Ω):p∈H2(Ωs),s=+,−}

equipped with the norm

 ∥p∥2˜H2(Ω):=∥p∥2H1(Ω)+∥p∥2H2(Ω+)+∥p∥2H2(Ω−),  ∀p∈˜H2(Ω),

where is the usual Sobolev space of order . By Sobolev embedding theorem, for any , we have . Then we have the following regularity theorem for the weak solution of the variational problem (2.3); see  and .

{theorem}

The variational problem (2.3) has a unique solution which satisfies for some constant

 ∥p∥˜H2(Ω)≤C∥f∥L2(Ω). (2.4)

We now describe an immersed interface finite element method with piecewise -nonconforming functions.

For the simplicity of presentation, we assume that is a rectangular domain. First we consider uniform rectangular partitions of mesh size . Then we obtain triangular partitions by cutting the elements along diagonals. Thus we allow the interface to cut through the elements. We assume the following situation: The interface

• meets the edges of an interface element at no more than two points.

• meets each edge at most once except possibly it passes through two vertices.

These assumptions are reasonable if we choose sufficiently small.

We call an element an interface element if the interface passes through the interior of , otherwise we call a non-interface element. (If one of the edges is part of the interface, then the element is a non-interface element.) Let be a collection of all edges of .

Let be the line segment connecting the intersections of the interface and the edges of a triangle . This line segment divides into two parts and with . There is a small region in T such that (see figure 2). Since can be considered as an approximation of the -curve , the interface is perturbed by a term. From [4, 13], one can see for the interpolation polynomial defined below, such a perturbation will only affect the interpolation error to the order of .

As usual, we want to construct local basis functions on each element of the partition . For a non-interface element , we simply use the standard linear shape functions on having degrees of freedom at the mid-points of the edges, and use to denote the linear spaces spanned by the three nodal basis functions on : Let be the midpoints of edges of . Then

 ¯¯¯¯Sh(T)=span{ϕi:ϕi is% linear on T and ϕi(mj)=δij,i,j=1,2,3} (2.5)

Alternatively, we can use average values along edges of as degrees of freedom, i.e., can defined by .

For this space, we have the following well-known approximation property [18, 20]:

 ∥p−Ihp∥L2(T)+h∥p−Ihp∥H1(T)≤Ch2∥p∥H2(T), (2.6)

where is the interpolation operator. Finally, we use to denote the space of the standard piecewise -nonconforming space with vanishing boundary nodal values.

### 2.1 Local basis functions on an interface element

We now consider a typical reference interface element whose geometric configuration is given in Fig. 2 in which the curve between points and is part of the interface. Let be the edges of . For , let denote the average of along , i.e., .

We construct a piecewise linear function of the form

 ϕ(X)={ϕ+(X)=a0+b0x+c0y, X=(x,y)∈T+,ϕ−(X)=a1+b1x+c1y, X=(x,y)∈T−, (2.7)

satisfying

 ¯ϕei=Vi,  i=1,2,3, (2.8) ϕ+(D)=ϕ−(D), ϕ+(E)=ϕ−(E), β+∂ϕ+∂n¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny{DE}=β−∂ϕ−∂n¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny% {DE}, (2.9)

where are given values and is the unit normal vector on the line segment . This is a piecewise linear function on that satisfies the homogeneous jump conditions along .

Suppose that a typical reference interface element has vertices at . We assume that the interface meets with the edges at and where . Then the unit normal vector to the interface is .

{theorem}

Given an reference interface triangle, the piecewise linear function defined by (2.7)-(2.9) is uniquely determined by three conditions

 ¯ϕei=Vi,i=1,2,3.
{proof}

Let . Since and are linear functions, we have

 ϕ(X)={ϕ+(X)=a0+b0x+c0y, X∈T+,ϕ−(X)=a1+b1x+c1y, X∈T−. (2.10)

The condition (2.8) gives the following three equations:

 ¯ϕe1 = a0+12b0+12c0=V1 (2.11)

and

 ¯ϕe2 = ∫e2ϕds=∫¯¯¯¯¯¯¯AEϕ−ds+∫¯¯¯¯¯¯¯ECϕ+ds (2.12) = (a1+y02c1)y0+(a0+y0+12c0)(1−y0)=V2,

where we used mid-point quadrature on and . Similarly, we have

 ¯ϕe3 = (a1+x02b1)x0+(a0+x0+12b0)(1−x0)=V3. (2.13)

From the continuity condition at and , we have

 a0+b0x0 = a1+b1x0, (2.14) a0+c0y0 = a1+c1y0 (2.15)

and the flux continuity condition along gives

 (b0,c0)⋅(y0,x0) = ρ(b1,c1)⋅(y0,x0), (2.16)

where and we have used that the normal direction of the line segment is .

Then the coefficient matrix of the above linear system for the unknowns and in this order is

 A=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝112120001−y0012(1−y20)y0012y201−x012(1−x20)0x012x200−1−x001x00−10−y010y00−y0−x00ρy0ρx0⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠ (2.17)

Tedious calculation(see appendix) shows that the determinant of the matrix is

 det(A)=14(x20+y20){ρ(x0y0−1)−x0y0}<0. (2.18)

Thus the coefficients of (2.10) are uniquely determined.

###### Remark 2.1

If and have the same value, then the piecewise linear function satisfying (2.8)-(2.9) reduces to a constant by uniqueness.

Now we can construct nodal basis functions on an interface element in general position through affine mapping. We let to denote the three-dimensional linear space spanned by these shape functions. We note that is a subspace of . Finally, we define the immersed interface finite element space as the collection of functions such that

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩  ϕ|T∈¯¯¯¯Sh(T), if T is a noninterface element,  ϕ|T∈ˆSh(T), if T is an % interface element,  ∫eϕ|T1ds=∫eϕ|T2ds, if T1,T2 are adjacent elements and e is a common edge of T1 and T2,  ∫eϕds=0, if e is part of the boundary ∂Ω.

Although for functions in the flux jump condition is enforced on line segments, they actually satisfy a weak flux jump condition along the interface. This is stated in the following lemma , whose proof is a simple application of the divergence theorem.

{lemma}

For an interface triangle T, every function satisfies the flux jump condition on in the following weak sense:

 ∫Γ∩T(β−∇ϕ−−β+∇ϕ+)⋅nΓds=0. (2.19)
{proof}

Let be any function in . By the divergence theorem, we have

 ∫Γ∩T(β−∇ϕ−−β+∇ϕ+)⋅nΓds+∫¯¯¯¯¯¯¯DE(β−∇ϕ−−β+∇ϕ+)⋅n¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny{DE}ds = ∫Trdiv(β−∇ϕ−−β+∇ϕ+)dx=0.

By the flux continuity of on ,

 ∫¯¯¯¯¯¯¯DE(β−∇ϕ−−β+∇ϕ+)⋅n¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\tiny{DE}ds=0,

which completes the proof.

### 2.2 Approximation property of nonconforming immersed interface space ˆSh(T)

In this subsection, we would like to study the approximation property of by defining an interpolation operator. The difficulty lies in the fact that does not belong to , the restriction of on , where (see Fig. 3). To overcome the difficulty, we introduce a bigger space which contains both of these spaces.

For a given interface element T, we consider a function space such that every satisfies

 ⎧⎪⎨⎪⎩p∈H1(T)∩H2(T+∩Ω+)∩H2(T−∩Ω−)∩H2(T+r)∩H2(T−r),∫Γ∩T(β−∇p−−β+∇p+)⋅nΓds=0,

where . For any , we define the following norms.

 |p|2X(T) = |p|22,T−∩Ω−+|p|22,T+∩Ω++|p|22,T−r+|p|22,T+r, ∥p∥2X(T) = ∥p∥21,T+|p|2X(T), |||[|||p]2,T = |p|X(T)+3∑i=1|¯pei|,

where are the average on each edge .

###### Remark 2.2

If , then and , where .

{lemma}

is a norm in the space which is equivalent to . {proof} Let . If , then and . Hence is linear on each of the four regions , , and . Since , is linear on each . Since satisfies flux continuity condition, . Now we apply Theorem 2.1 to conclude .

We now show the equivalence of and (cf. [3, p.77], ). First, note that by Sobolev embedding theorem, is compactly embedded in for any , where . So we see that . If , then is a continuous function on and thus

 |||[|||p]2,T ≤ C∥p∥X(T), (2.20)

where is independent of .

Now suppose that the converse

 ∥p∥X(T)≤C|||[|||p]2,T,  ∀p∈X(T)

fails for any . Then there exists a sequence in with

 ∥pk∥X(T)=1,  |||[|||pk]2,T≤1k,  k=1,2,⋯. (2.21)

Since is compactly imbedded in by Kondrasov theorem [18, p. 114], there exists a subsequence of which converges in . Without loss of generality, we can assume that the sequence itself converges. Then is a Cauchy sequence in . Noting that and , we see that is a Cauchy sequence in . By completeness, it converges to an element , and (2.20),(2.21) gives

 |||[|||p∗]2,T≤|||[|||p∗−pk]2,T+|||[|||pk]2,T≤C∥p∗−pk∥X(T)+1k→0.

But

 ∥p∗∥X(T)=1  and  |||[|||p∗]2,T=0.

This is a contradiction, since implies .

For any , we define using the average of on each edge by

 (¯¯¯¯¯¯¯¯Ihp)ei=¯pei,  i=1,2,3

and call the interpolant of in . We then define for by .

{lemma}

Let be an interface element. Then for any , we have

 ∥p−Ihp∥m,T≤Ch2−m∥p∥X(T),  m=0,1, (2.22)

where is the mesh size of . {proof} Let be a reference interface element. Then for any

 |||[|||^p−Ih^p]2,ˆT = |^p−Ih^p|X(ˆT)+3∑i=1|(¯^p−Ih^p)ei| = |^p−Ih^p|X(ˆT)=|^p|X(ˆT),

where we used the fact that on each edge and -seminorm of the piecewise linear function vanishes. Applying the scaling argument for , we have

 ∥p−Ihp∥m,T ≤ Ch1−m∥^p−Ih^p∥m,ˆT≤Ch1−m|||[|||^p−Ih^p]2,ˆT ≤ Ch1−m|^p|X(ˆT)≤Ch2−m|p|X(T).

By above lemma, Remark 2.2 and (2.6), we obtain the following interpolation estimate.

{theorem}

For any , there exists a constant such that

 ∥p−Ihp∥L2(Ω)+h∥p−Ihp∥1,h≤Ch2∥p∥˜H2(Ω), (2.23)

where .

## 3 Immersed interface FEM with ‘broken’ P1-nonconforming elements

We are now ready to define our immersed interface finite element method based on ‘broken’ -nonconforming element: Find such that

 ah(ph,ϕh) = (f,ϕh),  ∀ϕh∈ˆSh(Ω), (3.1)

where

 ah(p,ϕ) = ∑T∈Th∫Tβ∇p⋅∇ϕdx,   ∀p,ϕ∈Hh(Ω), Hh(Ω):=H10(Ω)+ˆSh(Ω).

Here, is endowed with the piecewise -norm . Note that if discrete Poincaré inequality holds, then noting that the bilinear operator is bounded and coercive on , the discrete problem (3.1) has a unique solution .

{lemma}

[Discrete Poincaré inequality] There exists a constant independent of such that for any

 C∥ϕ∥2L2(Ω)≤ah(ϕ,ϕ). (3.3)
{proof}

Let be the common edge of two adjacent elements and . Note that since , where , there exists a point such that . Then a slight modification of Lemma 2.1 in  proves the inequality.

For the energy-norm error estimate of the immersed interface finite element method, we need the well-known second Strang Lemma which is valid since is coercive.

{lemma}

[Second Strang Lemma] If are the solutions of (2.3) and (3.1) respectively, then there exists a constant such that

 ∥p−ph∥1,h≤C⎧⎨⎩infqh∈ˆSh(Ω)∥p−qh∥1,h+supϕh∈ˆSh(Ω)|ah(p,ϕh)−(f,ϕh)|∥ϕh∥1,h⎫⎬⎭. (3.4)

We shall need the following estimate; see Lemma 3 in .

{lemma}

Let be an edge of . Then there exists a constant such that for all

 ∣∣∣∫eϕ(v−¯¯¯ve)ds∣∣∣ ≤ Ch|ϕ|1,T|v|1,T,

where .

###### Remark 3.1

This lemma also holds when belongs to with understood as sum of piecewise norm .

{theorem}

Let be the solutions of (2.3) and (3.1) respectively. Then there exists a constant such that

 ∥p−ph∥1,h≤Ch∥p∥˜H2(Ω). (3.5)
{proof}

We use the second Strang Lemma. The first term is nothing but an approximation error. By Theorem 2.2, we have

 infqh∈ˆSh(Ω)∥p−qh∥1,h≤Ch∥p∥˜H2(Ω). (3.6)

For the consistency error, we have from the definition of and Green’s formula

 ah(p,ϕh)−(f,ϕh) = ∑T∈Th∫Tβ∇p⋅∇ϕhdx−∫Ωfϕhdx (3.7) = ∑T∈Th∫Tβ∇p⋅∇ϕhdx−(∑T∈Th∫Tβ∇p⋅∇ϕhdx−∑T∈Th<β∂p∂n,ϕh>∂T) = ∑T∈Th<β∂p∂n,ϕh>∂T,

where and is a unit outward normal vector on each . Since belongs to and has well-defined average value on the interior edges, and vanishing average on the boundary,we have by Lemma 3 and remark 3.1

 ∑T∈Th<β∂p∂n,ϕh>∂T = ∑T∈Th∑e⊂∂T<β∂p∂n−(¯¯¯¯¯¯¯¯¯¯¯¯β∂p∂n)e,ϕh>e (3.8) ≤ ∑T∈ThCh|β∂p∂n|1,T|ϕh|1,T ≤ Ch∥p∥˜H2(Ω)∥ϕh∥1,h.

This completes the proof.

## 4 L2-error estimate

We now apply the duality argument to obtain -norm estimate of the error. Let us consider an auxiliary problem: Given , find such that

 −div(β∇φ) = g   in Ω∖Γ, (4.1) φ = 0   on ∂Ω

with jump conditions across . Then we have

 ∥φ∥˜H2(Ω)≤C∥g∥L2(Ω). (4.2)

Let be the solution of the corresponding variational problem

 ah(vh,φh)=(vh,g),   ∀vh∈ˆSh(Ω). (4.3)

Then

 (p−ph,g) = ∑T∈Th∫Tβ∇(p