# A Numerical Method for Solving Elliptic Interface Problems Using Petrov-Galerkin Formulation with Adaptive Refinement

Authors: Liqun Wang; Songming Hou; Liwei Shi
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Elliptic interface problems have wide applications in engineering and science. Non-body-fitted grid has the advantage of saving the cost of mesh generation. In this paper, we propose a Petrov-Galerkin formulation using non-body-fitted grid for solving elliptic interface problems. In this method, adaptive mesh refinement is employed for cells with large errors. The new mesh still has all triangles being right triangles of the same shape. Numerical experiments show side-by-side comparison that to obtain the same accuracy, our new method has much less overall CPU time compared with the previous method even with some cost on mesh generation.


---

## Body

## 1. Introduction

Elliptic interface problems have wide applications in a variety of disciplines, such as electromagnetism, fluid dynamics, and material science. There are two types of grids used for solving such problems: the body-fitted grid and non-body-fitted grid. For the body-fitted grid, mesh generation has more computational cost than the non-body-fitted grid. The quality of the triangles is also an issue that needs special care. We focus on the discussion using non-body-fitted grid.In the past three decades, much attention has been paid to the numerical solution of elliptic interface problems on non-body-fitted regular Cartesian grids since the pioneering work of Peskin [1] on the first-order accurate immersed boundary method. Motivated by the immersed boundary method, to improve accuracy, in [2], the immersed interface method (IIM) was presented. The method achieves second-order accuracy by incorporating the interface conditions into the finite difference stencil in a way that preserves the interface conditions in both solution and its flux in the normal direction, [u]≠0 and [βun]≠0. The corresponding linear system is sparse but may not be symmetric or positive definite if there is a jump in the coefficient.Naturally, the standard finite element method has the property of symmetricity and positive definiteness. However, it used body-fitted grid. Can finite element method be applied using non-body-fitted grid? In [3], the immersed finite element (IFE) method is introduced for interface problems with homogeneous jump conditions. The idea is that the test and trial function basis are the same and are continuous but not smooth across the interface in an interface triangle. In [4], the method is generalized to deal with nonhomogeneous flux jump condition. In [5], the partially penalized IFE method is introduced to penalize the discontinuity of the IFE function on neighboring interface triangles. There are some work in the IFE formulation in three dimensions as well, such as [6, 7].As mentioned above, special treatment is required to use the IFE on non-homogeneous jump conditions. Is it possible for a finite element method to handle non-homogeneous jump conditions the same way as homogeneous jump conditions? In [8], a non-traditional finite element formulation for solving elliptic equations with smooth or sharp-edged interfaces was proposed with non-body-fitting grids for [u]≠0 and [βun]≠0. It achieved second-order accuracy in the L∞ norm for smooth interfaces and about 0.8th order for sharp-edged interfaces. In [9, 10], the method is analyzed and implemented in three dimensions. The resulting linear system is non-symmetric but positive definite.In [11], the matched interface and boundary (MIB) method was proposed to solve elliptic equations with smooth interfaces. In [12], the MIB method was generalized to treat sharp-edged interfaces. With an elegant treatment, second-order accuracy was achieved in the L∞ norm. In [13], the MIB method was extended to three-dimensional interface problems. Also, there has been a large body of work from the finite volume perspective for developing high order methods for elliptic equations in complex domains, such as [14, 15] for two-dimensional problems and [16] for three-dimensional problems. Another class of methods is the Boundary Condition Capturing Method [17–19]. In [20, 21], gradient recovery techniques were developed to improve the accuracy of gradient computation.In this paper, based on our early work on non-traditional finite element method using a non-symmetric weak formulation with uniform Cartesian grid, we improve the performance by using adaptive mesh refinement on the cells where the numerical error is large. The main idea is as follows: first, we use two different uniform Cartesian coarse grids so that we could compare their results to find at which cells the error is larger. Then, we do adaptive mesh refinement for these cells. Finally, we use our non-symmetric weak formulation on this non-uniform grid. For most problems the large errors occur only around the interface and therefore an alternative method refining around the interface is also proposed. We do present an example in which large errors occur away from the interface though. All triangles are right triangles of the same shape in the new grid. This grid is just a minor modification from the uniform Cartesian grid and the mesh generation cost is very low. However, since smaller triangles are used where the error is large, in the end theL∞ error is reduced. Numerical results show that to obtain the same L-infinity error, this new method needs much less overall CPU time compared with the previous method.


## 2. Equations and Weak Formulation

In this section, we briefly go through the equations and weak formulation in our previous work [8]. This is the foundation of our adaptive refinement method in this paper. Before adaptive refinement, the only difference between our previous work and the work in this paper is that the uniform triangular mesh is different.Consider an open bounded domainΩ⊂R2. Let Γ be an interface of co-dimension 1, which divides Ω into two disjoint open subdomains, Ω- and Ω+, hence Ω=Ω-∪Ω+∪Γ.We seek solutions of the variable coefficient elliptic equation away from the interfaceΓ given by(1)-▽·βx▽ux=fx,x∈Ω∖Γ,in which x=(x1,x2) denotes the spatial variables. The coefficient β(x) is assumed to be a 2×2 matrix that is uniformly elliptic on each disjoint subdomain, Ω- and Ω+, and its components are continuously differentiable on each disjoint subdomain, but they may be discontinuous across the interface Γ. The right-hand side f(x) is assumed to lie in L2(Ω).Given functionsa and b along the interface Γ, we prescribe the jump conditions(2)uΓx≡u+x-u-x=axβ▽u·nΓx≡n·β+x▽u+x-n·β-x▽u-x=bx.The “±” superscripts refer to limits taken from within the subdomains Ω±.Finally, the boundary conditions are given by(3)ux=gx,x∈∂Ω.The jump conditions are enforced strongly in the local system. Also, the flux jump condition appears in the weak formulation.We introduce the weak solution by the standard procedure of multiplying (1) by a test function ψ and integration by parts:(4)∫Ω+β▽u·▽ψ+∫Ω-β▽u·▽ψ=∫Ωfψ-∫Γbψ,where ψ is in H01(Ω). Note that although the test function ψ is the same as in the standard finite element method, the function u is not a linear combination of such basis functions. Instead, the jump conditions are enforced strongly.


## 3. Numerical Method

In this paper, we restrict ourselves to a rectangular domainΩ=(xmin,xmax)×(ymin,ymax) in the plane. Given positive integers I and J, set Δx=(xmax-xmin)/I and Δy=(ymax-ymin)/J. We define a uniform Cartesian grid (xi,yj)=(xmin+iΔx,ymin+jΔy) for i=0,…,I and j=0,…,J. Each (xi,yj) is called a grid point. For the case i=0,I or j=0,J, a grid point is called a boundary point; otherwise it is called an interior point. The grid size is defined as h=max⁡Δx,Δy>0.Two sets of grid functions are needed and they are denoted by(5)H1,h=ωh=ωi,j:0≤i≤I,0≤j≤Jand(6)H01,h=ωh=ωi,j∈H1,h:ωi,j=0ifi=0,Iorj=0,JCut every rectangular region[xi,xi+1]×[yj,yj+1] into four pieces of right triangular regions. Collecting all those triangular regions, we obtain a uniform triangulation Th:⋃K∈Th‍K; see Figure 1. This is slightly different from the triangulation in our previous work [8] for convenience of adaptive mesh refinement.Figure 1
A uniform triangulation.We call a cellK an interface cell if its vertices belong to different subdomains; we write K=K+∪K-. K+ and K- are separated by a straight interface segment ΓKh. We call a cell K a regular cell if all its vertices belong to the same subdomain, either Ω+ or Ω-; see Figure 2.Figure 2
Cells.
(a)
Regular cell (b)
Interface cellIn order to introduce our non-traditional finite element method, two finite element isomorphism mappings are needed from coefficient vector to finite element functions. The first one isTh:H1,h→H01(Ω). For any ψh∈H01,h, Th(ψh) is a standard continuous piecewise linear function, which is a linear function in every triangular cell and Th(ψh) matches ψh on grid points. The second one Uh is constructed as follows. For any uh∈H1,h with uh=gh at boundary points, Uh(uh) is a piecewise linear function and matches uh on grid points. It is a linear function in each regular cell, just like the first extension operator Uh(uh)=Th(uh) in a regular cell. In each interface cell, it consists of two pieces of linear functions; one is on K+ and the other is on K-.In each cell, we shall construct an approximate solutionuh to the interface problem taking into account the jump conditions.Case 1.
IfK is a regular cell (see Figure 2(a)), we can have the following equation:(7)∫Kβ∇Uhuh·∇Thψh=∫KfThψh.Case 2.
IfK is an interface cell (see Figure 2(b)), notice that points 2,3,4,5 are coplanar, and the value of point 5 can be denoted as a linear combination of the values of points 2,3,4: u5=c1u2+c2u3+c3u4. Then a local system defined on the interface cell K can be constructed as(8)u4=a4,u5=a5,c1u2+c2u3+c3u4=u5,β▽u·n6=b6,where point 6 is the midpoint of the line segment from 4 to 5.Solve this local system and get the values ofu4±,u5±, which are denoted by the linear combinations of u1,u2,u3. Then we have the following equation:(9)∫K+β∇Uhuh·∇Thψh+∫K-β∇Uhuh·∇Thψh=∫K+fThψh+∫K-fThψh-∫ΓKhbThψh.Putting all the cells together, we propose the following method.Method 1.
Find a discrete functionuh∈H1,h, such that uh=gh on boundary points so that, for all ψh∈H01,h, we have(10)∑K∈Th∫K+β∇Uhuh·∇Thψh+∫K-β∇Uhuh·∇Thψh=∑K∈Th∫K+fThψh+∫K-fThψh-∫ΓKhbThψh.
In order to improve the efficiency of our method, we introduce the adaptive refinement method to get a refined mesh. The idea is to use the error scale of the coarse grid to get the refinement triangulation set. First calculate the numerical result of the uniform coarse grids using(nx,ny) and (2nx,2ny) grid points in both x and y directions, then use this numerical result to get the numerical error enx,ny on the corresponding points. Based on the error scale on different grid points, we can decide the refinement scale of each triangulation and then get the refined mesh with our adaptive refinement method; we call this mesh adaptive mesh. The detailed algorithm is given in Algorithm 1.Algorithm 1:Adaptive mesh refinement.
Input: Triangulation set T, interior point set P Output: New triangulation set T, new interior point set P 1: Let tn be the length of triangulation T 2: for  i=1 to tn  do 3:if  Ti’s refine sign is equal to 1  then 4:s←1, it←1 5:while  s=1  do 6:Let (x,y) be the right angle point of Ti 7:Let (xm,ym) be the middle point of the hypotenuse of Ti 8:Let Tj be the triangle next to Ti and share the hypotenuse with Ti 9:(xn,yn)←2(xm,ym)-(x,y) 10:if  (xm,ym) is on the boundary  then 11:s←0 (see cell K10,8,13 in Figure 3) 12:else if  (xn,yn)∈P  then 13:s←0 (see cell K9,6,7 in Figure 3) 14:else 15:s←1, it←it+1, Ti←Tj, add Ti into a triangle set NT (see cell K15,4,7 in Figure 3) 16:end if 17:end while 18:for  j=it down to 1  do 19:Ti←NTj 20:Let (x,y) be the right angle point of Ti 21:Let (xm,ym) be the middle point of the hypotenuse of Ti 22:(xn,yn)←2(xm,ym)-(x,y) 23:if  (xm,ym) is on the boundary  then 24:Connect (x,y) and (xm,ym), separate the original cell into two new cells, delete the original cell from T and add two new cells into T 25:else if  (xn,yn)∈P  then 26:Connect (x,y) and (xm,ym), separate the two cells into four new cells, delete the two cells from T, add four new cells into T, and add the new point (xm,ym) into P 27:end if 28:end for 29:end if 30: end forFigure 3
Refined mesh.Remark. When it comes to numerical computation of interface problems, the major part of the error is from the region around the interface, so the refined mesh around this region needs to be the smallest triangulation. This is why we propose the adaptive mesh with further refinement on the interface, called adaptive interface mesh.


## 4. Numerical Experiments

In all numerical experiments below, the level-set functionϕ(x,y), the coefficients β±(x,y), and the solutions (11)u=u+x,y,inΩ+u-x,y,inΩ-are given. Hence (12)f=-∇·β∇u,a=u+-u-,b=β+∇u+·n-β-∇u-·non the whole domain Ω. g is obtained as a proper Dirichlet boundary condition, since the solutions are given.In Examples2, 4, and 5, we solve the problem in the rectangular domain Ω=[-1,1]×[-1,1]. In Example 3, we solve the problem in the rectangular domain Ω=[0,1]×[0,1]. N is the number of interior points, and N1/2 is the number of points in one dimension. The L∞ norm and the energy norm in the whole domain Ω are measured in the following ways: (13)eu∞=maxx∈Ω⁡u-uh,euβ=∫Ω▽u-uhTβ▽u-uh=∫Ω+▽u-uhTβ▽u-uh+∫Ω-▽u-uhTβ▽u-uh.Example 2.
The level-set functionϕ, the coefficients β±, and the solution u± are given as follows: (14)
ϕ
x
,
y
=
x
2
+
y
2
-
0.15
,
β
+
x
,
y
=
1000
x
2
+
3
cos
⁡
x
cos
⁡
x
y
2
+
2
,
β
-
x
,
y
=
x
2
+
3
c
o
s
x
c
o
s
x
y
2
+
2
,
u
+
x
,
y
=
s
i
n
5
x
y
+
1
,
u
-
x
,
y
=
s
i
n
10
x
y
.Table1 shows the error and CPU time on uniform triangular mesh, adaptive mesh, and adaptive interface mesh. From this table we can see that the adaptive mesh method has a better result than the uniform mesh method, and the method with adaptive interface mesh is significantly faster and more accurate than the other two methods. The left and right of Figure 4 are a comparison of the adaptive mesh, the numerical solution, and the numerical error using the adaptive mesh method with 3663 interior points (left) and the adaptive interface mesh method with 4167 interior points (right). The error plots do not include the outside boundary points where the solution is given with no error. From Figure 4(e) we can see that the maximum error comes from the area around the interface. Therefore we refine this place with the smallest triangulation. Figure 4(f) is the numerical error after the refinement; it is clear in this figure that the error around the interface has been decreased to the same scale as the error away from the interface.Table 1
Numerical results of Example2 with different mesh.
Mesh N e u ∞ CPU Time Uniform Mesh 404101 1.35e-3 6392.5s Adaptive Mesh 147664 1.32e-3 2218.3s Adaptive Interface Mesh 18249 1.28e-3 120.8sFigure 4
Numerical results of Example2.Example 3.
The level-set functionϕ, the coefficients β±, and the solution u± are given as follows: (15)ϕx,y=x+0.5,β+x,y=x2+100y2+1,β-x,y=1000x2+100y2+1,u+x,y=x4+x3-3x2+2y,x>0.5u+x,y=-74x-0.5+2y-916,x≤0.5u-x,y=2x+3y+3.The left and right of Table2 show the error on uniform triangular meshes and adaptive meshes with different grids, respectively. From this table we can see that the method with adaptive mesh is much faster than the method with uniform mesh. We plot the adaptive mesh using the adaptive mesh method with 5004 interior points in Figure 5(a). The numerical solution and numerical error are shown in (b) and (c) of Figure 5. After calculating the coarse grid error |e(nx,ny)|, we noticed that the error around the interface is smaller than other areas away from the interface, so in this example we do not need to use the adaptive interface method to get a further refinement mesh around the interface, as we did in Example 2. In Figure 5(c), it is clear that the maximum error lies on the line x=1, and from Figure 5(a) we can see that our method refines the grid around this line automatically and the smallest triangulation lies exactly on x=1. Figure 5(d) is a comparison of the L∞ error on uniform meshes and adaptive meshes with different grid, which shows that the convergency rate of the adaptive mesh is higher than the uniform mesh and that the adaptive mesh can get to higher than second-order accuracy in L∞ norm even when the coefficients are matrix and the ratio β-/β+=1000, which is more efficient than the uniform mesh.Table 2
Numerical results of Example3 with different mesh.
Uniform Mesh Adaptive Mesh N e u ∞ CPU Time N e u ∞ CPU Time 12641 3.02e-4 27.0s 4109 3.02e-4 21.3s 30505 1.27e-4 84.4s 8377 1.30e-4 37.4s 50881 7.69e-5 179.5s 16171 7.69e-5 73.0s 120541 3.27e-5 723.8s 30096 3.27e-5 155.7s 204161 1.94e-5 1789.4s 57391 1.94e-5 436.8s 479221 8.29e-6 18119.7s 106387 8.22e-6 1206.8s  Order: 1.98 Order: 2.17Figure 5
Numerical results of adaptive mesh for Example3.
(a)
Grids (b)
Numerical solution (c)
Numerical error (d)
L ∞ error for different gridsExample 4.
This example is taken from [22]. We use this example to compare our method with the adaptive immersed finite element method. The level-set function ϕ, the coefficients β±, and the solution u± are given as follows: (16)ϕx,y=x2+y2-0.25,β+x,y=1000,β-x,y=1,rx,y=x2+y2,u+x,y=r31000+1-110000.53,u-x,y=r3.Table3 is a comparison of the numerical results and CPU time of the adaptive IFEM mesh in [22], the adaptive mesh, and adaptive interface mesh in this paper. From this table we can see that our adaptive method with 933 interior points gives a better result than the adaptive IFEM mesh with 4021 interior points, and the result of our adaptive interface mesh is better than the adaptive mesh. The left and right of Figure 6 are a comparison of the adaptive mesh, the numerical solution, and the numerical error using the adaptive mesh method with 585 interior points (left) and the adaptive interface mesh method with 825 interior points (right). From this figure we can see that after the refinement using the adaptive method, the maximum error comes from area inside the interface but not exactly lies on the interface, so the difference of the results between the adaptive interface mesh and the adaptive mesh is not as large as it is in Example 2.Table 3
Numerical results of Example4 with different mesh.
Mesh N e u β e u ∞ CPU Time Adaptive IFEM Mesh in [22] 4021 4.11e-2  Adaptive Mesh 933 4.41e-2 1.32e-3 15.3s 3417 1.32e-2 3.75e-4 22.0s  Adaptive Interface Mesh 2017 1.28e-2 3.84e-4 17.6s 3453 8.86e-3 1.65e-4 26.0sFigure 6
Numerical results of Example4.Example 5.
This example is taken from [8]. It is an interface problem with a singular point. The level-set function ϕ, the coefficients β±, and the solution u± are given as follows. The interface is Lipschitz continuous and it has a kink at (0,0), and u is piecewise H2: (17)ϕx,y=y-2x,x+y>0ϕx,y=y+x2,x+y≤0β+x,y=1,β-x,y=2+sin⁡x+y,u+x,y=8,u-x,y=x2+y25/6+sin⁡x+yFigure7 shows the jump of solution and flux in the normal direction on the interface. Table 4 is a comparison of the numerical results and CPU time of the uniform mesh, adaptive mesh, and adaptive interface mesh with different grids. All three methods can get to higher than 1.5th order accuracy. The result of the method using the adaptive mesh is obviously better than the method using the uniform mesh, and the result of the method using the adaptive interface mesh is slightly better result than the method using the adaptive mesh, and the difference of the results of three meshes gets clearer as the triangulation becomes smaller. The left and right of Figure 8 are a comparison of the adaptive mesh, the numerical solution, and the numerical error using the adaptive mesh method with 2385 interior points (left) and the adaptive interface mesh method with 3037 interior points (right). In this example we can see that after the refinement using the adaptive method, the maximum error lies around the interface, so we further refined the mesh around the interface and get a better result as shown in Figure 8(f).Table 4
Numerical results of Example5 with different mesh.
Mesh N e u ∞ CPU Time Order Uniform Mesh 1105 6.65e-3 2.6s 1.58 8065 1.59e-3 18.5s 28561 5.61e-4 86.5s 99905 1.88e-4 561.3s 400513 6.59e-5 6448.1s  Adaptive Mesh 567 6.13e-3 3.6s 1.62 2385 1.56e-3 14.3s 9711 5.39e-4 66.2s 39038 1.82e-4 346.6s 155869 6.17e-5 3329.8s  Adaptive Interface Mesh 874 5.76e-3 5.8s 1.86 3037 1.42e-3 18.4s 11035 4.63e-4 72.0s 41714 1.35e-4 390.1s 161208 4.25e-5 3343.4sFigure 7
Jumps on the interface of Example5.
(a)
Jump of solution on the interface (b)
Jump of flux in the normal direction on the interfaceFigure 8
Numerical results of Example5.


## 5. Conclusion

In this paper, we proposed the adaptive refinement method to solve the elliptic interface problems using a non-symmetric weak formulation. The key idea is to use two uniform coarse grids to find the location of the maximum error. Sometimes this step can be omitted in practice because the cells with maximum error often locates near the interface, but we provided examples in which the cells with maximum errors are away from the interface. After fixing the location of the maximum error, we use the adaptive mesh refinement method to reduce the maximum error, so that the error of the whole area is of the same scale. Numerical experiments show that the adaptive refinement method has remarkably higher convergency and accuracy than the method using uniform Cartesian grid, and for the sameL∞ error the overall CPU time is highly reduced using our new method. The study of the three-dimensional form of this new method will be our future work, so that we can solve problems that are more practical.



---

*Source: 3721258-2018-09-10.xml*