Numerical Method for Solving Matrix Coefficient Elliptic Equation on Irregular Domains with Sharp-Edged Boundaries

Wepresent a new second-order accurate numericalmethod for solvingmatrix coefficient elliptic equation on irregular domainswith 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 L∞ norm in both two and three dimensions and numerically very stable.

We consider the variable coefficient elliptic equation where  = ( 1 , . . .,   ) refers to the spatial variable, ∇ is the gradient operator, and the right-hand side () is assumed to lie in  2 (Ω − ).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  − () =  () on Γ. ( 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][2][3][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][6][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][17][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 kernelfree 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][27][28][29][30][31][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.

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 Ω.
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.
In both two-and three-dimensional cases, a cell Δ  belongs to one of the following two different sets: If a cell belongs to Λ 1 , we call it a regular cell; otherwise we call it an interface cell, written as Δ −  are separated by the interface segment, denoted by Γ ℎ Δ  ; see Figures 2(b) and 4.
Two extension operators are needed.The first one is ) 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  1,ℎ 0 , is a finite dimensional subspace of  1 0 (Ω).The second extension operator  ℎ 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 Δ −  .

Construction of Jump Conditions and Weak Formulation
In [26][27][28][29][30][31][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 Δ  ∈ Λ 2 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.

Advances in Numerical Analysis
Let Then For three-dimensional problem, there are two kinds of interface cells.
Upon above discussion, the system described in ( 1)-( 2) is equivalent to the following system: Since  + = 0 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  1 (Ω) is used.For  1 0 (Ω), instead of the usual inner product, we choose one which is better suited to our problem: In this way, we have the following definition.

Advances in Numerical Analysis
Method 1. Find a discrete function  ℎ ∈  1,ℎ such that  ℎ =  ℎ on boundary points and so that for all  ℎ ∈  1,ℎ 0 , 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.

Numerical Experiments
Of all the numerical experiments below, in Ω − , the levelset 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 Ω.

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: The numerical result with our method using a 40 × 40 uniform Cartesian grid is shown in Figure 5. Table 1 shows the error on different grids.
The numerical result with our method using a 40 × 40 uniform Cartesian grid is shown in Figure 6.Table 2 shows the error on different grids.
The numerical result with our method using a 40×40 uniform Cartesian grid is shown in Figure 7. Table 3 shows the error on different grids.

Numerical Examples of Three-Dimensional Problems
Figure 9 shows the interface and boundary condition using a 40 × 40 × 40 uniform Cartesian grid.Table 5 shows the error on different grids.
Figure 10 shows the interface and boundary condition using a 40 × 40 × 40 uniform Cartesian grid.Table 6 shows the error on different grids.

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.

Figure 5 :
Figure 5: Numerical result of the problem with "star" shape.

Figure 6 :
Figure 6: Numerical result of the problem with "face" shape.

Figure 8 Example 5 .
Figure 8 shows the interface and boundary condition using a 40 × 40 × 40 uniform Cartesian grid.Table4shows the error on different grids.Example 5.The level-set function , the coefficients  − , and the solution  − are given as follows:

Table 1 :
Numerical error of the problem with "star" shape.

Table 2 :
Numerical error of the problem with "face" shape.

Table 3 :
Numerical error of the problem with "target" shape.