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

Authors: Liqun Wang; Liwei Shi
Journal: Advances in Numerical Analysis
 (2013)

---

## 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 theL∞ norm in both two and three dimensions and numerically very stable.


---

## Body

## 1. Introduction

LetΩ-⊂ℝd (d=2,3) be an open-bounded domain with a Lipschitz continuous boundary Γ.We consider the variable coefficient elliptic equation(1)-∇·(β-(x)∇u-(x))=f(x)inΩ-,
where x=(x1,…,xd) refers to the spatial variable, ∇ is the gradient operator, and the right-hand side f(x) is assumed to lie in L2(Ω-). The coefficient β-(x) is a d×d matrix that is uniformly elliptic, and its entries are continuously differentiable on Ω-.For a given functiong(x) on the boundary Γ, the Dirichlet boundary condition is prescribed as
(2)u-(x)=g(x)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–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 L∞ 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 L∞ 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 Section3, 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 Ω=(xmin,xmax)×(ymin,ymax) in two-dimensional space and a fictitious cube domain Ω=(xmin,xmax)×(ymin,ymax)×(zmin,zmax) 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 n=∇ϕ/|∇ϕ| can be obtained on Ω-, which is a unit normal vector of Γ at Γ pointing from Ω- to Ω+.The domain ofΩ=∂Ω⋃Ω+⋃Ω-⋃Γ.
(a)
The rectangular domain
(b)
The cube domain


### 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 I and J, set Δx=(xmax-xmin)/I and Δy=(ymax-ymin)/J. We define a uniform Cartesian grid with grid points given by (xi,yj)=(xmin+iΔx,ymin+jΔy) for i=0,…,I and j=0,…,J. A grid point is called a boundary point if i=0,I or j=0,J; otherwise, it is called an interior point. The grid size is defined as h=max(Δx,Δy)>0.Two sets of grid functions are needed and they are denoted by(3)H1,h={ωh=(ωi,j):0≤i≤I,0≤j≤J},H01,h={ωh=(ωi,j)∈H1,h:ωi,j=0ifi=0,Iorj=0,J}.To triangulate the domainΩ into a set of uniform triangular cells, we divide each rectangular region [xi,xi+1]×[yj,yj+1] into two right triangular cells: one is bounded by x=xi,y=yj, and y=((yj+1-yj)/(xi-xi+1))(x-xi+1)+yj and the other by x=xi+1,y=yj+1, and y=((yj+1-yj)/(xi-xi+1))(x-xi+1)+yj. 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 y=((yj+1-yj)/(xi+1-xi))(x-xi)+yj and get a different uniform triangulation from the same Cartesian grid. There is no conceptual difference for our method on these two triangulations.Discretization of two-dimensional problems.
(a)
A schematic illustration of triangulation
(b)
Interface cellsIfϕ(xi,yj)≤0, the grid point (xi,yj) 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 I, J, and K, set Δx=(xmax-xmin)/I, Δy=(ymax-ymin)/J, and Δz=(zmax-zmin)/K.We define a uniform Cartesian grid with grid points given by(xi,yj,zk)=(xmin+iΔx,ymin+jΔy,zmin+kΔz) for i=0,…,I, j=0,…,J, and k=0,…,K. A grid point is called a boundary point if i=0,I, j=0,J, or k=0,K; otherwise it is called an interior point. The grid size is defined as h=max(Δx,Δy,Δz)>0.The following two sets of grid functions are needed:(4)H1,h={ωh=(ωi,j,k):0≤i≤I,0≤j≤J,0≤k≤K},H01,h={ωh=(ωi,j,k)∈H1,h:ωi,j,k=0ifi=0,I,orj=0,J,ork=0,K}.Each cube cell region[xi,xi+1]×[yj,yj+1]×[zk,zk+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.Discretization of three-dimensional problems.
(a)
Cube cells of three-dimensional problems
(b)
Territorialization of three-dimensional problemsIfϕ(xi,yj,zk)≤0, the grid point (xi,yj,zk) is counted as in Ω-¯; otherwise it is counted as in Ω+.In both two- and three-dimensional cases, a cellΔL belongs to one of the following two different sets:
(5)Λ1={ΔL⊂Ω:AllthevertexesofcellΔLareinthesamedomain(Ω+orΩ-)},(6)Λ2={ΔL⊂Ω:ThevertexesofcellΔLareindifferentdomains(Ω+,Ω-)}.If a cell belongs toΛ1, we call it a regular cell; otherwise we call it an interface cell, written as ΔL=ΔL+⋃‍ΔL-. ΔL+ and ΔL- are separated by the interface segment, denoted by ΓΔLh; see Figures 2(b) and 4.The interface segment of three-dimensional problems.
(a)
Case 1: The interface segment is a triangle
(b)
Case 2: The interface segment is a polygonTwo extension operators are needed. The first one isTh:H1,h→H01(Ω). For any ψh∈H01,h, Th(ψh) is a standard continuous piecewise linear function, which is a linear function in every cell and Th(ψh) matches ψh on grid points. Clearly such a function set, denoted by H01,h, is a finite dimensional subspace of H01(Ω). The second extension operator Uh is constructed as follows. For any uh∈H1,h with uh=gh at boundary points, Uh(uh) is a piecewise linear function and matches uh on grid points. It is a linear function in each regular cell, just like the first extension operator Uh(uh)=Th(uh) in a regular cell. In each interface cell, it consists of two pieces of linear functions: one is on ΔL+ and the other is on ΔL-.


## 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 I and J, set Δx=(xmax-xmin)/I and Δy=(ymax-ymin)/J. We define a uniform Cartesian grid with grid points given by (xi,yj)=(xmin+iΔx,ymin+jΔy) for i=0,…,I and j=0,…,J. A grid point is called a boundary point if i=0,I or j=0,J; otherwise, it is called an interior point. The grid size is defined as h=max(Δx,Δy)>0.Two sets of grid functions are needed and they are denoted by(3)H1,h={ωh=(ωi,j):0≤i≤I,0≤j≤J},H01,h={ωh=(ωi,j)∈H1,h:ωi,j=0ifi=0,Iorj=0,J}.To triangulate the domainΩ into a set of uniform triangular cells, we divide each rectangular region [xi,xi+1]×[yj,yj+1] into two right triangular cells: one is bounded by x=xi,y=yj, and y=((yj+1-yj)/(xi-xi+1))(x-xi+1)+yj and the other by x=xi+1,y=yj+1, and y=((yj+1-yj)/(xi-xi+1))(x-xi+1)+yj. 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 y=((yj+1-yj)/(xi+1-xi))(x-xi)+yj and get a different uniform triangulation from the same Cartesian grid. There is no conceptual difference for our method on these two triangulations.Discretization of two-dimensional problems.
(a)
A schematic illustration of triangulation
(b)
Interface cellsIfϕ(xi,yj)≤0, the grid point (xi,yj) 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 I, J, and K, set Δx=(xmax-xmin)/I, Δy=(ymax-ymin)/J, and Δz=(zmax-zmin)/K.We define a uniform Cartesian grid with grid points given by(xi,yj,zk)=(xmin+iΔx,ymin+jΔy,zmin+kΔz) for i=0,…,I, j=0,…,J, and k=0,…,K. A grid point is called a boundary point if i=0,I, j=0,J, or k=0,K; otherwise it is called an interior point. The grid size is defined as h=max(Δx,Δy,Δz)>0.The following two sets of grid functions are needed:(4)H1,h={ωh=(ωi,j,k):0≤i≤I,0≤j≤J,0≤k≤K},H01,h={ωh=(ωi,j,k)∈H1,h:ωi,j,k=0ifi=0,I,orj=0,J,ork=0,K}.Each cube cell region[xi,xi+1]×[yj,yj+1]×[zk,zk+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.Discretization of three-dimensional problems.
(a)
Cube cells of three-dimensional problems
(b)
Territorialization of three-dimensional problemsIfϕ(xi,yj,zk)≤0, the grid point (xi,yj,zk) is counted as in Ω-¯; otherwise it is counted as in Ω+.In both two- and three-dimensional cases, a cellΔL belongs to one of the following two different sets:
(5)Λ1={ΔL⊂Ω:AllthevertexesofcellΔLareinthesamedomain(Ω+orΩ-)},(6)Λ2={ΔL⊂Ω:ThevertexesofcellΔLareindifferentdomains(Ω+,Ω-)}.If a cell belongs toΛ1, we call it a regular cell; otherwise we call it an interface cell, written as ΔL=ΔL+⋃‍ΔL-. ΔL+ and ΔL- are separated by the interface segment, denoted by ΓΔLh; see Figures 2(b) and 4.The interface segment of three-dimensional problems.
(a)
Case 1: The interface segment is a triangle
(b)
Case 2: The interface segment is a polygonTwo extension operators are needed. The first one isTh:H1,h→H01(Ω). For any ψh∈H01,h, Th(ψh) is a standard continuous piecewise linear function, which is a linear function in every cell and Th(ψh) matches ψh on grid points. Clearly such a function set, denoted by H01,h, is a finite dimensional subspace of H01(Ω). The second extension operator Uh is constructed as follows. For any uh∈H1,h with uh=gh at boundary points, Uh(uh) is a piecewise linear function and matches uh on grid points. It is a linear function in each regular cell, just like the first extension operator Uh(uh)=Th(uh) in a regular cell. In each interface cell, it consists of two pieces of linear functions: one is on ΔL+ and the other is on ΔL-.


## 3. Construction of Jump Conditions and Weak Formulation

In [26–32], jump conditions a and b 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 a and b for each interface cell ΔL∈Λ2 from the points with known value across the interface Γ and u- on Ω- and then calculate the linear integral
(7)[u]Γ(x)≡u+(x)-u-(x)=a(x),[(β∇u)·n]Γ(x)≡n·(β+(x)∇u+(x))-n·(β-(x)∇u-(x))=b(x).Hereβ is smooth on Ω and Γ is kept to be Lipschitz continuous. n is the unit vector defined in Section 2.For differentu+ and β+ we have different a and b. For simplicity, set u+ as 0 and β+=β- on Ω+; then a(x)=0-u-|Γ(x)=-g(x). For b(x), since ∇u+=0, the difficult part is to generate ∇u-. For accuracy consideration, we choose the position of the ghost point G=(xG,yG)∈ΔL as the midpoint of the grid points on ΔL⋂‍Ω-. Let uG be the value of point G; then uG-=(1/m)∑i=1m‍ui-, where ui- is the value of the grid point (xi,yi)/(xi,yi,zi) on ΔL⋂‍Ω-.For two-dimensional problem, cellΔ123 is cut into two parts, P145 and P4235; see Figure 2(b). If part P145∈Ω-, then G is chosen to be point 1, uG-=u1-; if part P4235∈Ω-, then G is chosen to be the midpoint of points 2 and 3, uG-=(u2-+u3-)/2. Since the values of points 4 and 5 are known as u4-=g(x4,y4) and u5-=g(x5,y5), the value of ∇u- can be determined by
(8)ux-=c1,xu4-+c2,xu5-+c3,xuG-,uy-=c1,yu4-+c2,yu5-+c3,yuG-.Let(9)n=(n1n2),β-=(β11-β12-β21-β22-).Then(10)b=β+∇u+·n-β-∇u-·n=-(β11-ux-n1+β12-uy-n1=-b+β21-ux-n2+β22-uy-n2).For three-dimensional problem, there are two kinds of interface cells.Case 1.
The interface cellΔ1234 is cut into a triangular pyramid P1567 and a pentahedron P567234; see Figure 4(a). If part P1567∈Ω-, then G is chosen to be point 1, uG-=u1-; if part P567234∈Ω-, then G is chosen to be the midpoint of points 2, 3, and 4, uG-=(u2-+u3-+u4-)/3.Case 2.
The interface cellΔ1234 is cut into a tetrahedron P125678 and a pentahedron P567834. We can always rotate the cell to make P567834∈Ω-, such that G is the midpoint of points 3 and 4, uG-=(u3-+u4-)/2.In Figure4(a), the values of points 5, 6, and 7 are known as u5-=g(x5,y5,z5), u6-=g(x6,y6,z6), and u7-=g(x7,y7,z7); therefore, the value of ∇u- can be calculated by
(11)ux-=c1,xu5-+c2,xu6-+c3,xu7-+c4,xuG-,uy-=c1,yu5-+c2,yu6-+c3,yu7-+c4,yuG-,uz-=c1,zu5-+c2,zu6-+c3,zu7-+c4,zuG-.In Figure4(b), the values of points 5, 6, 7, and 8 are known as u5-=g(x5,y5,z5), u6-=g(x6,y6,z6), u7-=g(x7,y7,z7), and u8-=g(x8,y8,z8); therefore, the value of ∇u- can be calculated by
(12)ux-=c1,xu5-+c2,xu6-+c3,xu7-+c4,xu8-+c5,xuG-,uy-=c1,yu5-+c2,yu6-+c3,yu7-+c4,yu8-+c5,yuG-,uz-=c1,zu5-+c2,zu6-+c3,zu7-+c4,zu8-+c5,zuG-.Let(13)n=(n1n2n3),β-=(β11-β12-β13-β21-β22-β23-β31-β32-β33-).Then(14)b=β+∇u+·n-β-∇u-·n=-(β11-ux-n1+β12-uy-n1+β13-uz-n1+β21-ux-n2+β22-uy-n2=-m+β23-uz-n2+β31-ux-n3+β32-uy-n3+β33-uz-n3).Lemma 1.
All coefficientsc in equations (8)–(12) are finite and independent of uh and a.Proof.
See Lemma 3.1 in [26].Upon above discussion, the system described in (1)-(2) is equivalent to the following system:
(15)-∇·(β∇u(x))=f(x),x∈Ω∖Γ,a(x)=[u](x)=-u-(x)=-g(x),x∈Γ,b(x)=[(β∇u)·n](x)=-n·β-(x)∇u-(x),x∈Γ,u(x)=0,x∈∂Ω.Sinceu+=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 H1(Ω) is used. For H01(Ω), instead of the usual inner product, we choose one which is better suited to our problem:
(16)∫Ω-‍β∇u·∇ψ+∫Γ‍bψ=∫Ω-‍fψ.In this way, we have the following definition.Definition 2.
A functionu with u|Ω-∈H1(Ω-) and u=g on Γ is a weak solution of system (1)-(2) if it satisfies
(17)∫Ω-‍β∇u·∇ψ+∫Γ‍bψ=∫Ω-‍fψ,∀ψ∈H01(Ω).A classical solution of system (1)-(2), u|Ω-∈C2(Ω-¯), is necessarily a weak solution.Theorem 3.
Iff∈L2(Ω), then there exists a unique weak solution of equations (1)-(2).Proof.
See Theorem 2.1 in [33].Theorem 4.
For alluh∈H1,h, Uh(uh) can be constructed uniquely, provided Th and ϕ are given.Proof.
See Theorem 3.1 in [26].Method 1.
Find a discrete functionuh∈H1,h such that uh=gh on boundary points and so that for all ψh∈H01,h, we have
(18)∑K∈Th‍(∫K-‍β∇Uh(uh)·∇Th(ψh)+∫ΓKh‍bTh(ψh))=∑K∈Th‍∫K-‍fTh(ψh).Theorem 5.
Ifβ- is positive definite, then the matrix A for the linear system generated by Method 1 is positive definite.Proof.
For any vectorc∈Rn,
(19)cTAc=∑i,j=1n‍aijcicj=B[∑i=1n‍ciui,∑i=1n‍ciψi],
where ui and ψi are basis functions for the solution and the test functions, respectively. According to Lemma 1, matrix A is independent of a; choose a such that ui=ψi everywhere. Let
(20)∑i=1n‍ciui=∑i=1n‍ciψi=w.
Considering the positive definiteness ofβ-, we have
(21)cTAc=B[w,w]=∫Ω-‍(β(x)∇w(x))T∇w(x)dx>0.
Therefore,A is positive definite.In our implementation, the integrals are computed with Gaussian quadrature rule. For each cell, we denote the midpoint of each edgeuiuj as pij. In numerical computation, we apply the average of all the f(pij) in each cell.


## 4. Numerical Experiments

Of all the numerical experiments below, inΩ-¯, the level-set function ϕ, the coefficients β-, and the solutions u- are given. Hence the source term f and the Dirichlet boundary condition g can be directly evaluated and the jump conditions a and b can be constructed. The L∞ 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 u- are given as follows:
(22)ϕ(r,θ)={-Rsin(θt/2)sin((θt/2)+θ-θr-(2π(i-1)/5))+r,θr+π(2i-2)5≤θ<θr+π(2i-1)5,-Rsin(θt/2)sin((θt/2)-θ+θr-(2π(i-1)/5))+r,θr+π(2i-3)5≤θ<θr+π(2i-2)5,
with θt=π/5, θr=π/7, R=6/7, and i=1,2,3,4,5(23)β-(x,y)=(x2+5xy+1x+2y2+3),u-(x,y)=6+sin(6πx)sin(6πy).The numerical result with our method using a40×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.
n
x
×
n
y
Error inu-
Order
20 × 20
0.51
40 × 40
0.20
1.32
80 × 80
0.061
1.73
160 × 160
0.016
1.90
320 × 320
0.0042
1.96Figure 5
Numerical result of the problem with “star” shape.Example 2.
This example has a “happy face” interface. The level-set functionϕ(x,y), the coefficients β-, and the solution u- are given as follows:
(2)ϕ(x,y)=max(min(ϕ1,ϕ2,ϕ3),ϕ4,ϕ5,ϕ6,min(ϕ7,ϕ8)),ϕ1(x,y)=x2+y2-0.752-0.152,ϕ2(x,y)=(x-0.75)2+y2-0.152,ϕ3(x,y)=(x+0.75)2+y2-0.152,ϕ4(x,y)=-0.10.12(x-0.2)2-0.120.1(y-0.22)2+0.12·0.1,ϕ5(x,y)=-0.10.12(x+0.2)2-0.120.1(y-0.22)2+0.12·0.1,ϕ6(x,y)=-x2-(y+0.08)2+0.122,ϕ7(x,y)=-x2-(y+0.625)2+0.4252,ϕ8(x,y)=-x2-(y+0.25)2+0.22,β-(x,y)=(x2-y2+3x2-y2+1x2-y2+1x2-y2+4),u-(x,y)=5-5x2-5y2.The numerical result with our method using a40×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.
n
x
×
n
y
Error inu-
Order
20 × 20
0.016
40 × 40
0.0041
1.97
80 × 80
0.0010
1.97
160 × 160
2.58
e
-
004
2.02Figure 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 u- are given as follows:
(25)ϕ(x,y)=min(min(min3,min4))ϕ1,min(min(min1,min2),min(min3,min4))),ϕ1(x,y)=x2+y2-0.02,ϕ2(x,y)=x2+y2-0.08,ϕ3(x,y)=x2+y2-0.15,ϕ4(x,y)=x2+y2-0.25,ϕ5(x,y)=x2+y2-0.35,ϕ6(x,y)=x2+y2-0.45,ϕ7(x,y)=x2+y2-0.55,ϕ8(x,y)=x2+y2-0.80,ϕ9(x,y)=x2+y2-1,min1=-min(ϕ2,-ϕ3),min2=-min(ϕ4,-ϕ5),min3=-min(ϕ6,-ϕ7),min4=-min(ϕ8,-ϕ9),β-(x,y)=(2+sin(x+y)002+sin(x+y)),u-(x,y)=x2+y2+sin(x+y).The numerical result with our method using a40×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.
n
x
×
n
y
Error inu-
Order
20 × 20
0.0034
40 × 40
8.44
e
-
00
4
1.99
80 × 80
2.14
e
-
00
4
1.98
160 × 160
5.35
e
-
00
5
2.00Figure 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 u- are given as follows:
(26)ϕ(x,y,z)=max(min3,ϕ6),ϕ1(x,y,z)=x2+y2-0.02,ϕ2(x,y,z)=x2+y2-0.1,ϕ3(x,y,z)=x2+y2-0.3,ϕ4(x,y,z)=x2+y2-0.6,ϕ5(x,y,z)=x2+y2-0.89,ϕ6(x,y,z)=x2+y2+z2-0.9,min1=-min(ϕ2,-ϕ3),min2=-min(ϕ4,-ϕ5),min3=min(ϕ1,min(min1,min2)),β-(x,y,z)=(4x2+6sin(x+y)xsin(x+y)2z2+30.5sin(x)x0.5sin(x)cos(z)+5),u-(x,y,z)=5-sin(x2)+3y2+z.Figure8 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 Example4.
n
x
×
n
y
×
n
z
Error inu-
Order
10 × 10 × 10
0.016
20 × 20 × 20
0.0050
1.70
40 × 40 × 40
0.0013
1.96
80 × 80 × 80
0.00032
2.02Figure 8
Numerical result of Example4.Example 5.
The level-set functionϕ, the coefficients β-, and the solution u- are given as follows:
(27)ϕ(x,y,z)=min(min(ϕ1,ϕ2),min(ϕ3,ϕ4)),ϕ1(x,y,z)=-0.25+x2+y2+(z+0.5)2,ϕ2(x,y,z)=-0.25+(x-0.4)2+y2+z2,ϕ3(x,y,z)=-0.25+(x+0.3)2+y2+z2,ϕ4(x,y,z)=-0.25+x2+(y+0.5)2+z2,β-(x,y,z)=(cos(x+y)2+3z0.2sin(z-x)zz2+5y0.2sin(z-x)ysin(z)2+2),u-(x,y,z)=exp(-(x2+y2+z2)20).Figure9 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 Example5.
n
x
×
n
y
×
n
z
Error inu-
Order
10 × 10 × 10
0.00090
20 × 20 × 20
0.00025
1.86
40 × 40 × 40
0.000062
2.01
80 × 80 × 80
0.000015
2.00Figure 9
Numerical result of Example5.Example 6.
The level-set functionϕ, the coefficients β-, and the solution u- are given as follows:
(28)ϕ(x,y,z)=min((x-0.39)2+y2+z2-0.25,(x+0.39)2+y2+z2-0.25),β-(x,y,z)=(4x2+6sin(x+y)xysin(x+y)2z2+30.5sin(x)xy0.5sin(x)cos(xy+z)2+5),u-(x,y,z)=10-2x3+3y2+sin(z-y).Figure10 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 Example6.
n
x
×
n
y
×
n
z
Error inu-
Order
10 × 10 × 10
0.041
20 × 20 × 20
0.011
1.92
40 × 40 × 40
0.0029
1.88
80 × 80 × 80
0.00076
1.95Figure 10
Numerical result of Example6.


## 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 u- are given as follows:
(22)ϕ(r,θ)={-Rsin(θt/2)sin((θt/2)+θ-θr-(2π(i-1)/5))+r,θr+π(2i-2)5≤θ<θr+π(2i-1)5,-Rsin(θt/2)sin((θt/2)-θ+θr-(2π(i-1)/5))+r,θr+π(2i-3)5≤θ<θr+π(2i-2)5,
with θt=π/5, θr=π/7, R=6/7, and i=1,2,3,4,5(23)β-(x,y)=(x2+5xy+1x+2y2+3),u-(x,y)=6+sin(6πx)sin(6πy).The numerical result with our method using a40×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.
n
x
×
n
y
Error inu-
Order
20 × 20
0.51
40 × 40
0.20
1.32
80 × 80
0.061
1.73
160 × 160
0.016
1.90
320 × 320
0.0042
1.96Figure 5
Numerical result of the problem with “star” shape.Example 2.
This example has a “happy face” interface. The level-set functionϕ(x,y), the coefficients β-, and the solution u- are given as follows:
(2)ϕ(x,y)=max(min(ϕ1,ϕ2,ϕ3),ϕ4,ϕ5,ϕ6,min(ϕ7,ϕ8)),ϕ1(x,y)=x2+y2-0.752-0.152,ϕ2(x,y)=(x-0.75)2+y2-0.152,ϕ3(x,y)=(x+0.75)2+y2-0.152,ϕ4(x,y)=-0.10.12(x-0.2)2-0.120.1(y-0.22)2+0.12·0.1,ϕ5(x,y)=-0.10.12(x+0.2)2-0.120.1(y-0.22)2+0.12·0.1,ϕ6(x,y)=-x2-(y+0.08)2+0.122,ϕ7(x,y)=-x2-(y+0.625)2+0.4252,ϕ8(x,y)=-x2-(y+0.25)2+0.22,β-(x,y)=(x2-y2+3x2-y2+1x2-y2+1x2-y2+4),u-(x,y)=5-5x2-5y2.The numerical result with our method using a40×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.
n
x
×
n
y
Error inu-
Order
20 × 20
0.016
40 × 40
0.0041
1.97
80 × 80
0.0010
1.97
160 × 160
2.58
e
-
004
2.02Figure 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 u- are given as follows:
(25)ϕ(x,y)=min(min(min3,min4))ϕ1,min(min(min1,min2),min(min3,min4))),ϕ1(x,y)=x2+y2-0.02,ϕ2(x,y)=x2+y2-0.08,ϕ3(x,y)=x2+y2-0.15,ϕ4(x,y)=x2+y2-0.25,ϕ5(x,y)=x2+y2-0.35,ϕ6(x,y)=x2+y2-0.45,ϕ7(x,y)=x2+y2-0.55,ϕ8(x,y)=x2+y2-0.80,ϕ9(x,y)=x2+y2-1,min1=-min(ϕ2,-ϕ3),min2=-min(ϕ4,-ϕ5),min3=-min(ϕ6,-ϕ7),min4=-min(ϕ8,-ϕ9),β-(x,y)=(2+sin(x+y)002+sin(x+y)),u-(x,y)=x2+y2+sin(x+y).The numerical result with our method using a40×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.
n
x
×
n
y
Error inu-
Order
20 × 20
0.0034
40 × 40
8.44
e
-
00
4
1.99
80 × 80
2.14
e
-
00
4
1.98
160 × 160
5.35
e
-
00
5
2.00Figure 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 u- are given as follows:
(26)ϕ(x,y,z)=max(min3,ϕ6),ϕ1(x,y,z)=x2+y2-0.02,ϕ2(x,y,z)=x2+y2-0.1,ϕ3(x,y,z)=x2+y2-0.3,ϕ4(x,y,z)=x2+y2-0.6,ϕ5(x,y,z)=x2+y2-0.89,ϕ6(x,y,z)=x2+y2+z2-0.9,min1=-min(ϕ2,-ϕ3),min2=-min(ϕ4,-ϕ5),min3=min(ϕ1,min(min1,min2)),β-(x,y,z)=(4x2+6sin(x+y)xsin(x+y)2z2+30.5sin(x)x0.5sin(x)cos(z)+5),u-(x,y,z)=5-sin(x2)+3y2+z.Figure8 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 Example4.
n
x
×
n
y
×
n
z
Error inu-
Order
10 × 10 × 10
0.016
20 × 20 × 20
0.0050
1.70
40 × 40 × 40
0.0013
1.96
80 × 80 × 80
0.00032
2.02Figure 8
Numerical result of Example4.Example 5.
The level-set functionϕ, the coefficients β-, and the solution u- are given as follows:
(27)ϕ(x,y,z)=min(min(ϕ1,ϕ2),min(ϕ3,ϕ4)),ϕ1(x,y,z)=-0.25+x2+y2+(z+0.5)2,ϕ2(x,y,z)=-0.25+(x-0.4)2+y2+z2,ϕ3(x,y,z)=-0.25+(x+0.3)2+y2+z2,ϕ4(x,y,z)=-0.25+x2+(y+0.5)2+z2,β-(x,y,z)=(cos(x+y)2+3z0.2sin(z-x)zz2+5y0.2sin(z-x)ysin(z)2+2),u-(x,y,z)=exp(-(x2+y2+z2)20).Figure9 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 Example5.
n
x
×
n
y
×
n
z
Error inu-
Order
10 × 10 × 10
0.00090
20 × 20 × 20
0.00025
1.86
40 × 40 × 40
0.000062
2.01
80 × 80 × 80
0.000015
2.00Figure 9
Numerical result of Example5.Example 6.
The level-set functionϕ, the coefficients β-, and the solution u- are given as follows:
(28)ϕ(x,y,z)=min((x-0.39)2+y2+z2-0.25,(x+0.39)2+y2+z2-0.25),β-(x,y,z)=(4x2+6sin(x+y)xysin(x+y)2z2+30.5sin(x)xy0.5sin(x)cos(xy+z)2+5),u-(x,y,z)=10-2x3+3y2+sin(z-y).Figure10 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 Example6.
n
x
×
n
y
×
n
z
Error inu-
Order
10 × 10 × 10
0.041
20 × 20 × 20
0.011
1.92
40 × 40 × 40
0.0029
1.88
80 × 80 × 80
0.00076
1.95Figure 10
Numerical result of Example6.


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



---

*Source: 967342-2013-12-24.xml*