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

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 Ω⁻⊂ℝ^𝑑 (𝑑=2,3) be an open-bounded domain with a Lipschitz continuous boundary Γ.

We consider the variable coefficient elliptic equation −∇⋅𝛽⁻(𝑥)∇𝑢⁻(𝑥)=𝑓(𝑥)inΩ⁻, (1) 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 𝑢⁻(𝑥)=𝑔(𝑥)onΓ. (2)

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 Ω=(𝑥_min,𝑥_max)×(𝑦_min,𝑦_max) in two-dimensional space and a fictitious cube domain Ω=(𝑥_min,𝑥_max)×(𝑦_min,𝑦_max)×(𝑧_min,𝑧_max) in three-dimensional space; see Figure 1. We assume that there is a Lipschitz continuous and piecewise smooth level-set function 𝜙 on Ω, where Γ={𝜙=0}, Ω⁻={𝜙<0}, and Ω⁺={𝜙>0}. 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

Figure 1 

The domain of Ω=𝜕Ω⋃Ω+⋃Ω−⋃Γ.

2.1. Discretization of Two-Dimensional Problems

In the two-dimensional case, 𝛽 is a 2×2 matrix coefficient that is uniformly elliptic in Ω⁻. Given two positive integers 𝐼 and 𝐽, set Δ𝑥=(𝑥_max−𝑥_min)/𝐼 and Δ𝑦=(𝑦_max−𝑦_min)/𝐽. We define a uniform Cartesian grid with grid points given by (𝑥_𝑖,𝑦_𝑗)=(𝑥_min+𝑖Δ𝑥,𝑦_min+𝑗Δ𝑦) for 𝑖=0,…,𝐼 and 𝑗=0,…,𝐽. A grid point is called a boundary point if 𝑖=0,𝐼 or 𝑗=0,𝐽; otherwise, it is called an interior point. The grid size is defined as ℎ=max(Δ𝑥,Δ𝑦)>0.

Two sets of grid functions are needed and they are denoted by 𝐻^1,ℎ=𝜔^ℎ=𝜔_𝑖,𝑗,𝐻∶0≤𝑖≤𝐼,0≤𝑗≤𝐽_01,ℎ=𝜔^ℎ=𝜔_𝑖,𝑗∈𝐻^1,ℎ∶𝜔_𝑖,𝑗=0if𝑖=0,𝐼or𝑗=0,𝐽. (3)

To triangulate the domain Ω into a set of uniform triangular cells, we divide each rectangular region [𝑥_𝑖,𝑥_𝑖+1]×[𝑦_𝑗,𝑦_𝑗+1] into two right triangular cells: one is bounded by 𝑥=𝑥_𝑖,𝑦=𝑦_𝑗, and 𝑦=((𝑦_𝑗+1−𝑦_𝑗)/(𝑥_𝑖−𝑥_𝑖+1))(𝑥−𝑥_𝑖+1)+𝑦_𝑗 and the other by 𝑥=𝑥_𝑖+1,𝑦=𝑦_𝑗+1, and 𝑦=((𝑦_𝑗+1−𝑦_𝑗)/(𝑥_𝑖−𝑥_𝑖+1))(𝑥−𝑥_𝑖+1)+𝑦_𝑗. 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 𝑦=((𝑦_𝑗+1−𝑦_𝑗)/(𝑥_𝑖+1−𝑥_𝑖))(𝑥−𝑥_𝑖)+𝑦_𝑗 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

Figure 2 

Discretization of two-dimensional problems.

If 𝜙(𝑥_𝑖,𝑦_𝑗)≤0, 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 3×3 matrix coefficient that is uniformly elliptic in Ω⁻. Given positive integers 𝐼, 𝐽, and 𝐾, set Δ𝑥=(𝑥_max−𝑥_min)/𝐼, Δ𝑦=(𝑦_max−𝑦_min)/𝐽, and Δ𝑧=(𝑧_max−𝑧_min)/𝐾.

We define a uniform Cartesian grid with grid points given by (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)=(𝑥_min+𝑖Δ𝑥,𝑦_min+𝑗Δ𝑦,𝑧_min+𝑘Δ𝑧) for 𝑖=0,…,𝐼, 𝑗=0,…,𝐽, and 𝑘=0,…,𝐾. A grid point is called a boundary point if 𝑖=0,𝐼, 𝑗=0,𝐽, or 𝑘=0,𝐾; otherwise it is called an interior point. The grid size is defined as ℎ=max(Δ𝑥,Δ𝑦,Δ𝑧)>0.

The following two sets of grid functions are needed: 𝐻^1,ℎ=𝜔^ℎ=𝜔_{𝑖,𝑗,𝑘},𝐻∶0≤𝑖≤𝐼,0≤𝑗≤𝐽,0≤𝑘≤𝐾_01,ℎ=𝜔^ℎ=𝜔_{𝑖,𝑗,𝑘}∈𝐻^1,ℎ∶𝜔_{𝑖,𝑗,𝑘}=0if𝑖=0,𝐼,or𝑗=0,𝐽,or𝑘=0,𝐾. (4)

Each cube cell region [𝑥_𝑖,𝑥_𝑖+1]×[𝑦_𝑗,𝑦_𝑗+1]×[𝑧_𝑘,𝑧_𝑘+1] 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

Figure 3 

Discretization of three-dimensional problems.

If 𝜙(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)≤0, 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: Λ₁=Δ_𝐿⊂Ω∶AllthevertexesofcellΔ_𝐿areinthesamedomainΩ⁺orΩ⁻,Λ₂=Δ_𝐿⊂Ω∶ThevertexesofcellΔ_𝐿areindifferentdomainsΩ⁺,Ω⁻. (5)(6)

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

Figure 4 

The interface segment of three-dimensional problems.

Two extension operators are needed. The first one is 𝑇^ℎ∶𝐻^1,ℎ→𝐻¹₀(Ω). For any 𝜓^ℎ∈𝐻₀^1,ℎ, 𝑇^ℎ(𝜓^ℎ) 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,ℎ, is a finite dimensional subspace of 𝐻¹₀(Ω). 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 Δ⁻_𝐿.

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 [𝑢]_Γ(𝑥)≡𝑢⁺(𝑥)−𝑢⁻(𝑥)=𝑎(𝑥),𝛽∇𝑢⋅𝑛_Γ(𝑥)≡𝑛⋅𝛽⁺(𝑥)∇𝑢⁺(𝑥)−𝑛⋅𝛽⁻(𝑥)∇𝑢⁻(𝑥)=𝑏(𝑥). (7)

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 𝑎(𝑥)=0−𝑢⁻|_Γ(𝑥)=−𝑔(𝑥). For 𝑏(𝑥), since ∇𝑢⁺=0, 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 𝑢⁻_𝐺=(1/𝑚)∑^𝑚_𝑖=1𝑢⁻_𝑖, 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 2 and 3, 𝑢⁻_𝐺=(𝑢⁻₂+𝑢⁻₃)/2. Since the values of points 4 and 5 are known as 𝑢⁻₄=𝑔(𝑥₄,𝑦₄) and 𝑢⁻₅=𝑔(𝑥₅,𝑦₅), the value of ∇𝑢⁻ can be determined by 𝑢⁻_𝑥=𝑐_1,𝑥𝑢⁻₄+𝑐_2,𝑥𝑢⁻₅+𝑐_3,𝑥𝑢⁻_𝐺,𝑢⁻_𝑦=𝑐_1,𝑦𝑢⁻₄+𝑐_2,𝑦𝑢⁻₅+𝑐_3,𝑦𝑢⁻_𝐺. (8)

Let 𝑛𝑛=¹𝑛₂,𝛽⁻=𝛽⁻₁₁𝛽⁻₁₂𝛽⁻₂₁𝛽⁻₂₂. (9)

Then 𝑏=𝛽⁺∇𝑢⁺⋅𝑛−𝛽⁻∇𝑢⁻𝛽⋅𝑛=−⁻₁₁𝑢⁻_𝑥𝑛₁+𝛽⁻₁₂𝑢⁻_𝑦𝑛₁+𝛽⁻₂₁𝑢⁻_𝑥𝑛₂+𝛽⁻₂₂𝑢⁻_𝑦𝑛₂. (10)

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, 𝑢⁻_𝐺=(𝑢⁻₂+𝑢⁻₃+𝑢⁻₄)/3.

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, 𝑢⁻_𝐺=(𝑢⁻₃+𝑢⁻₄)/2.

In Figure 4(a), the values of points 5, 6, and 7 are known as 𝑢⁻₅=𝑔(𝑥₅,𝑦₅,𝑧₅), 𝑢⁻₆=𝑔(𝑥₆,𝑦₆,𝑧₆), and 𝑢⁻₇=𝑔(𝑥₇,𝑦₇,𝑧₇); therefore, the value of ∇𝑢⁻ can be calculated by 𝑢⁻_𝑥=𝑐_1,𝑥𝑢⁻₅+𝑐_2,𝑥𝑢⁻₆+𝑐_3,𝑥𝑢⁻₇+𝑐_4,𝑥𝑢⁻_𝐺,𝑢⁻_𝑦=𝑐_1,𝑦𝑢⁻₅+𝑐_2,𝑦𝑢⁻₆+𝑐_3,𝑦𝑢⁻₇+𝑐_4,𝑦𝑢⁻_𝐺,𝑢⁻_𝑧=𝑐_1,𝑧𝑢⁻₅+𝑐_2,𝑧𝑢⁻₆+𝑐_3,𝑧𝑢⁻₇+𝑐_4,𝑧𝑢⁻_𝐺. (11)

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 𝑢⁻_𝑥=𝑐_1,𝑥𝑢⁻₅+𝑐_2,𝑥𝑢⁻₆+𝑐_3,𝑥𝑢⁻₇+𝑐_4,𝑥𝑢⁻₈+𝑐_5,𝑥𝑢⁻_𝐺,𝑢⁻_𝑦=𝑐_1,𝑦𝑢⁻₅+𝑐_2,𝑦𝑢⁻₆+𝑐_3,𝑦𝑢⁻₇+𝑐_4,𝑦𝑢⁻₈+𝑐_5,𝑦𝑢⁻_𝐺,𝑢⁻_𝑧=𝑐_1,𝑧𝑢⁻₅+𝑐_2,𝑧𝑢⁻₆+𝑐_3,𝑧𝑢⁻₇+𝑐_4,𝑧𝑢⁻₈+𝑐_5,𝑧𝑢⁻_𝐺. (12)

Let 𝑛𝑛=𝑛¹₂𝑛₃,𝛽⁻=𝛽⁻₁₁𝛽⁻₁₂𝛽⁻₁₃𝛽⁻₂₁𝛽⁻₂₂𝛽⁻₂₃𝛽⁻₃₁𝛽⁻₃₂𝛽⁻₃₃. (13)

Then 𝑏=𝛽⁺∇𝑢⁺⋅𝑛−𝛽⁻∇𝑢⁻𝛽⋅𝑛=−⁻₁₁𝑢⁻_𝑥𝑛₁+𝛽⁻₁₂𝑢⁻_𝑦𝑛₁+𝛽⁻₁₃𝑢⁻_𝑧𝑛₁+𝛽⁻₂₁𝑢⁻_𝑥𝑛₂+𝛽⁻₂₂𝑢⁻_𝑦𝑛₂+𝛽⁻₂₃𝑢⁻_𝑧𝑛₂+𝛽⁻₃₁𝑢⁻_𝑥𝑛₃+𝛽⁻₃₂𝑢⁻_𝑦𝑛₃+𝛽⁻₃₃𝑢⁻_𝑧𝑛₃. (14)

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: (𝑥)(𝑥)𝑎(𝑥)=[𝑢𝑥)−∇⋅𝛽∇𝑢=𝑓,𝑥∈Ω⧵Γ,](=−𝑢⁻(𝑥)(𝑥)𝑏(𝑥)=𝑥)=−𝑔,𝑥∈Γ,𝛽∇𝑢⋅𝑛(=−𝑛⋅𝛽⁻(𝑥)∇𝑢⁻(𝑥)𝑢(𝑥),𝑥∈Γ,=0,𝑥∈𝜕Ω. (15)

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 𝐻¹(Ω) is used. For 𝐻¹₀(Ω), instead of the usual inner product, we choose one which is better suited to our problem: _Ω−𝛽∇𝑢⋅∇𝜓+_Γ𝑏𝜓=_Ω−𝑓𝜓. (16)

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 _Ω−𝛽∇𝑢⋅∇𝜓+_Γ𝑏𝜓=_Ω−𝑓𝜓,∀𝜓∈𝐻¹₀(Ω). (17)

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 𝑢^ℎ∈𝐻^1,ℎ, 𝑈^ℎ(𝑢^ℎ) can be constructed uniquely, provided 𝑇^ℎ and 𝜙 are given.

Proof. See Theorem 3.1 in [26].

Method 1. Find a discrete function 𝑢^ℎ∈𝐻^1,ℎ such that 𝑢^ℎ=𝑔^ℎ on boundary points and so that for all 𝜓^ℎ∈𝐻₀^1,ℎ, we have _{𝐾∈𝑇ℎ𝐾−}𝛽∇𝑈^ℎ𝑢^ℎ⋅∇𝑇^ℎ𝜓^ℎ+_Γℎ𝐾𝑏𝑇^ℎ𝜓^ℎ=_{𝐾∈𝑇ℎ𝐾−}𝑓𝑇^ℎ𝜓^ℎ. (18)

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

Proof. For any vector 𝑐∈𝑅^𝑛, 𝑐^𝑇𝐴𝑐=^𝑛_𝑖,𝑗=1𝑎_𝑖𝑗𝑐_𝑖𝑐_𝑗=𝐵^𝑛_𝑖=1𝑐_𝑖𝑢^𝑖,^𝑛_𝑖=1𝑐_𝑖𝜓^𝑖, (19) 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 ^𝑛_𝑖=1𝑐_𝑖𝑢^𝑖=^𝑛_𝑖=1𝑐_𝑖𝜓^𝑖=𝑤. (20)
Considering the positive definiteness of 𝛽⁻, we have 𝑐^𝑇𝐴𝑐=𝐵[𝑤,𝑤]=_Ω−𝛽(𝑥)∇𝑤(𝑥)^𝑇∇𝑤(𝑥)𝑑𝑥>0. (21)
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: 𝜙()=−𝜃𝑟,𝜃𝑅sin_𝑡/2sin𝜃_𝑡/2+𝜃−𝜃_𝑟−(())𝜃2𝜋𝑖−1/5+𝑟,_𝑟+𝜋(2𝑖−2)5≤𝜃<𝜃_𝑟+𝜋(2𝑖−1)5,−𝑅sin𝜃_𝑡/2sin𝜃_𝑡/2−𝜃+𝜃_𝑟−(())𝜃2𝜋𝑖−1/5+𝑟,_𝑟+𝜋(2𝑖−3)5≤𝜃<𝜃_𝑟+𝜋(2𝑖−2)5, (22) with 𝜃_𝑡=𝜋/5, 𝜃_𝑟=𝜋/7, 𝑅=6/7, and 𝑖=1,2,3,4,5𝛽⁻=𝑥𝑥,𝑦²+5𝑥𝑦+1𝑥+2𝑦²,𝑢+3₋𝑥,𝑦=6+sin(6𝜋𝑥)sin6𝜋𝑦. (23)

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.

Table 1 

Numerical error of the problem with “star” shape.

Figure 5 

Numerical result of the problem with “star” shape.

Example 2. This example has a “happy face” interface. The level-set function 𝜙(𝑥,𝑦), the coefficients 𝛽⁻, and the solution 𝑢⁻ are given as follows: 𝜙𝑥,𝑦=maxmin𝜙₁,𝜙₂,𝜙₃,𝜙₄,𝜙₅,𝜙₆,min𝜙₇,𝜙₈,𝜙₁𝑥,𝑦=𝑥²+𝑦²−0.75²−0.15²,𝜙₂𝑥,𝑦=(𝑥−0.75)²+𝑦²−0.15²,𝜙₃𝑥,𝑦=(𝑥+0.75)²+𝑦²−0.15²,𝜙₄𝑥,𝑦=−0.1()0.12𝑥−0.2²−0.120.1𝑦−0.22²𝜙+0.12⋅0.1,₅𝑥,𝑦=−0.1()0.12𝑥+0.2²−0.120.1𝑦−0.22²𝜙+0.12⋅0.1,₆𝑥,𝑦=−𝑥²−𝑦+0.08²+0.12²,𝜙₇𝑥,𝑦=−𝑥²−𝑦+0.625²+0.425²,𝜙₈𝑥,𝑦=−𝑥²−𝑦+0.25²+0.2²,𝛽⁻=𝑥𝑥,𝑦²−𝑦²+3𝑥²−𝑦²𝑥+1₂−𝑦²+1𝑥²−𝑦²,𝑢+4₋𝑥,𝑦=5−5𝑥²−5𝑦². (2)

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.

Table 2 

Numerical error of the problem with “face” shape.

Figure 6 

Numerical result of the problem with “face” shape.

Example 3. This example has a “target" interface. The level-set function 𝜙, the coefficients 𝛽⁻, and the solution 𝑢⁻ are given as follows: 𝜙𝑥,𝑦=min𝜙₁,minminmin₁,min₂,minmin₃,min₄,𝜙₁𝑥,𝑦=𝑥²+𝑦²𝜙−0.02,₂𝑥,𝑦=𝑥²+𝑦²𝜙−0.08,₃𝑥,𝑦=𝑥²+𝑦²𝜙−0.15,₄𝑥,𝑦=𝑥²+𝑦²𝜙−0.25,₅𝑥,𝑦=𝑥²+𝑦²𝜙−0.35,₆𝑥,𝑦=𝑥²+𝑦²𝜙−0.45,₇𝑥,𝑦=𝑥²+𝑦²𝜙−0.55,₈𝑥,𝑦=𝑥²+𝑦²𝜙−0.80,₉𝑥,𝑦=𝑥²+𝑦²−1,min₁=−min𝜙₂,−𝜙₃,min₂=−min𝜙₄,−𝜙₅,min₃=−min𝜙₆,−𝜙₇,min₄=−min𝜙₈,−𝜙₉,𝛽⁻=0,𝑢𝑥,𝑦2+sin𝑥+𝑦02+sin𝑥+𝑦₋𝑥,𝑦=𝑥²+𝑦²+sin𝑥+𝑦. (25)

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.

Table 3 

Numerical error of the problem with “target” shape.

Figure 7 

Numerical result of the problem with “target” shape.

4.2. Numerical Examples of Three-Dimensional Problems

Example 4. The level-set function 𝜙, the coefficients 𝛽⁻, and the solution 𝑢⁻ are given as follows: 𝜙𝑥,𝑦,𝑧=maxmin₃,𝜙₆,𝜙₁𝑥,𝑦,𝑧=𝑥²+𝑦²𝜙−0.02,₂𝑥,𝑦,𝑧=𝑥²+𝑦²𝜙−0.1,₃𝑥,𝑦,𝑧=𝑥²+𝑦²𝜙−0.3,₄𝑥,𝑦,𝑧=𝑥²+𝑦²𝜙−0.6,₅𝑥,𝑦,𝑧=𝑥²+𝑦²𝜙−0.89,₆𝑥,𝑦,𝑧=𝑥²+𝑦²+𝑧²−0.9,min₁=−min𝜙₂,−𝜙₃,min₂=−min𝜙₄,−𝜙₅,min₃=min𝜙₁,minmin₁,min₂,𝛽⁻=𝑥,𝑦,𝑧4𝑥²+6sin𝑥+𝑦𝑥sin𝑥+𝑦2𝑧²(𝑥)(𝑥)(𝑧),𝑢+30.5sin𝑥0.5sincos+5₋𝑥,𝑦,𝑧=5−sin𝑥²+3𝑦²+𝑧. (26)

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

Table 4 

Numerical error of Example 4.

Figure 8 

Numerical result of Example 4.

Example 5. The level-set function 𝜙, the coefficients 𝛽⁻, and the solution 𝑢⁻ are given as follows: 𝜙𝑥,𝑦,𝑧=minmin𝜙₁,𝜙₂,min𝜙₃,𝜙₄,𝜙₁𝑥,𝑦,𝑧=−0.25+𝑥²+𝑦²+(𝑧+0.5)²,𝜙₂𝑥,𝑦,𝑧=−0.25+(𝑥−0.4)²+𝑦²+𝑧²,𝜙₃𝑥,𝑦,𝑧=−0.25+(𝑥+0.3)²+𝑦²+𝑧²,𝜙₄𝑥,𝑦,𝑧=−0.25+𝑥²+𝑦+0.5²+𝑧²,𝛽⁻=𝑥,𝑦,𝑧cos𝑥+𝑦²+3𝑧0.2sin(𝑧−𝑥)𝑧𝑧²()(𝑧)+5𝑦0.2sin𝑧−𝑥𝑦sin²,𝑢+2₋−𝑥𝑥,𝑦,𝑧=exp²+𝑦²+𝑧².20 (27)

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

Table 5 

Numerical error of Example 5.

Figure 9 

Numerical result of Example 5.

Example 6. The level-set function 𝜙, the coefficients 𝛽⁻, and the solution 𝑢⁻ are given as follows: 𝜙𝑥,𝑦,𝑧=min(𝑥−0.39)²+𝑦²+𝑧²−0.25,(𝑥+0.39)²+𝑦²+𝑧²,𝛽−0.25₋=𝑥,𝑦,𝑧4𝑥²+6sin𝑥+𝑦𝑥𝑦sin𝑥+𝑦2𝑧²(𝑥)(𝑥)+30.5sin𝑥𝑦0.5sincos𝑥𝑦+𝑧²,𝑢+5₋𝑥,𝑦,𝑧=10−2𝑥³+3𝑦²+sin𝑧−𝑦. (28)

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

Table 6 

Numerical error of Example 6.

Figure 10 

Numerical result of Example 6.

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