A Numerical Method for Solving 3D Elasticity Equations with Sharp-Edged Interfaces

Abstract

Interface problems occur frequently when two or more materials meet. Solving elasticity equations with sharp-edged interfaces in three dimensions is a very complicated and challenging problem for most existing methods. There are several difficulties: the coupled elliptic system, the matrix coefficients, the sharp-edged interface, and three dimensions. An accurate and efficient method is desired. In this paper, an efficient nontraditional finite element method with nonbody-fitting grids is proposed to solve elasticity equations with sharp-edged interfaces in three dimensions. The main idea is to choose the test function basis to be the standard finite element basis independent of the interface and to choose the solution basis to be piecewise linear satisfying the jump conditions across the interface. The resulting linear system of equations is shown to be positive definite under certain assumptions. Numerical experiments show that this method is second order accurate in the 𝐿^∞ norm for piecewise smooth solutions. More than 1.5th order accuracy is observed for solution with singularity (second derivative blows up).

1. Introduction

The macroscopic mechanical response of most elastic materials, such as wood, metal, rubber, rock, and bone, can be described by the elasticity theory. However, the numerical solution for an elasticity problem with interface is highly nontrivial, because the coefficients in the partial differential equation for elasticity system are discontinuous, and the solution to the PDE system also needs to satisfy the jump conditions across each material interface. Examples of elasticity interface problem include the minimum-compliance-design problems [1], the microstructural evolution [2], and the atomic interactions [3]. Designing an accurate and efficient method for these problems is a difficult job, especially when the interface is not smooth.

In this paper, we consider a three-dimensional elasticity problem with sharp-edged interfaces and matrix coefficients. The model of interest is as follows.

Consider an open bounded domain Ω⊂𝑅^𝑑. Let Γ be an interface of codimension 𝑑−1, which divides Ω into disjoint open subdomains, Ω⁻ and Ω⁺, hence Ω=Ω⁻⋃Ω⁺⋃Γ. Assume that the boundary 𝜕Ω and the boundary of each subdomain 𝜕Ω^± are Lipschitz continuous as submanifolds. Since 𝜕Ω^± are Lipschitz continuous, so is Γ. A unit normal vector of Γ can be defined as a.e. on Γ; see Section 1.5 in [4].

We seek solutions of the variable coefficient elliptic equation away from the interface Γ given by −∇⋅𝛽₁(𝑥)∇𝑢₁(𝑥)−∇⋅𝛽₂(𝑥)∇𝑢₂(𝑥)=𝑓₁(𝑥),−∇⋅𝛽₃(𝑥)∇𝑢₁(𝑥)−∇⋅𝛽₄(𝑥)∇𝑢₂(𝑥)=𝑓₂(𝑥),𝑥∈Ω⧵Γ, (1) in which 𝑥=(𝑥₁,…,𝑥_𝑑) denotes the spatial variables and ∇ is the gradient operator. 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 𝑓_𝑖(𝑥), 𝑖=1,2 is assumed to lie in 𝐿²(Ω).

Given functions 𝑎 and 𝑏 along the interface Γ, we prescribe the jump conditions: 𝑢_1Γ(𝑥)≡𝑢⁺₁(𝑥)−𝑢⁻₁(𝑥)=𝑎₁(𝑥),𝑢_2Γ(𝑥)≡𝑢⁺₂(𝑥)−𝑢⁻₂(𝑥)=𝑎₂(𝑥),𝑛⋅𝛽⁺₁(𝑥)∇𝑢⁺₁(𝑥)+𝛽⁺₂(𝑥)∇𝑢⁺₂(𝑥)−𝑛⋅𝛽⁻₁(𝑥)∇𝑢⁻₁(𝑥)+𝛽⁻₂(𝑥)∇𝑢⁻₂(𝑥)=𝑏₁(𝑥),𝑛⋅𝛽⁺₃(𝑥)∇𝑢⁺₁(𝑥)+𝛽⁺₄(𝑥)∇𝑢⁺₂(𝑥)−𝑛⋅𝛽⁻₃(𝑥)∇𝑢⁻₁(𝑥)+𝛽⁻₄(𝑥)∇𝑢⁻₂(𝑥)=𝑏₂(𝑥). (2)

The “±” superscripts refer to limits taken from within the subdomains Ω^±.

Finally, we prescribe boundary conditions 𝑢₁(𝑥)=𝑔₁(𝑥),𝑢₂(𝑥)=𝑔₂(𝑥),𝑥∈𝜕Ω, (3) for a given function 𝑔 on the boundary 𝜕Ω.

The setup of the problem is illustrated in Figure 1.

Figure 1 

Setup of the problem.

We consider the three-dimensional elasticity interface problem (𝑑=3) in this paper. Although the real problem would have a three-by-three coupled PDE system and six jump conditions, for simplicity of coding, we use two-by-two PDE system to demonstrate the effectiveness of our method. Note that the mathematical challenge is in the 3D geometry when generalized from the corresponding 2D problem. Changing the system from two-by-two to three-by-three does not bring in new mathematical challenge.

Interface problems have been studied by many researchers and several techniques have been developed. Finite element methods with body-fitted grid use a mesh triangulation to capture the boundary [5–8]. However, in many situations, such as when the boundary is moving, the mesh generation may be both computationally expensive and challenging. A uniform Cartesian grid is preferred.

In [9, 10], to simulate the flow patten of blood in the heart, Peskin proposed the “immersed boundary” method, which used numerical approximation of 𝛿-function for singular sources on the interface. In [11], in order to compute two-phase flow, a level-set method was combined with the “immersed boundary” method. These methods are simple to be used but difficult to achieve high-order accuracy.

In addition, 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. Using finite difference scheme typically requires taking high-order derivatives of jump conditions and interface in Taylor expansion. Also property of the discretized linear system is hard to analyze for interface problem with general jump condition. The “immersed interface” method was proposed in [12]. This method 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. Second-order differentiation of the interface is needed. The Matched Interface and Boundary (MIB) method was developed in [13] and improved to handle sharp-edged interfaces in two dimensions in [14] and in three dimensions in [15]. With an elegant treatment, second-order accuracy was achieved in the 𝐿^∞ norm.

The existing finite element schemes on unfitted meshes are usually designed by modifying the finite element basis near the interface. Examples are immersed finite element method [16, 17], adaptive immersed interface method [18], and extended finite element method [19–21].

An elasticity system can be solved by both finite difference and finite element methods. In [22, 23], a finite difference method called immersed interface method is proposed to solve elasticity problems with nonhomogeneous jump conditions. While second-order accuracy was achieved, the condition number of the discrete system is quite large especially in the nearly incompressible case (𝜆 is large) compared with that obtained from finite element formulations.

In [17], an immersed-interface finite element method was developed for scalar elliptic interface problems with nonhomogeneous jump conditions. In [24], the immersed-interface finite element method was developed for the elasticity system with homogeneous and nonhomogeneous jump conditions. The Soblev space theory provides strong theoretical foundations for convergence analysis for finite element methods.

There are some other approaches to solve the elliptic interface problems. In particular, some recent work [14, 25, 26] can handle sharp-edged interfaces. However, these have not been developed to solve the elasticity interface problems.

In [27], a nontraditional finite element method is proposed to solve the 2D elasticity interface problems with sharp-edged interfaces.

In this paper, based on the method in [27–30], we propose a numerical method for solving the 3D elasticity problem with sharp-edged interfaces. Compared with [27], this is not a trivial extension, as the geometry in three dimensions is more complicated. In two dimensions, the interface could only cut a triangle in one way (if not counting the case the interface hits a grid point as a separate case). In three dimension, There are two types of plane segments see Figure 3. We proved that the resulting linear system is nonsymmetric but positive definite under certain assumption. The method is simpler compared with that developed in [24] and can be applied for more general problems since we allow 𝛽_𝑖 to be matrices. The condition number grows with 𝑂(ℎ⁻²), the same as the case without interface. The condition number also grows only linearly with the ratio between 𝛽⁺ and 𝛽⁻.

The rest of the paper is organized as follows: Section 2 presents the weak form of the elasticity system. In Section 3, we discuss our new numerical method and some theoretical analysis. In Section 4, some examples are presented to demonstrate the performance of our method. We conclude in Section 5.

2. The Weak Formulation

We modify the weak formulation in [28, 30, 31]. We are going to use the usual Sobolev spaces 𝐻¹(Ω). For 𝐻¹₀(Ω), instead of the usual inner product, we choose one which is better suited to our problem: 𝐵[]=∫𝑢,𝑣^Ω+𝛽₁∇𝑢₁⋅∇𝑣₁+𝛽₂∇𝑢₂⋅∇𝑣₁+∫_Ω−𝛽₁∇𝑢₁⋅∇𝑣₁+𝛽₂∇𝑢₂⋅∇𝑣₁,∫_Ω+𝛽₃∇𝑢₁⋅∇𝑣₂+𝛽₄∇𝑢₂⋅∇𝑣₂+∫_Ω−𝛽₃∇𝑢₁⋅∇𝑣₂+𝛽₄∇𝑢₂⋅∇𝑣₂. (4)

Let Γ be any closed Lipschitz continuous hypersurface of dimension 𝑑−1 in Ω, where the overline denotes the closure of a set. Let 𝑅_𝜕Ω⁻ denote the restriction operator from 𝐻¹(Ω) to 𝐿²(𝜕Ω⁻). This restriction operator 𝑅_𝜕Ω⁻ is well defined and bounded, because it is closed Lipschitz continuous (see Theorem  2.4.2 in [32]) and 𝐶¹(Ω) is dense in 𝐻¹(Ω). Throughout this section, we shall always assume that our interface data 𝑎 and 𝑏 are the restrictions of functions 𝑎 and 𝑏̃∈𝐻¹(Ω) on 𝜕Ω⁻ and then limited on Γ, respectively. That is on Γ, 𝑎=𝑅_𝜕Ω⁻(𝑎),𝑏=𝑅_𝜕Ω⁻𝑏̃. (5) We shall always assume that our boundary data 𝑔 can be obtained: assume that there exists a function 𝑐̃∈𝐻¹(Ω) so that 𝑔 is given as, on 𝜕Ω, 𝑅𝑔=^𝜕Ω(𝑐̃−𝑎),on𝜕Ω⋂𝜕Ω⁻,𝑅_𝜕Ω(𝑐)̃,on𝜕Ω⧵𝜕Ω⁻. (6)

Equation (6) could be thought as a compatibility condition between 𝑎 and 𝑔. To simplify the notation, from now on, we will drop the tildes.

We will construct a unique solution of the problem in the class 𝐻(𝑎,𝑐)=𝑢∶𝑢−𝑐+𝑎𝜒Ω⁻∈𝐻¹₀(Ω). (7)

If 𝑢∈𝐻(𝑎,𝑐), then [𝑢]_Γ=𝑎 and 𝑢|_𝜕Ω=𝑔. Note that 𝐻¹₀(Ω) can be identified with 𝐻(0,0).

Definition 1. A function 𝑢∈𝐻(𝑎,𝑐) is a weak solution of (1)–(3), if 𝑣=𝑢−𝑐+𝑎𝜒(Ω⁻)∈𝐻¹₀(Ω) satisfies 𝐵𝑣,𝜓=𝐹𝜓, (8) for all 𝜓∈𝐻¹₀(Ω), where 𝐵𝑣,𝜓=_Ω+𝛽∇𝑣⋅∇𝜓+_Ω−𝐹𝜓=𝛽∇𝑣⋅∇𝜓,_Ω𝑓𝜓+_Ω𝛽∇𝑐⋅∇𝜓+_Ω−𝛽∇𝑎⋅∇𝜓+_Γ𝑏𝜓, (9) or equivalently.

Definition 2. A function 𝑢∈𝐻(𝑎,𝑐) is a weak solution of (1)–(3), if 𝑢 satisfies, for all 𝜓∈𝐻¹₀(Ω), _Ω+𝛽₁∇𝑢₁⋅∇𝜓₁+𝛽₂∇𝑢₂⋅∇𝜓₁+_Ω−𝛽₁∇𝑢₁⋅∇𝜓₁+𝛽₂∇𝑢₂⋅∇𝜓₁=_Ω𝑓₁𝜓₁+_Γ𝑏₁𝜓₁,_Ω+𝛽₃∇𝑢₁⋅∇𝜓₂+𝛽₄∇𝑢₂⋅∇𝜓₂+_Ω−𝛽₃∇𝑢₁⋅∇𝜓₂+𝛽₄∇𝑢₂⋅∇𝜓₂=_Ω𝑓₂𝜓₂+_Γ𝑏₂𝜓₂. (10)
A classical solution of (1)–(3), 𝑢|_Ω^±∈𝐶²(Ω^±), is necessarily a weak solution. Because all the subdomains’ boundaries 𝜕Ω^± are Lipschitz continuous, the integration by parts is legal in each subdomain, Ω^±; see Theorem  1.5.3.1 in [33].

We have the following theorem.

Theorem 3. If 𝑓∈𝐿²(Ω), and 𝑎, 𝑏, and 𝑐∈𝐻¹(Ω), then there exists a unique weak solution of (1)–(3) in 𝐻(𝑎,𝑐).

Proof. See Theorem 2.1 in [31].

3. Numerical Method

We define 𝑢𝑢≡¹𝑢₂𝑓,𝑓≡¹𝑓₂𝑔,𝑔≡¹𝑔₂,𝑎𝑎≡¹𝑎₂𝑏,𝑏≡¹𝑏₂𝛽,𝛽≡¹𝛽₂𝛽₃𝛽₄, (11) and choose test function 𝜓𝜓=¹00𝜓or₂ (12) and redefine the gradient and divergence operator 𝜕∇≡0𝜕𝜕𝑥0𝜕𝜕𝑦00𝜕𝜕𝑧0𝜕𝜕𝑥0𝜕𝜕𝑦,𝜕𝜕𝑧∇⋅≡𝜕𝜕𝑥𝜕𝜕𝑦𝜕𝜕𝑧000000𝜕𝜕𝑥𝜕𝜕𝑦.𝜕𝑧 (13)

Then (1) can be written as −∇⋅𝛽(𝑥)∇𝑢(𝑥)=𝑓(𝑥),𝑥∈Ω⧵Γ. (14) The jump condition equation (2) can be reformed as [𝑢]_Γ(𝑥)≡𝑢⁺(𝑥)−𝑢⁻(𝑥)(𝑥),𝛽=𝑎𝑛⋅⁺(𝑥)∇𝑢⁺(𝑥)−𝑛⋅𝛽⁻(𝑥)∇𝑢⁻(𝑥)=𝑏(𝑥), (15) and boundary condition is 𝑢(𝑥)=𝑔(𝑥)𝑥∈𝜕Ω. (16)

For simplicity, we will discuss the following properties under the form of (14), (15), and (16).

For ease of discussion, Section 3, and for accuracy testing in Section 4, we assume that 𝑎, 𝑏, and 𝑐 are smooth on the closure of Ω, 𝛽 and 𝑓 are smooth on Ω⁺ and Ω⁻, but they may be discontinuous across the interface Γ. However, 𝜕Ω, 𝜕Ω⁻, and 𝜕Ω⁺ are kept to be Lipschitz continuous. 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 Ω⁺.

In this paper, we restrict ourselves to a cube cell domain Ω=(𝑥_min,𝑥_max)×(𝑦_min,𝑦_max)×(𝑧_min,𝑧_max) in 3D space, and 𝛽 is a 3×3 matrix that is uniformly elliptic in each subdomain. Given positive integers 𝐼, 𝐽, and 𝐾, set Δ𝑥=(𝑥_max−𝑥_min)/𝐼, Δ𝑦=(𝑦_max−𝑦_min)/𝐽, and Δ𝑧=(𝑧_max−𝑧_min)/𝐾. We define a uniform Cartesian grid (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)=(𝑥_min+𝑖Δ𝑥,𝑦_min+𝑗Δ𝑦,𝑧_min+𝑘Δ𝑧) for 𝑖=0,…,𝐼, 𝑗=0,…,𝐽 and 𝑘=0,…,𝐾. Each (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) is called a grid point. For the case 𝑖=0,𝐼, 𝑗=0,𝐽, or 𝑘=0,𝐾, a grid point is called a boundary point, 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≤𝑗≤𝐽,0≤𝑘≤𝐾_01,ℎ=𝜔^ℎ=𝜔_{𝑖,𝑗,𝑘}∈𝐻_±^1,ℎ∶𝜔_{𝑖,𝑗,𝑘}.=0if𝑖=0,𝐼or𝑗=0,𝐽or𝑘=0,𝐾 (17)

We cut every cube cell region [𝑥_𝑖,𝑥_𝑖+1]×[𝑦_𝑗,𝑦_𝑗+1]×[𝑧_𝑘,𝑧_𝑘+1] into six tetrahedron regions. Collecting all those tetrahedron regions, we obtain a uniform tetrahedralization 𝑇^ℎ∶⋃_{𝐿∈𝑇ℎ}𝐿; see Figure 2.

Figure 2 

Cubic mesh and the tetrahedralization.

Figure 3 

Interface tetrahedron, left: Case 1, right: Case 2.

If 𝜙(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘)≤0, we count the grid point (𝑥_𝑖,𝑦_𝑗,𝑧_𝑘) as in Ω⁻; otherwise we count it as in Ω⁺. We call an edge (an edge of a tetrahedron in the tetrahedralization) an interface edge if two of its ends (vertices of tetrahedrons in the tetrahedralization) belong to different subdomains; otherwise we call it a regular edge.

We call a cell 𝐿 an interface cell if its vertices belong to different subdomains. In the interface cell, we write 𝐿=𝐿⁺⋃𝐿⁻. 𝐿⁺ and 𝐿⁻ are separated by a plane segment, denoted by Γ^ℎ_𝐿. There are two kinds of plane segment; see Figure 3. The vertices of the plane segment Γ^ℎ_𝐿 are located on the interface Γ, and their locations can be calculated from the linear interpolations of the discrete level-set functions 𝜙^ℎ=𝜙(𝑥_𝑖,𝑦_𝑗,𝑧_𝑘). The vertices of 𝐿⁺ are located in Ω⁺⋃Γ and the vertices of 𝐿⁻ are located in Ω⁻⋃Γ. 𝐿⁺ and 𝐿⁻are approximations of the regions of 𝐿⋂Ω⁺ and 𝐿⋂Ω⁻, respectively. We call a cell 𝐿 a regular cell if all its vertices belong to the same subdomain, either Ω⁺ or Ω⁻. For a regular cell, we also write 𝐿=𝐿⁺⋃𝐿⁻, where 𝐿⁻={} (empty set) if all vertices of 𝐿 are in Ω⁺ and 𝐿⁺={} if all vertices of 𝐿 are in Ω⁻. Clearly Γ^ℎ_𝐿={} in a regular cell, and 𝐿⁺ and 𝐿⁻ are approximations of the regions 𝐿⋂Ω⁺ and 𝐿⋂Ω⁻, respectively. We use |𝐿⁺| and |𝐿⁻| to represent the volumes of 𝐿⁺ and 𝐿⁻, respectively.

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 tetrahedron 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 𝐿⁻. The location of its discontinuity in the interface cell is the plane segment Γ^ℎ_𝐿. Note that the vertices of the plane segment are located on the interface Γ, and hence the interface condition [𝑢]=𝑎 is enforced exactly at the vertices of Γ^ℎ_𝐿. In each interface cell, the interface condition [𝛽∇𝑢⋅𝑛]=𝑏 is enforced with the value 𝑏 at the barycenter of Γ^ℎ_𝐿 for Case 1 and at the points labeled “9” and “10” in Figure 3 for Case 2. Clearly such a function is not continuous in general, and neither is the set of such functions a linear space. We denote the set of such functions as 𝐻^1,ℎ_𝑎,𝑐, which should be thought as an approximation of the solution class 𝐻(𝑎,𝑐) (see (7)) plus the restriction of [𝛽∇𝑢⋅𝑛]=𝑏. Similar versions of such extension can be found in the literature [16, 34]. In Case 1, the number of unknowns matches with the number of the equations and the extension 𝑈^ℎ is unique. In Case 2, the local construction problem is solved in a least square sense. The values of the points on the interface can be represented as a linear combination of a constant and the 8 unknown value of the vertices on the tetrahedron cell.

We propose the following method.

Method 4. Find a discrete function 𝑢^ℎ∈𝐻^1,ℎ such that 𝑢^ℎ=𝑔^ℎ on boundary points and so that for all 𝜓^ℎ∈𝐻₀^1,ℎ, we have _{𝐿∈𝑇ℎ𝐿+}𝛽∇𝑈^ℎ𝑢^ℎ⋅∇𝑇^ℎ𝜓^ℎ+_𝐿−𝛽∇𝑈^ℎ𝑢^ℎ⋅∇𝑇^ℎ𝜓^ℎ=_{𝐿∈𝑇ℎ𝐿+}𝑓𝑇^ℎ𝜓^ℎ+_𝐿−𝑓𝑇^ℎ𝜓^ℎ+_Γℎ𝐿𝑏𝑇^ℎ𝜓^ℎ. (18) Note that 𝑢=𝑔 on the boundary is the same as 𝑢−𝑐+𝑎𝜒(Ω⁻)=0 on the boundary.

Since our solution bases and test function bases are different, the matrix 𝐴 for the linear system generated by Method 4 is not symmetric in the presence of an interface. However, we can prove that it is positive definite under certain assumptions. In all our numerical examples, 𝐴 is positive definite numerically.

Theorem 5. If 𝛽 is continuous and positive definite, then the matrix 𝐴 for the linear system generated by Method 4 is positive definite.

Proof. For vector 𝑐∈𝑅^2𝑛 we have 𝑐^𝑇𝐴𝑐=𝐵^2𝑛_𝑖=1𝑐_𝑖𝑢_𝑖,^2𝑛_𝑖=1𝑐_𝑖𝜓_𝑖. (19) Here we can choose specific jump conditions 𝑎₁=𝑎₂=𝑏₁=𝑏₂=0 since changing 𝑎_𝑖, 𝑏_𝑖 will not affect the matrix 𝐴. In this case we have ^2𝑛_𝑖=1𝑐_𝑖𝑢_𝑖=^2𝑛_𝑖=1𝑐_𝑖𝜓_𝑖=𝑤, (20) since 𝐵[𝑤,𝑤]=_Ω𝛽(𝑥)∇𝑤(𝑥)^𝑇∇𝑤(𝑥)𝑑𝑥. (21) So if 𝛽 is positive definite, then ∇𝑤(𝑥)^𝑇𝛽^𝑇∇𝑤(𝑥)>0, for all 𝑥∈Ω. Therefore 𝐵[𝑤,𝑤]>0, which implies 𝑐^𝑇𝐴𝑐>0. Thus, 𝐴 is positive definite.

In some applications [24], the matrix 𝛽 is only semipositive definite with zero determinant. The previous Theorem does not apply. We prove later that when the matrix 𝛽 is of certain form that frequently appeared in applications and semipositive definite, the matrix 𝐴 generated by our method is still positive definite.

Theorem 6. If 𝜆>0, 𝜇>0 and 𝛽₁=^{𝜆+2𝜇000𝜇0}_00𝜇, 𝛽₂=^00𝜆0𝜇0_𝜇00, 𝛽₃=^00𝜇0𝜇0_𝜆00, and 𝛽₄=^𝜇000𝜇0_00𝜆+2𝜇, then the matrix 𝐴 for the linear system generated by Method 4 is positive definite.

Proof. Suppose for a contradiction that 𝐴 is not positive definite. Then there is a vector 𝑐∈𝑅^2𝑛 and 𝑐≠0 such that 𝑐^𝑇𝐴𝑐≤0. Let 𝑤𝑤=¹𝑤₂=^2𝑛_𝑖=1𝑐_𝑖𝜓_𝑖=^2𝑛_𝑖=1𝑐_𝑖𝑢_𝑖. (22) Let 𝐷_𝑤=𝜕𝑤₁𝜕𝑥𝜕𝑤₁𝜕𝑦𝜕𝑤₁𝜕𝑧𝜕𝑤₂𝜕𝑥𝜕𝑤₂𝜕𝑦𝜕𝑤₂,𝛽𝜕𝑧₁₂₃₄=,𝐺𝜆+2𝜇0000𝜆0𝜇00𝜇000𝜇𝜇0000𝜇𝜇000𝜇00𝜇0𝜆0000𝜆+2𝜇_𝑤=𝜕𝑤₁𝜕𝑥𝜕𝑤₁𝜕𝑦𝜕𝑤₁𝜕𝑧𝜕𝑤₂𝜕𝑥𝜕𝑤₂𝜕𝑦𝜕𝑤₂.𝜕𝑧 (23) Then 𝐵[]⟹𝑤,𝑤≤0_Ω𝛽∇𝑤(𝑥)^𝑇(𝑥)𝑑⟹∇𝑤𝑥≤0_Ω𝐷_𝑤𝛽₁₂₃₄𝐺_𝑤𝑑𝑥≤0. (24)
Since for all 𝑎=[𝑎₁,𝑎₂,𝑎₃,𝑎₄,𝑎₅,𝑎₆]^𝑇∈𝑅⁶, 𝑎^𝑇𝛽𝑎=𝑎₁+𝑎₆²𝜆+2𝑎²₁+𝑎²₆𝜇+𝑎₂+𝑎₅²𝜇+𝑎₃+𝑎₄²𝜇≥0. (25)
So 𝑎^𝑇𝛽𝑎=0 if and only if 𝑎₁=𝑎₆=0, 𝑎₂=−𝑎₅, and 𝑎₃=−𝑎₄. Then (𝜕𝑤₁/𝜕𝑥)(𝑥)=𝑎₁=0, for all 𝑥∈Ω.
However, 𝑤₁=∑^𝑛_𝑖=1𝑐_𝑖𝜓¹_𝑖 implies that 𝜕𝑤₁/𝜕𝑥=∑^𝑛_𝑖=1𝑐_𝑖(𝜕𝜓¹_𝑖/𝜕𝑥).
Since 𝑐=[𝑐₁,𝑐₂,…,𝑐_2𝑛]^𝑇≠0, without loss of generality, we assume that 𝑐₁≠0. If we choose point 𝑥∈Ω such that 𝜕𝜓¹₁(𝑥)/𝜕𝑥≠0 and 𝜕𝜓¹_𝑖(𝑥)/𝜕𝑥=0, 𝑖=2,3,…,𝑛, then ∑^𝑛_𝑖=1𝑐_𝑖(𝜕𝜓¹_𝑖/𝜕𝑥)≠0, a contradiction.
Therefore 𝑐^𝑇𝐴𝑐>0 for all 𝑐≠0; thus, 𝐴 is positive definite.

Remark 7. A positive definite matrix has positive determinant and is therefore invertible. It also has an 𝐿𝐷𝑀^𝑇 factorization where 𝐷=diag(𝑑_𝑖) and 𝑑_𝑖>0, and 𝐿,𝑀 are lower triangular. The linear system 𝐴𝑥=𝑏 can be solved efficiently.

4. Numerical Experiments

Consider the problem −∇⋅𝛽₁∇𝑢₁−∇⋅𝛽₂∇𝑢₂=𝑓₁,inΩ^±,−∇⋅𝛽₃∇𝑢₁−∇⋅𝛽₄∇𝑢₂=𝑓₂,inΩ^±,𝑢₁=𝑎₁𝑢,onΓ,₂=𝑎₂𝛽,onΓ,₁∇𝑢₁+𝛽₂∇𝑢₂⋅𝑛=𝑏₁𝛽,onΓ,₃∇𝑢₁+𝛽₄∇𝑢₂⋅𝑛=𝑏₂𝑢,onΓ,₁=𝑔₁𝑢,on𝜕Ω,₂=𝑔₂,on𝜕Ω, (26) on the paralle lepiped domain Ω=(𝑥_min,𝑥_max)×(𝑦_min,𝑦_max)×(𝑧_min,𝑧_max). Γ is an interface and prescribed by the zero level-set {(𝑥,𝑦,𝑧)∈Ω∣𝜙(𝑥,𝑦,𝑧)=0} of a level-set function 𝜙(𝑥,𝑦,𝑧). The unit normal vector of Γ is 𝑛=∇𝜙/|∇𝜙| pointing from Ω⁻={(𝑥,𝑦,𝑧)∈Ω∣𝜙(𝑥,𝑦,𝑧)≤0} to Ω⁺={(𝑥,𝑦,𝑧)∈Ω∣𝜙(𝑥,𝑦,𝑧)≥0}.

In all numerical experiments later, the level-set function 𝜙(𝑥,𝑦,𝑧), the coefficients 𝛽^±(𝑥,𝑦,𝑧), 𝑝^±(𝑥,𝑦,𝑧), and 𝑞^±(𝑥,𝑦,𝑧), and the solutions, 𝑢₁=𝑢⁺₁𝑥,𝑦,𝑧,inΩ⁺,𝑢₂=𝑢⁺₂𝑥,𝑦,𝑧,inΩ⁺,𝑢₁=𝑢⁻₁𝑥,𝑦,𝑧,inΩ⁻,𝑢₂=𝑢⁻₂𝑥,𝑦,𝑧,inΩ⁻, (27) are given. Hence 𝑓₁=−∇⋅𝛽₁∇𝑢₁−∇⋅𝛽₂∇𝑢₂,inΩ^±,𝑓₂=−∇⋅𝛽₃∇𝑢₁−∇⋅𝛽₄∇𝑢₂,inΩ^±,𝑎₁=𝑢⁺₁−𝑢⁻₁𝑎,onΓ,₂=𝑢⁺₂−𝑢⁻₂𝑏,onΓ,₁=𝛽⁺₁∇𝑢⁺₁+𝛽⁺₂∇𝑢⁺₂⋅𝑛−𝛽⁻₁∇𝑢⁻₁+𝛽⁻₂∇𝑢⁻₂𝑏⋅𝑛,onΓ,₂=𝛽⁺₃∇𝑢⁺₁+𝛽⁺₄∇𝑢⁺₂⋅𝑛−𝛽⁻₃∇𝑢⁻₁+𝛽⁻₄∇𝑢⁻₂⋅𝑛,onΓ. (28)𝑔 is obtained as a proper Dirichlet boundary condition, since the solutions are given.

All errors in solutions are measured in the 𝐿^∞ norm in the whole domain Ω.

We present three numerical examples to demonstrate the effectiveness of our method.

Example 1. The level-set function 𝜙(𝑥,𝑦,𝑧), the coefficients 𝛽^±₁(𝑥,𝑦,𝑧), 𝛽^±₂(𝑥,𝑦,𝑧), 𝛽^±₃(𝑥,𝑦,𝑧), and 𝛽^±₄(𝑥,𝑦,𝑧), and the solution 𝑢^±₁(𝑥,𝑦,𝑧), 𝑢^±₂(𝑥,𝑦,𝑧) in Example 1 are given as follows: 𝜙𝑥,𝑦,𝑧=minmin(𝑥−0.4)²+𝑦²+𝑧²()−0.25,𝑥+0.3²+𝑦²+𝑧²,𝑥−0.25₂+𝑦+0.5²+𝑧²,𝛽−0.25₊₁=(𝑥)𝑥,𝑦,𝑧12+sin²+𝑦cos(𝑧)²𝑥²sin𝑦+1𝑧²+9sin2𝑦²𝑥𝑧+1cos𝑦+1𝑧²,𝛽+𝑥+10₋₁=𝑥,𝑦,𝑧8𝑥²()()+7sin𝑥+𝑧+1𝑦𝑧+1𝑧+𝑥²12+cos𝑥+𝑦²0𝑥²𝑥+𝑦²5𝑦,𝛽8+cos₊₂=𝑥,𝑦,𝑧0.5𝑥²sin𝑦²(𝑧)()5+cos𝑥𝑧²(𝑧)(𝑥)3+𝑥𝑦𝑦3−𝑦cos+1sin²,𝛽⁻₂=(𝑥)𝑦𝑥,𝑦,𝑧0.2sin+0.3cos𝑦+𝑧+11.8+𝑥𝑧cos²3+sin(𝑧)²𝑦𝑦0.6𝑦+0.64+2cos0.5+sin2,𝛽0.3𝑥+0.3₊₃=𝑥,𝑦,𝑧sin𝑦+𝑧²1−cos(𝑥)3𝑦(𝑥)1.5+coscos²2+𝑥𝑦𝑦𝑧5(𝑥)𝑦+0.52.3+sinsin²,𝛽𝑥𝑧+1₋₃=𝑥,𝑦,𝑧𝑥+𝑦6(𝑧)(𝑥)+0.502+sincos²()(𝑧)2−𝑦𝑧1−sin𝑥+𝑧2+cos²sin𝑦5,𝛽+0.2cos𝑦−𝑧+1₊₄=(𝑥)𝑥,𝑦,𝑧8+sin𝑦+𝑧5+0.2sin𝑦²cos(𝑧)3+0.413−𝑥𝑦cos(𝑧)²,𝛽sin𝑦+𝑧+1𝑥+119+𝑦𝑧₋₄=𝑦𝑥𝑥,𝑦,𝑧12+cos5()(𝑧)+0.2sin𝑥+𝑧+1cos²10+𝑥²𝑥²sin(𝑥+𝑧)+1cos𝑦²,𝑢15+𝑥₊₁(𝑥)𝑦𝑧𝑥,𝑦,𝑧=sinsin2,𝑢⁻₁𝑥,𝑦,𝑧=𝑦²(𝑧),𝑢+cos₊₂𝑥,𝑦,𝑧=𝑥²𝑢+𝑦𝑧,₋₂𝑥,𝑦,𝑧=𝑦−cos(𝑥+𝑧). (29)
Table 1 shows the error on different grids. The numerical result shows second-order accuracy in the 𝐿^∞ norm for the solution. Table 2 includes the condition number, which grows with order 𝑂(ℎ⁻²), the same as the case without interface.

Table 1 

Example 1.

Table 2 

Example 1 with condition number.

Example 2. The level-set function 𝜙(𝑥,𝑦,𝑧), the coefficients 𝛽^±₁(𝑥,𝑦,𝑧), 𝛽^±₂(𝑥,𝑦,𝑧), 𝛽^±₃(𝑥,𝑦,𝑧), and 𝛽^±₄(𝑥,𝑦,𝑧), and the solution 𝑢^±₁(𝑥,𝑦,𝑧), 𝑢^±₂(𝑥,𝑦,𝑧) in Example 2 are given as follows: 𝜙𝑥,𝑦,𝑧=(𝑥−0.4)²+𝑦²+𝑧²𝛽−0.16,₊₁=(𝑥)𝑥,𝑦,𝑧12+sin²+𝑦cos(𝑧)²𝑥²sin𝑦+1𝑧²+9sin2𝑦²𝑥𝑧+1cos𝑦+1𝑧²,𝛽+𝑥+10₋₁=𝑥,𝑦,𝑧8𝑥²()()+7sin𝑥+𝑧+1𝑦𝑧+1𝑧+𝑥²12+cos𝑥+𝑦²0𝑥²𝑥+𝑦²5𝑦,𝛽8+cos₊₂=𝑥,𝑦,𝑧0.5𝑥²sin𝑦²(𝑧)()5+cos𝑥𝑧²(𝑧)(𝑥)3+𝑥𝑦𝑦3−𝑦cos+1sin²,𝛽⁻₂=(𝑥)𝑦𝑥,𝑦,𝑧0.2sin+0.3cos𝑦+𝑧+11.8+𝑥𝑧cos²3+sin(𝑧)²𝑦𝑦0.6𝑦+0.64+2cos0.5+sin2,𝛽0.3𝑥+0.3₊₃=𝑥,𝑦,𝑧sin𝑦+𝑧²(1−cos(𝑥))3𝑦(𝑥)1.5+coscos²2+𝑥𝑦𝑦𝑧5(𝑥)𝑦+0.52.3+sinsin²,𝛽𝑥𝑧+1₋₃=𝑥,𝑦,𝑧𝑥+𝑦6(𝑧)(𝑥)+0.502+sincos²()(𝑧)2−𝑦𝑧1−sin𝑥+𝑧2+cos²sin𝑦5,𝛽+0.2cos𝑦−𝑧+1₊₄=(𝑥)𝑥,𝑦,𝑧8+sin𝑦+𝑧5+0.2sin𝑦²cos(𝑧)3+0.413−𝑥𝑦cos(𝑧)²,𝛽sin𝑦+𝑧+1𝑥+119+𝑦𝑧₋₄=𝑦𝑥𝑥,𝑦,𝑧12+cos5()(𝑧)+0.2sin𝑥+𝑧+1cos²10+𝑥²𝑥²sin(𝑥+𝑧)+1cos𝑦²,𝑢15+𝑥₊₁(𝑥)𝑦𝑧𝑥,𝑦,𝑧=sinsin2,𝑢⁻₁𝑥,𝑦,𝑧=𝑦²(𝑧),𝑢+cos₊₂𝑥,𝑦,𝑧=𝑥²+𝑦²+𝑧^25/6,𝑢⁻₂𝑥,𝑦,𝑧=𝑦−cos(𝑥+𝑧). (30)
Table 3 shows the error on different grids. The numerical result shows more than 1.5th order accuracy in the 𝐿^∞ norm for the solution. The degeneracy is due to the poor regularity of the solution.

Table 3 

Example 2.

Example 3. The level-set function 𝜙(𝑥,𝑦,𝑧), the coefficients 𝛽^±₁(𝑥,𝑦,𝑧), 𝛽^±₂(𝑥,𝑦,𝑧), 𝛽^±₃(𝑥,𝑦,𝑧), and 𝛽^±₄(𝑥,𝑦,𝑧), and the solution 𝑢^±₁(𝑥,𝑦,𝑧), 𝑢^±₂(𝑥,𝑦,𝑧) in Example 3 are given as follows: 𝜙𝑥,𝑦,𝑧=minmin(𝑥−0.4)²+𝑦²+𝑧²()−0.25,𝑥+0.3²+𝑦²+𝑧²,𝑥−0.25₂+𝑦+0.5²+𝑧²,𝛽−0.25₊₁=,𝛽𝑥,𝑦,𝑧1000030003₋₁=,𝛽𝑥,𝑦,𝑧700020002₊₂=,𝛽𝑥,𝑦,𝑧004030300₋₂=,𝛽𝑥,𝑦,𝑧003020200₊₃=,𝛽𝑥,𝑦,𝑧003030400₋₃=,𝛽𝑥,𝑦,𝑧002020300₊₄=,𝛽𝑥,𝑦,𝑧3000300010₋₄=,𝑢𝑥,𝑦,𝑧200020007₊₁𝑥,𝑦,𝑧=exp𝑥²+𝑦²+𝑧²,𝑢⁻₁𝑥,𝑦,𝑧=cos(𝑧)+𝑦𝑥²,𝑢⁺₂𝑥,𝑦,𝑧=log𝑥²+𝑦²+𝑧²,𝑢⁻₂𝑥,𝑦,𝑧=𝑦²+𝑥𝑧. (31)
Table 4 shows the error on different grids. The numerical result shows second-order accuracy in the 𝐿^∞ norm for the solution. This type of examples is also presented in [22, 27] in two dimensions. The coefficient matrices for the PDE meet the assumption in Theorem 6.

Table 4 

Example 3.

Example 4. The level-set function 𝜙(𝑥,𝑦,𝑧), the coefficients 𝛽^±₁(𝑥,𝑦,𝑧), 𝛽^±₂(𝑥,𝑦,𝑧), 𝛽^±₃(𝑥,𝑦,𝑧), and 𝛽^±₄(𝑥,𝑦,𝑧), and the solution 𝑢^±₁(𝑥,𝑦,𝑧), 𝑢^±₂(𝑥,𝑦,𝑧) in Example 4 are given as follows:𝜙𝑥,𝑦,𝑧=min𝑥²+𝑦²+(𝑧+0.5)²−0.25,minmin(𝑥−0.4)²+𝑦²+𝑧²()−0.25,𝑥+0.3²+𝑦²+𝑧²,𝑥−0.25₂+𝑦+0.5²+𝑧²,𝛽−0.25₊₁=𝑥𝑥,𝑦,𝑧²+𝑧²+2sin𝑥+𝑦+𝑧²+5000sin𝑥+𝑦+𝑧²+2000sin𝑥+𝑦+𝑧²,𝛽+2₋₁=𝑦𝑥,𝑦,𝑧²+2+2cos𝑦𝑧²+2𝑥²000cos𝑦𝑧²+2𝑥²000cos𝑦𝑧²+2𝑥²,𝛽⁺₂𝑥,𝑦,𝑧=00𝑥²+𝑧²+30sin𝑥+𝑦+𝑧²+20sin𝑥+𝑦+𝑧²,𝛽+200₋₂𝑥,𝑦,𝑧=00𝑦²+20cos𝑦𝑧²+2𝑥²0cos𝑦𝑧²+2𝑥²,𝛽00₊₃𝑥,𝑦,𝑧=00sin𝑥+𝑦+𝑧²+20sin𝑥+𝑦+𝑧²𝑥+20₂+𝑧²,𝛽+300₋₃𝑥,𝑦,𝑧=00cos𝑦𝑧²+2𝑥²0cos𝑦𝑧²+2𝑥²0𝑦²,𝛽+200₊₄𝑥,𝑦,𝑧=sin𝑥+𝑦+𝑧²+2000sin𝑥+𝑦+𝑧²+2000𝑥²+𝑧²+2sin𝑥+𝑦+𝑧²,𝛽+5₋₄𝑥,𝑦,𝑧=cos𝑦𝑧²+2𝑥²000cos𝑦𝑧²+2𝑥²000𝑦²+2+2cos𝑦𝑧²+2𝑥²,𝑢⁺₁𝑥,𝑦,𝑧=exp𝑥²+𝑦²+𝑧²,𝑢⁻₁𝑥,𝑦,𝑧=cos(𝑧)+𝑦𝑥²,𝑢⁺₂𝑥,𝑦,𝑧=log𝑥²+𝑦²+𝑧²,𝑢⁻₂𝑥,𝑦,𝑧=𝑦²+𝑥𝑧. (32)

Table 5 shows the error on different grids. The numerical result shows second-order accuracy in the 𝐿^∞ norm for the solution. Compared with Example 3, the difference is that we have variable coefficient matrices of the same structure as those constant coefficient matrices in Example 3.

Table 5 

Example 4.

Example 5. The level-set function 𝜙(𝑥,𝑦,𝑧), the coefficients 𝛽^±₁(𝑥,𝑦,𝑧), 𝛽^±₂(𝑥,𝑦,𝑧), 𝛽^±₃(𝑥,𝑦,𝑧), and 𝛽^±₄(𝑥,𝑦,𝑧), and the solution 𝑢^±₁(𝑥,𝑦,𝑧), 𝑢^±₂(𝑥,𝑦,𝑧) in Example 5 are given as follows: 𝜙𝑥,𝑦,𝑧=min𝑥²+𝑦²+(𝑧+0.5)²()−0.25,minmin𝑥−0.4²+𝑦²+𝑧²()−0.25,𝑥+0.3²+𝑦²+𝑧²,𝑥−0.25₂+𝑦+0.5²+𝑧²,𝛽−0.25₊₁=(𝑥)𝑥,𝑦,𝑧12+sin²+𝑦cos(𝑧)²𝑥²sin𝑦+1𝑧²+9sin2𝑦²𝑥𝑧+1cos𝑦+1𝑧²,𝛽+𝑥+10₋₁=𝑥,𝑦,𝑧8𝑥²()()+7sin𝑥+𝑧+1𝑦𝑧+1𝑧+𝑥²12+cos𝑥+𝑦²0𝑥²𝑥+𝑦²5𝑦,𝛽8+cos₊₂=𝑥,𝑦,𝑧0.5𝑥²sin𝑦²(𝑧)()5+cos𝑥𝑧²3+𝑥𝑦𝑦²3−𝑦cos(𝑧)+1sin(𝑥)²,𝛽⁻₂=(𝑥)𝑦𝑥,𝑦,𝑧0.2sin+0.3cos𝑦+𝑧+11.8+𝑥𝑧cos²3+sin(𝑧)²𝑦𝑦0.6𝑦+0.63+2cos0.5+sin2,𝛽0.3𝑥+0.3₊₃=𝑥,𝑦,𝑧sin𝑦+𝑧²(1−cos(𝑥))3𝑦(𝑥)1.5+coscos²2+𝑥𝑦𝑦𝑧5(𝑥)𝑦+12.3+sinsin²,𝛽𝑥𝑧+1₋₃=𝑥,𝑦,𝑧𝑥+𝑦6(𝑧)(𝑥)+0.502+sincos²()(𝑧)2−𝑦𝑧1−sin𝑥+𝑧2+cos²sin𝑦5,𝛽+0.2cos𝑦−𝑧+1₊₄=(𝑥)𝑥,𝑦,𝑧8+sin𝑦+𝑧5+1sin𝑦²cos(𝑧)3+0.413−𝑥𝑦cos(𝑧)²,𝛽sin𝑦+𝑧+1𝑥+119+𝑦𝑧₋₄=𝑦𝑥𝑥,𝑦,𝑧12+cos5()(𝑧)+0.2sin𝑥+𝑧+1cos²10+𝑥²𝑥²sin(𝑥+𝑧)+1cos𝑦²,𝑢15+𝑥₊₁𝑥,𝑦,𝑧=exp𝑥²+𝑦²+𝑧²,𝑢⁻₁𝑥,𝑦,𝑧=cos(𝑧)+𝑥³,𝑢⁺₂𝑥,𝑦,𝑧=log𝑥²+𝑦²+𝑧²,𝑢⁻₂𝑥,𝑦,𝑧=𝑦³+𝑥𝑧. (33)
Table 6 shows the error on different grids. The numerical result shows second-order accuracy in the 𝐿^∞ norm for the solution. Compared with Example 4, the difference is that we used general variable coefficient matrices without the special structure in Example 4 to demonstrate that our method can handle more general problems.

Table 6 

Example 5.

Example 6. The level-set function 𝜙(𝑥,𝑦,𝑧), the coefficients 𝛽₁(𝑥,𝑦,𝑧), 𝛽₂(𝑥,𝑦,𝑧), 𝛽₃(𝑥,𝑦,𝑧), and 𝛽₄(𝑥,𝑦,𝑧), and the solution 𝑢^±₁(𝑥,𝑦,𝑧), 𝑢^±₂(𝑥,𝑦,𝑧) in this example are given as follows: 𝜙𝑥,𝑦,𝑧=𝑥²+𝑦²+𝑧²𝛽−0.25,₊₁=,𝛽𝑥,𝑦,𝑧100001000010₊₂=,𝛽𝑥,𝑦,𝑧100010001₊₃=,𝛽𝑥,𝑦,𝑧200020002₊₄=,𝑢𝑥,𝑦,𝑧300003000030₊₁𝑥,𝑦,𝑧=exp𝑥²+𝑦²+𝑧²,𝑢⁻₁𝑥,𝑦,𝑧=cos(𝑧)+𝑦𝑥²,𝑢⁺₂𝑥,𝑦,𝑧=log𝑥²+𝑦²+𝑧²,𝑢⁻₂𝑥,𝑦,𝑧=𝑦²+𝑥𝑧. (34)

When it is 11×11×11 grids, we get the result of Table 7. When it is 22×22×22 grids, we get the result of Table 8. We plot the condition number versus the ratio between 𝛽⁻ and 𝛽⁺ for 11-by-11-by-11 grid in Figure 4 and 22-by-22-by-22 grid in Figure 5. The correlation coefficients for 11-by-11-by-11, 22-by-22-by-22 grids are 0.999990 and 0.999987, which are clearly linear relations.

Table 7 

Example 6. 11×11×11 grids.

Table 8 

Example 6. 22×22×22 grids.

Figure 4 

Condition number versus the ratio between 𝛽− and 𝛽+ for 11-by-11-by-11 grid.

Figure 5 

Condition number versus the ratio between 𝛽− and 𝛽+ for 22-by-22-by-22 grid.

5. Conclusion

In this paper, we modified the method in [27] to solve the 3D elasticity problem with sharp-edged interfaces. Due to the more complicated geometry in 3D, this is not a trivial extension. We proved that the matrix for the linear system generated by our method is positive definite (but not symmetric) under certain assumption. Through numerical experiments, our method achieved close to second-order accuracy in the 𝐿^∞ norm for piecewise smooth solutions, and we can handle the difficulties of sharp-edged interfaces and singular solutions. The condition number of the large sparse linear system grows with order 𝑂(ℎ⁻²), the same as the case without interface. The condition number also grows only linearly with the ratio between 𝛽⁻ and 𝛽⁺.

Acknowledgment

L. Wang’s research is supported by the Science Foundation of China University of Petroleum, Beijing (no. YJRC-2013-48). S. Hou’s research is supported by NSF grant DMS-1317994.