Mathematical Problems in Engineering

Volume 2018, Article ID 3721258, 12 pages

https://doi.org/10.1155/2018/3721258

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

^{1}Department of Mathematics, College of Science, China University of Petroleum (Beijing), 102249, China^{2}Department of Mathematics and Statistics and Center of Applied Physics, Louisiana Tech University, Ruston, LA 71272, USA^{3}Department of Science and Technology Teaching, China University of Political Science and Law, Beijing, China

Correspondence should be addressed to Liwei Shi; moc.liamg@ylimhiewils

Received 29 May 2018; Accepted 30 July 2018; Published 10 September 2018

Academic Editor: Nunzio Salerno

Copyright © 2018 Liqun Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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, and . 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 and . It achieved second-order accuracy in the norm for smooth interfaces and about th 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 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 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.

#### 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 . Let be an interface of co-dimension , which divides into two disjoint open subdomains, and , hence .

We seek solutions of the variable coefficient elliptic equation away from the interface given byin which denotes the spatial variables. The coefficient is assumed to be a 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 is assumed to lie in .

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

#### 3. Numerical Method

In this paper, we restrict ourselves to a rectangular domain in the plane. Given positive integers and , set and . We define a uniform Cartesian grid for and . Each is called a grid point. For the case or , a grid point is called a boundary point; otherwise it is called an interior point. The grid size is defined as

Two sets of grid functions are needed and they are denoted by and

Cut every rectangular region 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.