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

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 equationswith sharp-edged interfaces in three dimensions.Themain 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 L∞ norm for piecewise smooth solutions. More than 1.5th order accuracy is observed for solution with singularity (second derivative blows up).


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

2
International Journal of Partial Differential Equations Given functions and along the interface Γ, we prescribe the jump conditions: The "±" superscripts refer to limits taken from within the subdomains Ω ± .
Finally, we prescribe boundary conditions for a given function on the boundary Ω.
The setup of the problem is illustrated in Figure 1.
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-bytwo 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-bythree 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][6][7][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 twophase 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. Secondorder 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][20][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 immersedinterface 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][28][29][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 International Journal of Partial Differential Equations 3 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 (ℎ −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: 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.

The Weak Formulation
We modify the weak formulation in [28,30,31]. We are going to use the usual Sobolev spaces 1 (Ω). For 1 0 (Ω), instead of the usual inner product, we choose one which is better suited to our problem: 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 1 (Ω) to 2 ( Ω − ). This restriction operator Ω − is well defined and bounded, because it is closed Lipschitz continuous (see Theorem 2.4.2 in [32]) and 1 (Ω) is dense in 1 (Ω). Throughout this section, we shall always assume that our interface data and are the restrictions of functionsã nd̃∈ 1 (Ω) on Ω − and then limited on Γ, respectively. That is on Γ, We shall always assume that our boundary data can be obtained: assume that there exists a functioñ∈ 1 (Ω) so that is given as, on Ω, Equation (6) could be thought as a compatibility condition between and . To simplify the notation, from now on, we will drop the tildes.
or equivalently.
We have the following theorem.

Numerical Method
We define International Journal of Partial Differential Equations and choose test function and redefine the gradient and divergence operator ] .
Then (1) can be written as The jump condition equation (2) can be reformed as and boundary condition is 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 levelset 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 Ω + .
Two sets of grid functions are needed, and they are denoted by We cut every cube cell region ] into six tetrahedron regions. Collecting all those tetrahedron regions, we obtain a uniform tetrahedralization ℎ : ⋃ ∈ ℎ ; see Figure 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 Two extension operators are needed. The first one is ℎ : 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,ℎ 0 , is a finite dimensional subspace of 1 0 (Ω). The second extension operator ℎ is constructed as follows. For any ℎ ∈ 1,ℎ with ℎ = ℎ at boundary points, ℎ ( ℎ ) is a piecewise linear function and matches ℎ on grid points. It is a linear function in each regular cell, just like the first extension operator ℎ ( ℎ ) = ℎ ( ℎ ) in a regular cell. In each interface cell, it consists of two pieces of linear functions, one is on + and the other is on − . 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,ℎ 0 , we have 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. 6 International Journal of Partial Differential Equations

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 Here we can choose specific jump conditions 1 = 2 = 1 = 2 = 0 since changing , will not affect the matrix . In this case we have since 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. Proof. Suppose for a contradiction that is not positive definite. Then there is a vector ∈ 2 and ̸ = 0 such that ≤ 0. Let Then Since for all = [ 1 , 2 , 3 , 4 , 5 , 6 ] ∈ 6 ,

Numerical Experiments
Consider the problem International Journal of Partial Differential Equations 7 In all numerical experiments later, the level-set function ( , , ), the coefficients ± ( , , ), ± ( , , ), and ± ( , , ), and the solutions, are given. Hence 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.  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 (ℎ −2 ), the same as the case without interface.      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. (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 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 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.   solution ± 1 ( , , ), ± 2 ( , , ) in this example are given as follows: ( , , ) = 2 + 2 + 2 − 0.25, 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.

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 (ℎ −2 ), the same as the case without interface. The condition number also grows only linearly with the ratio between − and + .