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

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.


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-bodyfitted 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 secondorder 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, [ ] ̸ = 0 and [ ] ̸ = 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 [ ] ̸ = 0 and [ ] ̸ = 0. It achieved second-order accuracy in the ∞ 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, secondorder accuracy was achieved in the ∞ 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 threedimensional problems. Another class of methods is the Boundary Condition Capturing Method [17][18][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 the ∞ 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.

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.
We seek solutions of the variable coefficient elliptic equation away from the interface Γ given by in which x = ( 1 , 2 ) 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 (x) is assumed to lie in 2 (Ω). Given functions and along the interface Γ, we prescribe the jump conditions The "±" superscripts refer to limits taken from within the subdomains Ω ± .
Finally, the boundary conditions are given by 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: where is in 1 0 (Ω). Note that although the test function is the same as in the standard finite element method, the function is not a linear combination of such basis functions. Instead, the jump conditions are enforced strongly.
Two sets of grid functions are needed and they are denoted by and 1,ℎ Cut every rectangular region [ , +1 ] × [ , +1 ] into four pieces of right triangular regions. Collecting all those triangular regions, we obtain a uniform triangulation ℎ : ⋃ ∈ ℎ ; see Figure 1. This is slightly different from the triangulation in our previous work [8] for convenience of adaptive mesh refinement.
We call a cell an interface cell if its vertices belong to different subdomains; we write = + ∪ − . + and − are   separated by a straight interface segment Γ ℎ . We call a cell a regular cell if all its vertices belong to the same subdomain, either Ω + or Ω − ; see Figure 2.
In 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 is ℎ : 1,ℎ → 1 0 (Ω). For any ℎ ∈ 1,ℎ 0 , ℎ ( ℎ ) is a standard continuous piecewise linear function, which is a linear function in every triangular cell and ℎ ( ℎ ) matches ℎ on grid points. The second one ℎ is constructed as follows. For any ℎ ∈ 1,ℎ with ℎ = ℎ at boundary points, ℎ ( ℎ ) is a piecewise linear function and matches ℎ on grid points. It is a linear function in each regular cell, just like the first extension operator ℎ ( ℎ ) = ℎ ( ℎ ) in a regular cell. In each interface cell, it consists of two pieces of linear functions; one is on + and the other is on − .
In each cell, we shall construct an approximate solution ℎ to the interface problem taking into account the jump conditions. Case 1. If is a regular cell (see Figure 2(a)), we can have the following equation: Case 2. If 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: . Then a local system defined on the interface cell can be constructed as where point 6 is the midpoint of the line segment from 4 to 5.
Solve this local system and get the values of ± 4 , ± 5 , which are denoted by the linear combinations of 1 , 2 , 3 . Then we have the following equation: Input: Triangulation set , interior point set Output: New triangulation set , new interior point set 1: Let be the length of triangulation 2: for = 1 to do 3: if 's refine sign is equal to 1 then 4: ← 1, ← 1 5: while s=1 do 6: Let ( , ) be the right angle point of 7: Let ( , ) be the middle point of the hypotenuse of 8: Let be the triangle next to and share the hypotenuse with 9: ( , ) ← 2( , ) − ( , ) 10: if ( , ) is on the boundary then 11: ← 0 (see cell 10,8,13 in Figure 3) 12: else if ( , ) ∈ then 13: ← 0 (see cell 9,6,7 in Figure 3 Putting all the cells together, we propose the following method.
Method 1. Find a discrete function ℎ ∈ 1,ℎ , such that ℎ = ℎ on boundary points so that, for all ℎ ∈ 1,ℎ 0 , we have 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 ( , ) and (2 , 2 ) grid points in both and directions, then use this numerical result to get the numerical error | ( , ) | 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.
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.

Numerical Experiments
In all numerical experiments below, the level-set function ( , ), the coefficients ± ( , ), and the solutions  is the number of interior points, and 1/2 is the number of points in one dimension. The ∞ norm and the energy norm in the whole domain Ω are measured in the following ways: Example 2. The level-set function , the coefficients ± , and the solution ± are given as follows:  Table 1 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.
The left and right of Table 2 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 | ( , ) |, 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 = 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 = 1. Figure 5(d) is a comparison of the ∞ 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 ∞ norm even when the coefficients are matrix and the ratio − / + = 1000, which is more efficient than the uniform mesh.   Example 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 ± are given as follows:  Table 3 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. 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 ± are given as follows. The interface is Lipschitz continuous and it has a kink at (0, 0), and is piecewise 2 :  Figure 7 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).

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 same ∞ 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.

Data Availability
All data can be accessed in the Numerical Experiments section of this article.

Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.