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

Authors: Liqun Wang; Songming Hou; Liwei Shi
Journal: International Journal of Partial Differential Equations
 (2013)

---

## 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 theL∞ norm for piecewise smooth solutions. More than 1.5th order accuracy is observed for solution with singularity (second derivative blows up).


---

## Body

## 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Ω⊂Rd. Let Γ be an interface of codimension d-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)-∇·(β1(x)∇u1(x))-∇·(β2(x)∇u2(x))=f1(x),-∇·(β3(x)∇u1(x))-∇·(β4(x)∇u2(x))=f2(x),x∈Ω∖Γ,
in which x=(x1,…,xd) denotes the spatial variables and ∇ is the gradient operator. The coefficient β(x) is assumed to be a d×d 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 fi(x), i=1,2 is assumed to lie in L2(Ω).Given functionsa and b along the interface Γ, we prescribe the jump conditions:
(2)[u1]Γ(x)≡u1+(x)-u1-(x)=a1(x),[u2]Γ(x)≡u2+(x)-u2-(x)=a2(x),n·(β1+(x)∇u1+(x)+β2+(x)∇u2+(x))-n·(β1-(x)∇u1-(x)+β2-(x)∇u2-(x))=b1(x),n·(β3+(x)∇u1+(x)+β4+(x)∇u2+(x))-n·(β3-(x)∇u1-(x)+β4-(x)∇u2-(x))=b2(x).The “±” superscripts refer to limits taken from within the subdomains Ω±.Finally, we prescribe boundary conditions(3)u1(x)=g1(x),u2(x)=g2(x),x∈∂Ω,
for a given function g on the boundary ∂Ω.The setup of the problem is illustrated in Figure1.Figure 1
Setup of the problem.We consider the three-dimensional elasticity interface problem (d=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 L∞ 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 βi to be matrices. The condition number grows with O(h-2), 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: Section2 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 H1(Ω). For H01(Ω), instead of the usual inner product, we choose one which is better suited to our problem:
(4)B[u,v]={∫Ω+‍(β1∇u1·∇v1+β2∇u2·∇v1)+∫Ω-‍(β1∇u1·∇v1+β2∇u2·∇v1),∫Ω+‍(β3∇u1·∇v2+β4∇u2·∇v2)+∫Ω-‍(β3∇u1·∇v2+β4∇u2·∇v2).LetΓ be any closed Lipschitz continuous hypersurface of dimension d-1 in Ω¯, where the overline denotes the closure of a set. Let R∂Ω- denote the restriction operator from H1(Ω) to L2(∂Ω-). This restriction operator R∂Ω- is well defined and bounded, because it is closed Lipschitz continuous (see Theorem  2.4.2 in [32]) and C1(Ω¯) is dense in H1(Ω). Throughout this section, we shall always assume that our interface data a and b are the restrictions of functions a~ and b~∈H1(Ω) on ∂Ω- and then limited on Γ, respectively. That is on Γ,
(5)a=R∂Ω-(a~),b=R∂Ω-(b~).
We shall always assume that our boundary data g can be obtained: assume that there exists a function c~∈H1(Ω) so that g is given as, on ∂Ω,
(6)g={R∂Ω(c~-a~),on∂Ω⋂‍∂Ω-,R∂Ω(c~),on∂Ω∖∂Ω-.Equation (6) could be thought as a compatibility condition between a and g. 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)H(a,c)={u:u-c+aχ(Ω-¯)∈H01(Ω)}.Ifu∈H(a,c), then [u]Γ=a and u|∂Ω=g. Note that H01(Ω) can be identified with H(0,0).Definition 1.
A functionu∈H(a,c) is a weak solution of (1)–(3), if v=u-c+aχ(Ω-¯)∈H01(Ω) satisfies
(8)B[v,ψ]=F(ψ),
for all ψ∈H01(Ω), where
(9)B[v,ψ]=∫Ω+‍β∇v·∇ψ+∫Ω-‍β∇v·∇ψ,F(ψ)=∫Ω‍fψ+∫Ω‍β∇c·∇ψ+∫Ω-‍β∇a·∇ψ+∫Γ‍bψ,
or equivalently.Definition 2.
A functionu∈H(a,c) is a weak solution of (1)–(3), if u satisfies, for all ψ∈H01(Ω),
(10)∫Ω+‍(β1∇u1·∇ψ1+β2∇u2·∇ψ1)+∫Ω-‍(β1∇u1·∇ψ1+β2∇u2·∇ψ1)=∫Ω‍f1ψ1+∫Γ‍b1ψ1,∫Ω+‍(β3∇u1·∇ψ2+β4∇u2·∇ψ2)+∫Ω-‍(β3∇u1·∇ψ2+β4∇u2·∇ψ2)=∫Ω‍f2ψ2+∫Γ‍b2ψ2.
A classical solution of (1)–(3), u|Ω±∈C2(Ω±¯), 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.
Iff∈L2(Ω), and a, b, and c∈H1(Ω), then there exists a unique weak solution of (1)–(3) in H(a,c).Proof.
See Theorem 2.1 in [31].


## 3. Numerical Method

We define(11)u≡[u1u2],f≡[f1f2],g≡[g1g2],a≡[a1a2],b≡[b1b2],β≡[β1β2β3β4],
and choose test function
(12)ψ=[ψ10]or[0ψ2]
and redefine the gradient and divergence operator
(13)∇≡[∂∂x0∂∂y0∂∂z00∂∂x0∂∂y0∂∂z],∇·≡[∂∂x∂∂y∂∂z000000∂∂x∂∂y∂∂z].Then (1) can be written as
(14)-∇·(β(x)∇u(x))=f(x),x∈Ω∖Γ.
The jump condition equation (2) can be reformed as
(15)[u]Γ(x)≡u+(x)-u-(x)=a(x),n·(β+(x)∇u+(x))-n·(β-(x)∇u-(x))=b(x),
and boundary condition is
(16)u(x)=g(x)x∈∂Ω.For simplicity, we will discuss the following properties under the form of (14), (15), and (16).For ease of discussion, Section3, and for accuracy testing in Section 4, we assume that a, b, and c are smooth on the closure of Ω, β and f 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 n=∇ϕ/|∇ϕ| 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Ω=(xmin,xmax)×(ymin,ymax)×(zmin,zmax) in 3D space, and β is a 3×3 matrix that is uniformly elliptic in each subdomain. 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 (xi,yj,zk)=(xmin+iΔx,ymin+jΔy,zmin+kΔz) for i=0,…,I, j=0,…,J and k=0,…,K. Each (xi,yj,zk) is called a grid point. For the case i=0,I, j=0,J, or k=0,K, a grid point is called a boundary point, otherwise it is called an interior point. The grid size is defined as h=max(Δx,Δy,Δz)>0.Two sets of grid functions are needed, and they are denoted by(17)H±1,h={ωh=(ωi,j,k):0≤i≤I,0≤j≤J,0≤k≤K},H01,h={ωh=(ωi,j,k)∈H±1,h:ωi,j,k=0ifi=0,H01,h=vvIorj=0,Jork=0,K}.We cut every cube cell region[xi,xi+1]×[yj,yj+1]×[zk,zk+1] into six tetrahedron regions. Collecting all those tetrahedron regions, we obtain a uniform tetrahedralization Th:⋃L∈Th‍L; see Figure 2.Figure 2
Cubic mesh and the tetrahedralization.Figure 3
Interface tetrahedron, left: Case 1, right: Case 2.Ifϕ(xi,yj,zk)≤0, we count the grid point (xi,yj,zk) 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 cellL an interface cell if its vertices belong to different subdomains. In the interface cell, we write L=L+⋃L-. L+ and L- are separated by a plane segment, denoted by ΓLh. There are two kinds of plane segment; see Figure 3. The vertices of the plane segment ΓLh are located on the interface Γ, and their locations can be calculated from the linear interpolations of the discrete level-set functions ϕh=ϕ(xi,yj,zk). The vertices of L+ are located in Ω+⋃Γ and the vertices of L- are located in Ω-⋃Γ. L+ and L-are approximations of the regions of L⋂Ω+ and L⋂Ω-, respectively. We call a cell L a regular cell if all its vertices belong to the same subdomain, either Ω+ or Ω-. For a regular cell, we also write L=L+⋃L-, where L-={} (empty set) if all vertices of L are in Ω+ and L+={} if all vertices of L are in Ω-. Clearly ΓLh={} in a regular cell, and L+ and L- are approximations of the regions L⋂Ω+ and L⋂Ω-, respectively. We use |L+| and |L-| to represent the volumes of L+ and L-, respectively.Two 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 tetrahedron 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-. The location of its discontinuity in the interface cell is the plane segment ΓLh. Note that the vertices of the plane segment are located on the interface Γ, and hence the interface condition [u]=a is enforced exactly at the vertices of ΓLh. In each interface cell, the interface condition [β∇u·n]=b is enforced with the value b at the barycenter of ΓLh 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 Ha,c1,h, which should be thought as an approximation of the solution class H(a,c) (see (7)) plus the restriction of [β∇u·n]=b. 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 Uh 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 functionuh∈H1,h such that uh=gh on boundary points and so that for all ψh∈H01,h, we have
(18)∑L∈Th‍(∫L+‍β∇Uh(uh)·∇Th(ψh)+∫L-‍β∇Uh(uh)·∇Th(ψh))=∑L∈Th‍(∫L+‍fTh(ψh)+∫L-‍fTh(ψh)+∫ΓLh‍bTh(ψh)).
Note that u=g on the boundary is the same as u-c+aχ(Ω-¯)=0 on the boundary.Since our solution bases and test function bases are different, the matrixA 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, A is positive definite numerically.Theorem 5.
Ifβ is continuous and positive definite, then the matrix A for the linear system generated by Method 4 is positive definite.Proof.
For vectorc∈R2n we have
(19)cTAc=B[∑i=12n‍ciui,∑i=12n‍ciψi].
Here we can choose specific jump conditions a1=a2=b1=b2=0 since changing ai, bi will not affect the matrix A. In this case we have
(20)∑i=12n‍ciui=∑i=12n‍ciψi=w,
since
(21)B[w,w]=∫Ω‍(β(x)∇w(x))T∇w(x)dx.
So if β is positive definite, then ∇w(x)TβT∇w(x)>0, for all x∈Ω. Therefore B[w,w]>0, which implies cTAc>0. Thus, A 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 A generated by our method is still positive definite.Theorem 6.
Ifλ>0, μ>0 and β1=[λ+2μ000μ000μ], β2=[00λ0μ0μ00], β3=[00μ0μ0λ00], and β4=[μ000μ000λ+2μ], then the matrix A for the linear system generated by Method 4 is positive definite.Proof.
Suppose for a contradiction thatA is not positive definite. Then there is a vector c∈R2n and c≠0 such that cTAc≤0. Let
(22)w=[w1w2]=∑i=12n‍ciψi=∑i=12n‍ciui.
Let
(23)Dw=[∂w1∂x∂w1∂y∂w1∂z∂w2∂x∂w2∂y∂w2∂z],β1234=[λ+2μ0000λ0μ00μ000μμ0000μμ000μ00μ0λ0000λ+2μ],Gw=[∂w1∂x∂w1∂y∂w1∂z∂w2∂x∂w2∂y∂w2∂z].
Then
(24)B[w,w]≤0⟹∫Ω‍(β∇w(x→))T∇w(x)dx→≤0⟹∫Ω‍Dwβ1234Gwdx→≤0.
Since for alla=[a1,a2,a3,a4,a5,a6]T∈R6,
(25)aTβa=(a1+a6)2λ+2(a12+a62)μ+(a2+a5)2μ+(a3+a4)2μ≥0.
SoaTβa=0 if and only if a1=a6=0, a2=-a5, and a3=-a4. Then (∂w1/∂x)(x→)=a1=0, for all x→∈Ω.
However,w1=∑i=1n‍ciψi1 implies that ∂w1/∂x=∑i=1n‍ci(∂ψi1/∂x).
Sincec=[c1,c2,…,c2n]T≠0, without loss of generality, we assume that c1≠0. If we choose point x→∈Ω such that ∂ψ11(x→)/∂x≠0 and ∂ψi1(x→)/∂x=0, i=2,3,…,n, then ∑i=1n‍ci(∂ψi1/∂x)≠0, a contradiction.
ThereforecTAc>0 for all c≠0; thus, A is positive definite.Remark 7.
A positive definite matrix has positive determinant and is therefore invertible. It also has anLDMT factorization where D=diag(di) and di>0, and L,M are lower triangular. The linear system Ax=b can be solved efficiently.


## 4. Numerical Experiments

Consider the problem(26)-∇·(β1∇u1)-∇·(β2∇u2)=f1,inΩ±,-∇·(β3∇u1)-∇·(β4∇u2)=f2,inΩ±,[u1]=a1,onΓ,[u2]=a2,onΓ,[(β1∇u1+β2∇u2)·n]=b1,onΓ,[(β3∇u1+β4∇u2)·n]=b2,onΓ,u1=g1,on∂Ω,u2=g2,on∂Ω,
on the paralle lepiped domain Ω=(xmin,xmax)×(ymin,ymax)×(zmin,zmax). Γ is an interface and prescribed by the zero level-set {(x,y,z)∈Ω∣ϕ(x,y,z)=0} of a level-set function ϕ(x,y,z). The unit normal vector of Γ is n=∇ϕ/|∇ϕ| pointing from Ω-={(x,y,z)∈Ω∣ϕ(x,y,z)≤0} to Ω+={(x,y,z)∈Ω∣ϕ(x,y,z)≥0}.In all numerical experiments later, the level-set functionϕ(x,y,z), the coefficients β±(x,y,z), p±(x,y,z), and q±(x,y,z), and the solutions,
(27)u1=u1+(x,y,z),inΩ+,u2=u2+(x,y,z),inΩ+,u1=u1-(x,y,z),inΩ-,u2=u2-(x,y,z),inΩ-,
are given. Hence
(28)f1=-∇·(β1∇u1)-∇·(β2∇u2),inΩ±,f2=-∇·(β3∇u1)-∇·(β4∇u2),inΩ±,a1=u1+-u1-,onΓ,a2=u2+-u2-,onΓ,b1=(β1+∇u1++β2+∇u2+)·n-(β1-∇u1-+β2-∇u2-)·n,b1=bbbbbbbbbbbb(β1+∇u1++β2+∇u2+)·n-onΓ,b2=(β3+∇u1++β4+∇u2+)·n-(β3-∇u1-+β4-∇u2-)·n,b2=bbbbbbbbbbbb(β3+∇u1++β4+∇u2+)·n-onΓ.g is obtained as a proper Dirichlet boundary condition, since the solutions are given.All errors in solutions are measured in theL∞ norm in the whole domain Ω.We present three numerical examples to demonstrate the effectiveness of our method.Example 1.
The level-set functionϕ(x,y,z), the coefficients β1±(x,y,z), β2±(x,y,z), β3±(x,y,z), and β4±(x,y,z), and the solution u1±(x,y,z), u2±(x,y,z) in Example 1 are given as follows:
(29)ϕ(x,y,z)=min(min((x-0.4)2+y2+z2-0.25,y2+gz2f-025(x+0.3)2+y2+z2-0.25),fufz2-0x2+(y+0.5)2+z2-0.25),β1+(x,y,z)=(12+sin(x)2+ycos(z)2x2sin(y)+1z2+9sin(2y)2xz+1cos(y)+1z2+x+10),β1-(x,y,z)=(8x2+7sin(x+z)+1yz+1(z+x)212+cos(x+y)20x2(x+y)258+cos(y)),β2+(x,y,z)=(0.5x2sin(y)25+cos(z)(xz)23+xyy3-ycos(z)+1sin(x)2),β2-(x,y,z)=(0.2sin(x)+0.3cos(y+z)+11.8+xz(cos(y))23+sin(z)20.6y+0.64+2cos(y)0.5+sin(y)20.3x+0.3),β3+(x,y,z)=(sin(y+z)21-cos(x)31.5+cos(y)cos(x)22+xyyz5+0.52.3+sin(x)sin(y)2xz+1),β3-(x,y,z)=((x+y)6+0.502+sin(z)cos(x)22-yz1-sin(x+z)2+cos(z)2sin(y)5+0.2cos(y-z)+1),β4+(x,y,z)=(8+sin(x)(y+z)5+0.2sin(y)2cos(z)3+0.413-xycos(z)2sin(y+z)+1x+119+yz),β4-(x,y,z)=(12+cos(y)x5+0.2sin(x+z)+1cos(z)210+x2x2sin(x+z)+1cos(y)215+x),u1+(x,y,z)=sin(x)sin(y)z2,u1-(x,y,z)=y2+cos(z),u2+(x,y,z)=x2+yz,u2-(x,y,z)=y-cos(x+z).
Table1 shows the error on different grids. The numerical result shows second-order accuracy in the L∞ norm for the solution. Table 2 includes the condition number, which grows with order O(h-2), the same as the case without interface.Table 1
Example1.
n
x
×
n
y
×
n
z
Err inU
Order
11
×
11
×
11
0.00664
22
×
22
×
22
0.00174
1.93
44
×
44
×
44
0.00046
1.90
88
×
88
×
88
0.00012
1.90Table 2
Example1 with condition number.
n
x
×
n
y
×
n
z
Err inU
Order
Condition number
8
×
8
×
8
0.00944
80.49853
16
×
16
×
16
0.00302
1.64
3.35416
e
+
002
32
×
32
×
32
0.00080
1.92
1.56186
e
+
003Example 2.
The level-set functionϕ(x,y,z), the coefficients β1±(x,y,z), β2±(x,y,z), β3±(x,y,z), and β4±(x,y,z), and the solution u1±(x,y,z), u2±(x,y,z) in Example 2 are given as follows:
(30)ϕ(x,y,z)=(x-0.4)2+y2+z2-0.16,β1+(x,y,z)=(12+sin(x)2+ycos(z)2x2sin(y)+1z2+9sin(2y)2xz+1cos(y)+1z2+x+10),β1-(x,y,z)=(8x2+7sin(x+z)+1yz+1(z+x)212+cos(x+y)20x2(x+y)258+cos(y)),β2+(x,y,z)=(0.5x2sin(y)25+cos(z)(xz)23+xyy3-ycos(z)+1sin(x)2),β2-(x,y,z)=(0.2sin(x)+0.3cos(y+z)+11.8+xz(cos(y))23+sin(z)20.6y+0.64+2cos(y)0.5+sin(y)20.3x+0.3),β3+(x,y,z)=(sin(y+z)2(1-cos(x))31.5+cos(y)cos(x)22+xyyz5+0.52.3+sin(x)sin(y)2xz+1),β3-(x,y,z)=((x+y)6+0.502+sin(z)cos(x)22-yz1-sin(x+z)2+cos(z)2sin(y)5+0.2cos(y-z)+1),β4+(x,y,z)=(8+sin(x)(y+z)5+0.2sin(y)2cos(z)3+0.413-xycos(z)2sin(y+z)+1x+119+yz),β4-(x,y,z)=(12+cos(y)x5+0.2sin(x+z)+1cos(z)210+x2x2sin(x+z)+1cos(y)215+x),u1+(x,y,z)=sin(x)sin(y)z2,u1-(x,y,z)=y2+cos(z),u2+(x,y,z)=(x2+y2+z2)5/6,u2-(x,y,z)=y-cos(x+z).
Table3 shows the error on different grids. The numerical result shows more than 1.5th order accuracy in the L∞ norm for the solution. The degeneracy is due to the poor regularity of the solution.Table 3
Example2.
n
x
×
n
y
×
n
z
Err inU
Order
11
×
11
×
11
0.00641
22
×
22
×
22
0.00210
1.61
44
×
44
×
44
0.00069
1.61
88
×
88
×
88
0.00022
1.63Example 3.
The level-set functionϕ(x,y,z), the coefficients β1±(x,y,z), β2±(x,y,z), β3±(x,y,z), and β4±(x,y,z), and the solution u1±(x,y,z), u2±(x,y,z) in Example 3 are given as follows:
(31)ϕ(x,y,z)=min(min((x-0.4)2+y2+z2-0.25,bbbbbbbbbbbbb(x+0.3)2+y2+z2-0.25),x2+(y+0.5)2+z2-0.25),β1+(x,y,z)=(1000030003),β1-(x,y,z)=(700020002),β2+(x,y,z)=(004030300),β2-(x,y,z)=(003020200),β3+(x,y,z)=(003030400),β3-(x,y,z)=(002020300),β4+(x,y,z)=(3000300010),β4-(x,y,z)=(200020007),u1+(x,y,z)=exp(x2+y2+z2),u1-(x,y,z)=cos(z)+yx2,u2+(x,y,z)=log(x2+y2+z2),u2-(x,y,z)=y2+xz.
Table4 shows the error on different grids. The numerical result shows second-order accuracy in the L∞ 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
Example3.
n
x
×
n
y
×
n
z
Err inU
Order
11
×
11
×
11
0.11372
22
×
22
×
22
0.03207
1.83
44
×
44
×
44
0.00839
1.93
88
×
88
×
88
0.00215
1.96Example 4.
The level-set functionϕ(x,y,z), the coefficients β1±(x,y,z), β2±(x,y,z), β3±(x,y,z), and β4±(x,y,z), and the solution u1±(x,y,z), u2±(x,y,z) in Example 4 are given as follows:(32)ϕ(x,y,z)=min(x2+y2+(z+0.5)2-0.25,min(min((x-0.4)2+y2+z2-0.25,bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbvvbb(x+0.3)2+y2+z2-0.25),bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbx2+(y+0.5)2+z2-0.25)),β1+(x,y,z)=(x2+z2+2(sin(x+y+z)2+5)000sin(x+y+z)2+2000sin(x+y+z)2+2),β1-(x,y,z)=(y2+2+2(cos(yz)2+2x2)000cos(yz)2+2x2000cos(yz)2+2x2),β2+(x,y,z)=(00x2+z2+30sin(x+y+z)2+20sin(x+y+z)2+200),β2-(x,y,z)=(00y2+20cos(yz)2+2x20cos(yz)2+2x200),β3+(x,y,z)=(00sin(x+y+z)2+20sin(x+y+z)2+20x2+z2+300),β3-(x,y,z)=(00cos(yz)2+2x20cos(yz)2+2x20y2+200),β4+(x,y,z)=(sin(x+y+z)2+2000sin(x+y+z)2+2000x2+z2+2(sin(x+y+z)2+5)),β4-(x,y,z)=(cos(yz)2+2x2000cos(yz)2+2x2000y2+2+2(cos(yz)2+2x2)),u1+(x,y,z)=exp(x2+y2+z2),u1-(x,y,z)=cos(z)+yx2,u2+(x,y,z)=log(x2+y2+z2),u2-(x,y,z)=y2+xz.Table5 shows the error on different grids. The numerical result shows second-order accuracy in the L∞ 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
Example4.
n
x
×
n
y
×
n
z
Err inU
Order
11
×
11
×
11
0.13994
22
×
22
×
22
0.04011
1.80
44
×
44
×
44
0.01036
1.95
88
×
88
×
88
0.00263
1.98Example 5.
The level-set functionϕ(x,y,z), the coefficients β1±(x,y,z), β2±(x,y,z), β3±(x,y,z), and β4±(x,y,z), and the solution u1±(x,y,z), u2±(x,y,z) in Example 5 are given as follows:
(33)ϕ(x,y,z)=min(x2+y2+(z+0.5)2-0.25,bbbbbbbmin(min((x-0.4)2+y2+z2-0.25,bbbbbbbbbbbbbbbbb(x+0.3)2+y2+z2-0.25),bbbbbbbbbbbbbx2+(y+0.5)2+z2-0.25)),β1+(x,y,z)=(12+sin(x)2+ycos(z)2x2sin(y)+1z2+9sin(2y)2xz+1cos(y)+1z2+x+10),β1-(x,y,z)=(8x2+7sin(x+z)+1yz+1(z+x)212+cos(x+y)20x2(x+y)258+cos(y)),β2+(x,y,z)=(0.5x2sin(y)25+cos(z)(xz)23+xyy23-ycos(z)+1sin(x)2),β2-(x,y,z)=(0.2sin(x)+0.3cos(y+z)+11.8+xz(cos(y))23+sin(z)20.6y+0.63+2cos(y)0.5+sin(y)20.3x+0.3),β3+(x,y,z)=(sin(y+z)2(1-cos(x))31.5+cos(y)cos(x)22+xyyz5+12.3+sin(x)sin(y)2xz+1),β3-(x,y,z)=((x+y)6+0.502+sin(z)cos(x)22-yz1-sin(x+z)2+cos(z)2sin(y)5+0.2cos(y-z)+1),β4+(x,y,z)=(8+sin(x)(y+z)5+1sin(y)2cos(z)3+0.413-xycos(z)2sin(y+z)+1x+119+yz),β4-(x,y,z)=(12+cos(y)x5+0.2sin(x+z)+1cos(z)210+x2x2sin(x+z)+1cos(y)215+x),u1+(x,y,z)=exp(x2+y2+z2),u1-(x,y,z)=cos(z)+x3,u2+(x,y,z)=log(x2+y2+z2),u2-(x,y,z)=y3+xz.
Table6 shows the error on different grids. The numerical result shows second-order accuracy in the L∞ 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
Example5.
n
x
×
n
y
×
n
z
Err inU
Order
11
×
11
×
11
0.09240
22
×
22
×
22
0.02498
1.89
44
×
44
×
44
0.00635
1.96
88
×
88
×
88
0.00160
1.99Example 6.
The level-set functionϕ(x,y,z), the coefficients β1(x,y,z), β2(x,y,z), β3(x,y,z), and β4(x,y,z), and the solution u1±(x,y,z), u2±(x,y,z) in this example are given as follows:
(34)ϕ(x,y,z)=x2+y2+z2-0.25,β1+(x,y,z)=(100001000010),β2+(x,y,z)=(100010001),β3+(x,y,z)=(200020002),β4+(x,y,z)=(300003000030),u1+(x,y,z)=exp(x2+y2+z2),u1-(x,y,z)=cos(z)+yx2,u2+(x,y,z)=log(x2+y2+z2),u2-(x,y,z)=y2+xz.When it is11×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
Example6. 11×11×11 grids.
β
-
=
β
+
×
n
Err inU
Condition number
β
-
=
β
+
×
1
0.07981
2.63433
e
+
002
β
-
=
β
+
×
50
0.20388
1.31957
e
+
004
β
-
=
β
+
×
100
0.39149
2.67980
e
+
004
β
-
=
β
+
×
150
0.54457
4.04338
e
+
004
β
-
=
β
+
×
200
0.67088
5.40725
e
+
004
β
-
=
β
+
×
250
0.77679
6.77095
e
+
004
β
-
=
β
+
×
300
0.86828
8.13431
e
+
004Table 8
Example6. 22×22×22 grids.
β
-
=
β
+
×
n
Err inU
Condition number
β
-
=
β
+
×
1
0.02173
1.08341
e
+
003
β
-
=
β
+
×
50
0.08743
5.29426
e
+
004
β
-
=
β
+
×
100
0.15283
1.07574
e
+
005
β
-
=
β
+
×
150
0.21488
1.62334
e
+
005
β
-
=
β
+
×
200
0.27731
2.17177
e
+
005
β
-
=
β
+
×
250
0.33511
2.72071
e
+
005
β
-
=
β
+
×
300
0.38877
3.26954
e
+
005Figure 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 L∞ 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 O(h-2), the same as the case without interface. The condition number also grows only linearly with the ratio between β- and β+.



---

*Source: 476873-2013-07-21.xml*