#### Abstract

We present a new second-order accurate numerical method for solving matrix coefficient elliptic equation on irregular domains with sharp-edged boundaries. Nontraditional finite element method with non-body-fitting grids is implemented on a fictitious domain in which the irregular domains are embedded. First we set the function and coefficient in the fictitious part, and the nonsmooth boundary is then treated as an interface. The emphasis is on the construction of jump conditions on the interface; a special position for the ghost point is chosen so that the method is more accurate. The test function basis is chosen to be the standard finite element basis independent of the interface, and the solution basis is chosen to be piecewise linear satisfying the jump conditions across the interface. This is an efficient method for dealing with elliptic equations in irregular domains with non-smooth boundaries, and it is able to treat the general case of matrix coefficient. The complexity and computational expense in mesh generation is highly decreased, especially for moving boundaries, while robustness, efficiency, and accuracy are promised. Extensive numerical experiments indicate that this method is second-order accurate in the norm in both two and three dimensions and numerically very stable.

#### 1. Introduction

Let () be an open-bounded domain with a Lipschitz continuous boundary .

We consider the variable coefficient elliptic equation where refers to the spatial variable, is the gradient operator, and the right-hand side is assumed to lie in . The coefficient is a matrix that is uniformly elliptic, and its entries are continuously differentiable on .

For a given function on the boundary , the Dirichlet boundary condition is prescribed as

Elliptic partial differential equations are often used to construct models of the most basic theories underlying physics and engineering, such as electromagnetism, material science, and fluid dynamics. Different kinds of boundary conditions arise with the wide range of applications, such as Dirichlet boundary condition, Neumann boundary condition, and Robin boundary condition.

Elliptic equation on irregular domains has been studied by many researchers and several techniques have been developed. Finite element methods use a mesh triangulation to capture the boundary [1–4]. However, in many situations, such as when the boundary is moving, the mesh generation may be both computational expensive and challenging. A more preferred method is to combine the Cartesian grid method with level-set approach [5–7] to capture the boundary. Higher-order accuracy can be achieved by modifying standard difference formulas; examples are Shortley-Weller discretization [8] and the fast solution proposed by Mayo in [9] for solving Poisson or biharmonic equation on irregular regions with smooth boundary.

Since the pioneer work of Peskin [10] in 1977, much attention has been paid to the numerical solution of elliptic equations with discontinuous coefficients and singular sources on regular Cartesian grids. A large class of finite difference methods have been proposed. The main idea is to use difference scheme and stencils carefully near the interface to incorporate jump conditions and achieve high-order local truncation error using Taylor expansion. The immersed interface method [10, 11] incorporates the interface conditions into the finite difference scheme near the interface to achieve second-order accuracy based on a Taylor expansion in a local coordinate system. The boundary condition capturing method [12] uses the ghost fluid method [13] to capture the boundary conditions. The method extends the solution from one side across the interface using the jump conditions. In [14], the matched interface and boundary method was proposed to solve elliptic equations with smooth interfaces. In [15], the matched interface and boundary method was generalized to treat sharp-edged interfaces. With an elegant treatment, second-order accuracy was achieved in the norm. However, for oscillatory solutions, the errors degenerate.

The existing finite element schemes on non-body-fitted meshes are usually designed by modifying the finite element basis near the interface. The extended finite element method [16–18] was proposed by Ted Belytschko and collaborators to ease difficulties in solving problems with localized features that are not efficiently resolved by mesh refinement. In immersed finite element method [19], a Lagrangian solid mesh moves on top of a background Eulerian fluid mesh which spans the entire computational domain. The penalty finite element method [20, 21] modifies the bilinear form near the interface by penalizing the jump of the solution value (with no general flux jump) across the interface.

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 [22, 23] for two-dimensional problems and [24] for three-dimensional problems. Another recent work in this area is a class of kernel-free boundary integral (KFBI) methods for solving elliptic BVPs, presented in [25].

In this paper, we implemented the nontraditional finite element method into matrix coefficient elliptic equation on irregular domains. The idea of applying interface problem solvers into noninterface problems is completely new. The nontraditional finite element method is originally developed in [26–32] for solving elliptic or elasticity equations with sharp-edged interfaces. This method uses non-body-fitting Cartesian grids and uses different basis for the solution and test function; its linear system is independent of jump condition. In this way, it is easy to implement, especially for moving interface. It is capable of treating the more general case of variable matrix coefficient and can deliver (almost) second-order accuracy in norm. Since the method is based on the assumption of discontinuous solution, in this work, we embed the research domain with irregular boundary into a fictitious rectangular/cube. After setting the function and coefficient in the fictitious part, we construct the jump conditions and calculate the integration on the boundary of the research domain (we call it interface in the following discussion). By implementing the nontraditional finite element method on the fictitious domain containing the irregular domains which are under research, we overcome the difficulty of capturing the sharp-edged boundary with high accuracy. Both two- and three-dimensional models are studied. Extensive numerical experiments show the efficiency of this method.

The rest of the paper is organized as follows. The next section presents the discretization of two- and three-dimensional irregular domains. In Section 3, we construct the jump conditions and develop the weak formulation of the problem. In Section 4, six extensive numerical experiments are presented to show the second-order accuracy of our method. Finally some conclusions are drawn in Section 5.

#### 2. Discretization of the Domain

For simplicity, we first embed the irregular domain into a fictitious rectangular/cube ; the boundary of is denoted by , the domain between and is denoted by , and let . For ease of discussion, from now on we call the interface and the boundary of . In this paper, we restrict ourselves to a fictitious rectangular domain in two-dimensional space and a fictitious cube domain in three-dimensional space; see Figure 1. We assume that there is a Lipschitz continuous and piecewise smooth level-set function on , where , , and . A unit vector can be obtained on , which is a unit normal vector of at pointing from to .

**(a) The rectangular domain**

**(b) The cube domain**

##### 2.1. Discretization of Two-Dimensional Problems

In the two-dimensional case, is a matrix coefficient that is uniformly elliptic in . Given two positive integers and , set and . We define a uniform Cartesian grid with grid points given by for and . A grid point is called a boundary point if or ; 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

To triangulate the domain into a set of uniform triangular cells, we divide each rectangular region into two right triangular cells: one is bounded by , and and the other by , and . Collecting all these triangular cells, we obtain a uniform triangulation of ; see Figure 2(a). Note that we can also choose the hypotenuse to be and get a different uniform triangulation from the same Cartesian grid. There is no conceptual difference for our method on these two triangulations.

**(a) A schematic illustration of triangulation**

**(b) Interface cells**

If , the grid point is counted as in ; otherwise, it is counted as in .

##### 2.2. Discretization of Three-Dimensional Problems

For the three-dimensional problem, is a matrix coefficient that is uniformly elliptic in . Given positive integers , , and , set , , and .

We define a uniform Cartesian grid with grid points given by for , , and . A grid point is called a boundary point if , , or ; otherwise it is called an interior point. The grid size is defined as .

The following two sets of grid functions are needed:

Each cube cell region is cut into six tetrahedron cells. Upon collecting all these cells, we obtain a uniform territorialization of the cube cell region ; see Figure 3.

**(a) Cube cells of three-dimensional problems**

**(b) Territorialization of three-dimensional problems**

If , the grid point is counted as in ; otherwise it is counted as in .

In both two- and three-dimensional cases, a cell belongs to one of the following two different sets:

If a cell belongs to , we call it a regular cell; otherwise we call it an interface cell, written as . and are separated by the interface segment, denoted by ; see Figures 2(b) and 4.

**(a) Case 1: The interface segment is a triangle**

**(b) Case 2: The interface segment is a polygon**

Two extension operators are needed. The first one is . For any , is a standard continuous piecewise linear function, which is a linear function in every cell and matches on grid points. Clearly such a function set, denoted by , is a finite dimensional subspace of . The second extension operator is constructed as follows. For any 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 .

#### 3. Construction of Jump Conditions and Weak Formulation

In [26–32], jump conditions and are known and are used to derive the linear combination of the values of the points across the interface . However, in this paper, we need to construct and for each interface cell from the points with known value across the interface and on and then calculate the linear integral

Here is smooth on and is kept to be Lipschitz continuous. is the unit vector defined in Section 2.

For different and we have different and . For simplicity, set as 0 and on ; then . For , since , the difficult part is to generate . For accuracy consideration, we choose the position of the ghost point as the midpoint of the grid points on . Let be the value of point ; then , where is the value of the grid point / on .

For two-dimensional problem, cell is cut into two parts, and ; see Figure 2(b). If part , then is chosen to be point 1, ; if part , then is chosen to be the midpoint of points and , . Since the values of points 4 and 5 are known as and , the value of can be determined by

Let

Then

For three-dimensional problem, there are two kinds of interface cells.

*Case 1. *The interface cell is cut into a triangular pyramid and a pentahedron ; see Figure 4(a). If part , then is chosen to be point 1, ; if part , then is chosen to be the midpoint of points 2, 3, and 4, .

*Case 2. *The interface cell is cut into a tetrahedron and a pentahedron . We can always rotate the cell to make , such that is the midpoint of points 3 and 4, .

In Figure 4(a), the values of points 5, 6, and 7 are known as , , and ; therefore, the value of can be calculated by

In Figure 4(b), the values of points 5, 6, 7, and 8 are known as , , , and ; therefore, the value of can be calculated by

Let

Then

Lemma 1. *All coefficients in equations (8)–(12) are finite and independent of and .*

*Proof. *See Lemma 3.1 in [26].

Upon above discussion, the system described in (1)-(2) is equivalent to the following system:

Since is known, we only care about the value on . We generalize the weak formulation for the elliptic equation with matrix coefficient. The usual Sobolev space is used. For , instead of the usual inner product, we choose one which is better suited to our problem:

In this way, we have the following definition.

*Definition 2. *A function with and on is a weak solution of system (1)-(2) if it satisfies

A classical solution of system (1)-(2), , is necessarily a weak solution.

Theorem 3. *If , then there exists a unique weak solution of equations (1)-(2).*

*Proof. *See Theorem 2.1 in [33].

Theorem 4. *For all , can be constructed uniquely, provided and are given.*

*Proof. *See Theorem 3.1 in [26].

*Method 1. *Find a discrete function such that on boundary points and so that for all , we have

Theorem 5. *If is positive definite, then the matrix for the linear system generated by Method 1 is positive definite.*

*Proof. *For any vector ,
where and are basis functions for the solution and the test functions, respectively. According to Lemma 1, matrix is independent of ; choose such that everywhere. Let

Considering the positive definiteness of , we have

Therefore, is positive definite.

In our implementation, the integrals are computed with Gaussian quadrature rule. For each cell, we denote the midpoint of each edge as . In numerical computation, we apply the average of all the in each cell.

#### 4. Numerical Experiments

Of all the numerical experiments below, in , the level-set function , the coefficients , and the solutions are given. Hence the source term and the Dirichlet boundary condition can be directly evaluated and the jump conditions and can be constructed. The norm of the error in solution is measured over the domain . The uniform Cartesian grid is restricted to the whole domain .

##### 4.1. Numerical Examples of Two-Dimensional Problems

*Example 1. *This example has a “star” interface. The level-set function , the coefficients , and the solution are given as follows:
with , , , and

The numerical result with our method using a uniform Cartesian grid is shown in Figure 5. Table 1 shows the error on different grids.

*Example 2. *This example has a “happy face” interface. The level-set function , the coefficients , and the solution are given as follows:

The numerical result with our method using a uniform Cartesian grid is shown in Figure 6. Table 2 shows the error on different grids.

*Example 3. *This example has a “target" interface. The level-set function , the coefficients , and the solution are given as follows:

The numerical result with our method using a uniform Cartesian grid is shown in Figure 7. Table 3 shows the error on different grids.

##### 4.2. Numerical Examples of Three-Dimensional Problems

*Example 4. *The level-set function , the coefficients , and the solution are given as follows:

Figure 8 shows the interface and boundary condition using a uniform Cartesian grid. Table 4 shows the error on different grids.

*Example 5. *The level-set function , the coefficients , and the solution are given as follows:

Figure 9 shows the interface and boundary condition using a uniform Cartesian grid. Table 5 shows the error on different grids.

*Example 6. *The level-set function , the coefficients , and the solution are given as follows:

Figure 10 shows the interface and boundary condition using a uniform Cartesian grid. Table 6 shows the error on different grids.

#### 5. Conclusion

We developed a novel numerical method for solving matrix valued coefficient elliptic equation on irregular domains. The research domain with sharp-edged boundary is embedded into a fictitious rectangular/cube and ghost point is used to construct the jump conditions on the interface. Nontraditional finite element method, which employs different basis for the solution and the test function, is implemented in the whole domain. Both two- and three-dimensional problems are discussed in this paper. Extensive numerical experiments show that this method is second-order accurate in the norm in two and three dimensions. The work shows the robustness and efficiency of the method and indicates that this method can be used as an alternative way to solve elliptic equations; it lays a foundation for interface problems on irregular domains with various kinds of boundary conditions; as such, the variety of problems that are solvable is largely expanded.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The authors are grateful to the editor and referees for their valuable comments and suggestions. L. Wang’s research is supported by Science Foundation of China University of Petroleum Beijing (no. YJRC-2013-48).