The Least Squares Finite Element Method for Elasticity Interface Problem on Unfitted Mesh

In this paper, we propose and analyze the least squares finite element methods for the linear elasticity interface problem in the stress-displacement system on unfitted meshes. We consider the cases that the interface is $C^2$ or polygonal, and the exact solution $(\sigma,u)$ belongs to $H^s(div; \Omega_0 \cup \Omega_1) \times $H^{1+s}(\Omega_0 \cup \Omega_1)$ with $s>1/2$. Two types of least squares functionals are defined to seek the numerical solution. The first is defined by simply applying the $L^2$ norm least squares principle, and requires the condition $s \geq 1$. The second is defined with a discrete minus norm, which is related to the inner product in $H^{-1/2}(\Gamma)$. The use of this discrete minus norm results in a method of optimal convergence rates and allows the exact solution has the regularity of any $s>1/2$. The stability near the interface for both methods is guaranteed by the ghost penalty bilinear forms and we can derive the robust condition number estimates. The convergence rates under $L^2$ norm and the energy norm are derived for both methods. We illustrate the accuracy and the robustness of the proposed methods by a series of numerical experiments for test problems in two and three dimensions.


Introduction
In this paper, we develop the least squares finite element methods (LSFEMs) for linear elasticity interface problems, which model the elasticity structure with different or even singular material properties, and have many applications in fields of materials science and continuum mechanics [2,4,25,26,39].For such problems, the governing equations usually have discontinuous coefficients and involve the inhomogeneous jump conditions.Because of the discontinuity near the interface and the irregular geometry of the interface, it is still challenging to design efficient numerical methods for such equations.
The finite element method is an important numerical method for solving interface problems.In the last decades, various numerical schemes have been developed for the elliptic interface problem, and we refer to [6,10,11,16,33,36,40,53] for some typical methods.The finite element methods can be roughly classified into fitted and unfitted methods based on types of grids.The body-fitted method requires the mesh to be aligned with the interface for representing the geometry of the interface accurately.For complex geometries, it is a challenging and time-consuming task to generate a high quality body-fitted mesh especially in high dimensions [36].In the unfitted method, the interface description is decoupled from the generation of the mesh, which provides a good flexibility when handling the problem with complex geometries.Examples of such methods are the cut finite element method [6,10,11,30,33,58], the immersed finite element method [28,40,41] and the aggregated finite element method [3,17,37].
Recently, the unfitted finite element methods are also been applied to solve the linear elasticity interface problem.In [34], Hansbo and Hansbo proposed a linear finite element method and derived the optimal convergence rates in error measurements under the assumption that the exact solution  is piecewise  2 .In [4], Becker et al. developed a mixed finite element method with the linear accuracy for the displacement-pressure formulation.The inf-sup condition and the optimal error estimates are verified.Both methods are also called Nitsche extended finite element methods (Nitsche-XFEM), where the jump conditions are weakly imposed by the Nitsche penalty method [33].In [58], Zhang followed the interface-penalty idea and presented a high-order unfitted method for the elasticity interface problem.Combining with the penalty method and the hybridizable discontinuous Galerkin approximation, Han et al. developed an X-HDG method for this problem [32].Another type of unfitted finite element methods is the immersed finite element method [40], which mainly modifies the basis functions near the interface to capture the jump of the solution.In [45,46], the authors presented the nonconforming immersed finite element methods for the linear elasticity interface problem.Other types of immersed finite element methods can be found in [28,38].We note that all above mentioned methods derive the error estimates under the assumption that the exact solution  has at least piecewise  2 regularity.Many analysis techniques that are developed for piecewise  2 solutions in the elliptic interface problem can be used in this case.But these techniques may be unavailable for the solution that has only piecewise  1+ ( < 1) regularity.In addition, for piecewise  1 (div) or piecewise  1 (curl) functions, applying these techniques to estimate the errors will result in a suboptimal convergence rate, see [43,48] for unfitted methods in solving (curl)-and (div)-interface problems.To our best knowledge, there are few works on the interface problem of low regularity.In [29], the author presented an immersed finite element method for (curl)-interface problem with the optimal convergence rates.Only piecewise  1 (curl) regularity of the exact solution is required in the analysis.
In this paper, we develop least squares finite element methods on unfitted meshes for the linear elasticity interface problem, based on the stress-displacement formulation.For traditional linear elasticity problems, the LSFEMs have been investigated in [8, 12-14, 42, 52].The LSFEM can offer the advantage of circumventing the inf-sup condition arising in mixed methods and ensure the resulting linear system is always symmetric positive definite.As the standard LSFEM, we define a least squares functional and seek the numerical solution by minimizing the functional over finite element spaces.Two types of least squares functionals are used in this paper.The first is defined by simply applying the  2 norm least squares principle to the stress-displacement.The defined functional only involves the  2 norms and the jump conditions are also enforced in the sense of  2 norms.This method requires the exact solution (, ) has the regularity   (div; Ω 0 ∪ Ω 1 ) ×  1+ (Ω 0 ∪ Ω 1 ) with  ≥ 1.The convergence rate in the (div) ×  1 norm is shown to be half order lower than the optimal rate.From the embedding theory, we know that any  ∈ (div; Ω 0 ∪ Ω 1 ) has the normal trace  •  ∈  −1/2 (Γ), but the stronger  2 norm is applied to handle the normal trace in this method.This is the reason that the condition  ≥ 1 is required and the convergence rate is not optimal.To overcome this difficulty, we define another least squares functional with a discrete minus inner product.This method follows the ideas in [7,8], where the discrete minus norms corresponding to  −1 (Ω) are used.In this paper, we define a discrete minus inner product that is related to the inner product in the space  −1/2 (Γ).The use of this inner product allows us to relax the regularity condition as  > 1/2, and gives the optimal convergence rate under the error measurement with respect to the required regularity.It is noticeable that the optimal convergence rate under the (div) norm is achieved for functions in   (div; Ω 0 ∪ Ω 1 ).We also point out that in the unfitted methods, the  1 trace estimate is usually the main tool to estimate the numerical error on the interface.For the low regularity case  < 1, this estimate is unavailable, and we use the embedding theory in the error estimation instead.
Another important issue for unfitted methods is the presence of small cuts near the interface.In our method, we employ the ghost penalty method [9] to cure the effects bringing by small cuts.The ghost penalty bilinear forms also correspond to  2 norms, and they can be naturally added in the least squares functional.We can prove a uniform upper bound of the condition number to the resulting linear system with proper penalty forms.We give a suitable penalty bilinear form with the polynomial local extension, and also the standard penalty forms given in [9,30] can be used in our methods.
The rest of this article is organized as follows.In Section 2, we give the basic notation and introduce the stress-displacement formulation for the linear elasticity interface problem.The ghost penalty forms are also presented in this section.In Section 3, we define the associated least squares functional of the interface problem.Section 4 develops the numerical schemes.The least squares finite element methods with  2 norms and with the discrete minus norm are established in Sections 4.1 and 4.2, respectively.The error estimations are also included.Finally, numerical results of test problems in two and three dimensions are presented in Section 5.

Problem setting and preliminaries
Let Ω ⊂ R  ( = 2, 3) be a convex polygonal (polyhedral) domain with the boundary Ω.Let Ω 0 Ω be a polygonal (polyhedral) subdomain or a subdomain with a  2 -smooth boundary.We denote by Γ := Ω 0 the topological boundary, which can be regarded as an interface dividing Ω into two disjoint domains Ω 0 and Ω 1 , where Ω 1 := Ω∖Ω 0 , Ω 0 ∩ Ω 1 = ∅ and Ω 0 ∪ Ω 1 = Ω.The model problem considered in this paper is the linear elasticity interface problem defined on Ω, which reads: seek the stress  = (  ) × and the displacement where  is the source term, and ,  are the jump conditions on the interface.The jump operators are defined as (6).The Lamé parameters ,  are assumed to be piecewise positive constant functions, ((), ()) := The constitutive law is expressed by the linear operator  : R × → R × : where tr(•) denotes the trace operator, and I := (  ) × is the identity tensor.The function () denotes the symmetric strain tensor: We assume that the interface problem (1) admits a unique solution (, ) ∈ Σ  × V +1 with  > 1/2, where For  = 0, 1, we let   := | Ω ,   := | Ω .We further assume that   and   can be extended to the whole domain Ω in the sense that there exist ̃︀   ∈ H  (div; Ω) and ̃︀ . Consequently,  and  can be decomposed as where   is the characteristic function corresponding to the domain Ω  .We refer to [1,20,35] for more details about the extension of Sobolev spaces.From (4), we formally introduce two projection operators   ( = 0, 1) that    := ̃︀   ∈ H  (div; Ω) and    := ̃︀   ∈ H +1 (Ω).Let us introduce the notation required in the definition to the scheme.We denote by  ℎ a quasi-uniform partition of Ω into triangle (tetrahedron) elements.The mesh  ℎ is unfitted that the element faces in  ℎ are not required to be aligned with the interface Γ.For any element  ∈  ℎ , we denote by ℎ  its diameter and by   the radius of the largest disk (ball) inscribed in .Let ℎ := max ∈ ℎ ℎ  be the mesh size, and let  := min ∈ ℎ   .The mesh  ℎ is quasi-uniform in the sense that there exists a constant   independent of ℎ such that ℎ ≤   .
Throughout this paper,  and  with subscripts are denoted to be generic positive constants that may vary in different lines, but are always independent of the mesh size ℎ, how the interface Γ cuts the mesh  ℎ , and the Lamé parameter  defined in (2).
To cure the effect bringing by small cuts near the interface, we follow the idea in the ghost penalty method [9,30], which uses the data from the interior domain to ensure the stability near the interface.For this goal, we assume that we can construct two bilinear forms   ℎ,0 : ℎ, (, )( = 0, 1) for ∀ ∈  2 (Ω ℎ, ).In our method, we assume that the forms satisfy the following two properties: P1: the  2 norm extension property: P2: the weak consistency: The suitable penalty forms   ℎ, (•, •) can be constructed by the face-based penalties and the projection-based penalties, see [10,30].In Section 2.7 of [30], the penalty forms constructed for elliptic interface problems satisfy P1 and P2, which can be used in our method.We also note that extra properties of the penalty form in solving elliptic interface problems are required, as the  1 seminorm extension property and the inverse estimate, see EP1-EP4 of [30].In our method, P1 and P2 are enough to ensure the stability near the interface.Thus, the required ghost penalty may be more simple.
Here, we outline a method to construct the penalty bilinear forms by the local polynomial extension, and the implementation is easy.The idea of the local extension has also been widely used in unfitted methods, see [3,11,17,36,37,57].
We close this section by giving the  1 trace estimate [31,33,53] on the interface.
Lemma 1.There exists a constant  such that We refer to [31] for the proof only assuming Γ is Lipschitz.The  1 trace estimate is fundamental in the penalty-type unfitted finite element methods, as the main tool to handle the numerical error on the interface, such as [3,10,30,33,36,37,49,53].By (13), the error on the interface can be bounded by the estimates on elements.However, this trace estimate (13) requires the  1 regularity.For the problem with the exact solution in   (div; Ω 0 ∪Ω 1 ) or   (curl; Ω 0 ∪Ω 1 ), applying (13) to estimate the numerical errors on the interface will lead to a suboptimal convergence rate.The optimal convergence needs a higher regularity assumption that the exact solution is piecewise  +1 -smooth.see [43,48] for unfitted methods on (div)-and (curl)-interface problems.In addition, the trace estimate ( 13) is not suitable for the solution of low regularity, such as the solution belongs to the space   (Ω 0 ∪ Ω 1 ) with  < 1.

Least squares functional for the interface problem
In this section, we introduce an associated least squares functional to the interface problem (1).Let Σ := Σ 0 and V := V 1 be the spaces coinciding with  = 0 in (3).We define the quadratic functional  (•; •) by where ( , ; The trace terms   (•; •) and   (•; the trace theory.The exact solution (, ) clearly minimizes the functional  (•; •) since  (, ;  , , ) = 0.In fact, we can prove that (, ) is also the unique solution to the minimization problem inf ( ,)∈Σ×V  ( , ;  , , ).For this purpose, we will give a norm equivalence property of  (•, •), and the norm equivalence is also crucial for the error analysis in the least squares finite element method [5].The pair of spaces Σ × V can be naturally equipped with the norm The equivalence between  (•; •) and ‖ • ‖ e is given in the following lemma.
Lemma 2. There exist constants  such that Proof.By the definition of , we have that F. YANG The term ‖ − ()‖  2 (Ω0∪Ω1) can be bounded by ‖( , )‖ e , following from the triangle inequality and (17).
From the embedding theory, we know that Similarly, there holds . The upper bound in ( 16) is reached.From Theorem 3.1 of [13], the lower bound in ( 16) holds for  ∈ H(div; Ω) and  ∈ H 1 (Ω), i.e. the lower bound holds for  and  satisfying For the general case, we construct two auxiliary functions ︀  and ̃︀  to prove (16).Consider the elliptic problems and . We then extend  0 and  1 to the domain Ω by zero.Let ̃︀  :=  + ∇ 0 and ̃︀ Combining the above two estimates leads to the lower bound in (16), which completes the proof.
It is noted that the generic constants in ( 16) are independent of , which ensure the proposed methods are robust when  → ∞.

Least squares finite element methods
In this section, we present the numerical schemes for solving the elasticity interface problem (1).We begin by introducing the approximation finite element spaces.For the mesh  ℎ, ( = 0, 1), we define Σ  ℎ, ⊂ H(div; Ω ℎ, ) as the (div)-conforming tensor-valued piecewise polynomial space of degree .As the definition of H(div; Ω ℎ, ), each column of functions in Σ  ℎ, belongs to (div; Ω ℎ, ), i.e. the BDM  space or the RT  space.In this paper, the schemes are established on Σ  ℎ, is the BDM  space.The methods and the analysis can be extended to the RT  space without any difficulty.We define V  ℎ,0 ⊂ H 1 (Ω ℎ,0 ) as the vector-valued  0 finite element space of degree , and define V  ℎ,1 ⊂ H 1 (Ω ℎ,1 ) as the  0 finite element space with zero trace on Ω, i.e.V  ℎ,1 | Ω = 0.
We define the extended approximation spaces •  1 for the stress and the displacement, respectively, where the characteristic function   is defined in (4).Then, any  ℎ ∈ Σ  ℎ and any  ℎ ∈ V  ℎ admit a unique decomposition that where  ℎ, ∈ Σ  ℎ, ,  ℎ, ∈ V  ℎ, .As the decomposition (4), we also formally let   ( = 0, 1) be the projection operator such that . The approximation spaces are conforming in the sense that Σ  ℎ ⊂ Σ and V  ℎ ⊂ V.A natural idea to seek numerical solutions is minimizing the functional  (•; •) over the discrete conforming spaces, i.e. inf ( ℎ , ℎ )∈Σ  ℎ ×V  ℎ  ( ℎ ,  ℎ ;  , , ), which, however, is not an easy task.It is impossible to directly compute the trace terms   (•; •) and   (•; •) in (15), and the minimization problem cannot be readily rewritten into a variational problem by the means of Euler-Lagrange equation because of the presence of . The main idea of the proposed method is to apply some computational trace terms to replace   (•; •) and   (•; •) in  (•; •) to define new functionals.The numerical approximations are then obtained by minimizing the new functionals.
For the convergence analysis, we define the spaces Σ ℎ := {} ∪ Σ  ℎ and V ℎ := {} ∪ V  ℎ , which contain the exact solution and all finite element functions, respectively.Here, we introduce the ghost penalty forms for the spaces Σ ℎ and V ℎ from the forms (12).By the definition of   ( = 0, 1), there holds 12), we define the forms From properties P1 and P2, we have the following estimates, and 4.1. 2 norm least squares finite element method for  ≥ 1 We present the numerical method that defines the least squares functional using only  2 norms.This method requires the exact solution has the regularity  ≥ 1.As mentioned earlier, we will bound the trace terms   (•; •) and   (•; •) by applying stronger and easily-computational norms.In this case, the space For the displacement variable, we also give a stronger norm for ‖•‖  1/2 (Γ) .We first consider the case that Γ is  2smooth.By the embedding relationship , where ∇ Γ denotes the tangential gradient on the interface [23].The tangential gradient ∇ Γ () for ∀ ∈ H 1 (Γ) only depends on values of  on Γ ∩  , where  is a neighbourhood of .For  = 0, 1, we note that any finite element function  ℎ, ∈ V  ℎ, ( = 0, 1) is continuous and piecewise smooth on Γ, i.e.  ℎ, | Γ ∈  0 (Γ) and By the Sobolev interpolation inequality ( [47], Thm.7.4) and the Cauchy-Schwarz inequality, we have that The right hand side of ( 22) only involves the  2 norms, and can be easily computed.
For the case that Γ is polygonal (polyhedral), we let Γ  (1 ≤  ≤ ) denote the sides of Γ, where every side Γ  is a line segment (polygon).Since Γ has corners, the estimate ( 22) will be modified in a piecewise manner.From the embedding theory ( [27], Thm.1.5.2.1), we have that ) is continuous and piecewise smooth on each Γ  .We can know that ( The second inequality follows from the Sobolev interpolation inequality on each Γ  .Compared with (22), the norm for the tangential gradient on Γ in (23) can be understood in a piecewise manner, i.e. we can formally let for the case Γ is polygonal.Then, the upper bound of the estimate (23) will seem the same as the upper bound of (22).To simplify the representation, we use the notation consistent with the case that Γ is  2 .Now, we define the discrete least squares functional  ℎ (•; •) that only involves the  2 norms by where (•; •) is defined in (15), and The last two seminorms in (24) are applied to guarantee the uniform upper bound of the condition number of the resulting linear system.The numerical solutions are sought by solving the minimization problem inf Since all terms in  ℎ (•; •) are defined with  2 norms, the problem (26) can be solved by writing the corresponding Euler-Lagrange equation, which reads: seek where the forms   ℎ (•; •) and   ℎ (•) are defined as and For the convergence analysis, we are aiming to show that   ℎ (•; •) is bounded and coercive, and satisfies a weakened Galerkin orthogonality property.Then, the error estimate can be reached in a standard argument.
We introduce the energy norm ‖ • ‖ e ℎ by By ( 21)-( 23), we have that and Combining the definitions to  (•; •) and  ℎ (•; •) and the above two estimates, there holds From Lemma 2, we find that The boundedness and the coercivity of the bilinear form   ℎ (•; •) directly follows from the Cauchy-Schwarz inequality, the estimate (17) and the estimate (30).
Lemma 3.There exist constants  such that Bringing the exact solution (, ) into the discrete variational problem (27) yields that which leads to a weakened Galerkin orthogonality property of   ℎ (•; •).
The error estimation to the numerical solution of (24) follows from Lemmas 3 to 5.
We estimate the condition number for the linear system (27), which is especially desired in the unfitted method.
Theorem 2. There exists a constant  such that where   ℎ is the linear system with respect to   ℎ (•; •).
Proof.Since the bilinear form   ℎ (•; •) is bounded and coercive with respect to the norm ‖•‖ e ℎ .From Section 3.2 of [24], the main step is to show the relationship between the energy norm ‖ • ‖ e ℎ and the  2 norm.Note that the finite element spaces Σ  ℎ and V  ℎ are the combinations of Σ  ℎ,0 and Σ  ℎ,1 and of V  ℎ,0 and V  ℎ,1 , respectively, where the spaces Σ  ℎ, and V  ℎ, are defined on the domain Ω ℎ, .Hence, our goal is to show that Since ‖ • ‖ e ℎ is stronger than ‖ • ‖ e , we have that The lower bound in (42) then follows from the estimates (19) and (20).For the upper bound, we apply the triangle inequality and the inverse estimate to find that )︁ .

F. YANG
We have completed the error analysis for the discrete variational form (27).The scheme is robust in the sense that the constants appearing in the error bounds (39) and ( 40) are independent of , and also are independent of how the interface Γ cuts the mesh.From Theorem 1, the convergence rate of the numerical error under the energy norm is half order lower than the optimal rate without the higher regularity assumption.The major reason is that we use the stronger  2 norm to replace the minus norm ‖ • ‖  −1/2 (Γ) to define the new quadratic functional, and the errors on the interface are essentially established by the  1 trace estimate as (37) and (38).In addition, both estimates require the exact solution (, ) has at least  1 (div) ×  2 regularity.Consequently, this method and the analysis cannot be extended to the case of low regularity that  < 1.

Least squares finite element method with the discrete minus norm for 𝑠 > 1/2
In this subsection, we give the numerical method that has the optimal convergence speed for  > 1/2.As stated before, the replacement of the  2 norm cannot work and the  1 trace estimate is also inappropriate for the case of low regularity.The natural choice is to involve the ‖ • ‖  −1/2 (Γ) in the discrete least squares functional, but the minus norm is not easy to compute.Here we follow the idea in [7,8] to employ a discrete inner product, which is related to the minus norm ‖ • ‖  −1/2 (Γ) .We first give an inner product in H −1/2 (Γ).
Roughly speaking, we have defined an equivalent norm and an inner product for the space H −1/2 (Γ), although the operator  and the inner product (•, •) −1/2,Γ still cannot be directly computed.In the numerical scheme, we introduce a discrete operator  ℎ to replace  to make the method computationally feasible.The discrete operator  ℎ is established with the space V 1 ℎ,0 , which is the  0 linear finite element space defined on the partition  ℎ,0 .Given any  ∈ H −1/2 (Γ), we define the following discrete weak form: seek Remark 2. The discrete variational form (45) can be regarded as an unfitted scheme for the elliptic system (43) based on the mesh  ℎ,0 .The ghost penalty form  1 ℎ,0 (•, •) ensures the uniform upper bound of the condition number to the corresponding linear system.From Theorem 2.16 of [30], the upper bound is (ℎ −2 ) for (45).
It can be easily seen that the problem (45) admits a unique solution in V 1 ℎ,0 .Similarly, this property allows us to define an operator  ℎ : ℎ,0 that  ℎ  is the solution to the discrete problem (45) for ∀ ∈ H −1/2 (Γ).As (44), we define the discrete inner product (•, •) −1/2,ℎ and the induced norm Let  ℎ =  ℎ in ( 45), together with the property ( 8), it can be seen that and Then, we give the relationship between  and  ℎ .
From the Sobolev extension theory [1], there exists a linear extension operator  0 : ℎ,0 be the Scott-Zhang interpolation operator of the space V 1 ℎ,0 , which satisfies the approximation estimates From the trace estimate (13), it is quite standard to derive that ‖ −  1 ℎ,0 ‖  2 (Γ) ≤ ℎ 1/2 ‖‖  1 (Ω ℎ,0 ) for ∀ ∈  1 (Ω ℎ,0 ).Then, we deduce that The second term can be bounded by the approximation property, that is We let  ℎ =  1 ℎ,0 ( 0 ( )) in (45) and apply (48) to bound the first term, Collecting all above estimates yields the second estimate (50), which completes the proof.Now, let us define the least squares functional ̃︀  ℎ (•; •) by where   ℎ (•; •) is defined as in (25) and ̃︀   ℎ (•; •) is defined as The numerical solutions are sought by minimizing the functional ̃︀  ℎ (•; •) over the finite element spaces Σ  ℎ ×V  ℎ .This minimization problem is equivalent to a variational problem by writing the Euler-Lagrange equation, which reads: seek where the forms  ̃︀  ℎ (•; •) and  ̃︀  ℎ (•) are defined as and The convergence analysis of the solution in ( 53) is analogous to the discrete variational form (27).We introduce the energy norm ‖ • ‖ ̃︀ e ℎ by From Lemma 6, there holds Together with the estimate (29), we can know that As Lemma 3, it is similar to get that  ̃︀  ℎ (•; •) is bounded and coercive under the energy norm ‖•‖ ̃︀ e ℎ , and satisfies a weakened Galerkin orthogonality property.Lemma 7.There exist constants  such that Lemma 8. Let (, ) ∈ Σ  × V +1 be the exact solution to (1), let ( ℎ ,  ℎ ) ∈ Σ  ℎ × V  ℎ be the numerical solution to (27), there holds Proof.The proofs are similar to Lemmas 3 and 4.
Moreover, we give the approximation estimate under the norm ‖ • ‖ ̃︀ e ℎ .
, which is the same as in the proof to Lemma 5.
The rest is to bound the terms ℎ‖[ . For the case  ≥ 1, both errors can be bounded by the  1 trace estimate, and we derive that and similarly, there holds For the case 1/2 <  < 1, the error ℎ‖[ [ −   ] ]  ‖ 2  2 (Γ) cannot be bounded as (58) because the  1 trace estimate is now unavailable.In this case, the embedding theory will be the main tool to estimate the errors.We let V  ℎ ⊂ H 1 (Ω) be the tensor-valued  0 finite element space of degree  on the mesh  ℎ .Since    ∈ H  (Ω), we let  s   be its Scott-Zhang interpolant into the space V  ℎ [18], which satisfies that ‖ s   −   ‖  2 (Ω) ≤ ℎ 0 ‖  ‖   (Ω) .For arbitrarily small  > 0, the space H 1/2+ (Ω  ) continuously embeds into L 2 (Γ).Notice that  s  ,    ∈ H 1/2+ (Ω  ) for small enough , and we apply the inverse estimate and the triangle inequality to derive that where  1 = min(, ) − .From the embedding H 3/2+ (Ω  ) ˓→ H 1 (Γ) [21], it is similar to get that Collecting all above estimates leads to the approximation property (57), which completes the proof.
The error estimation for the numerical solution of ( 53) can be reached.The proof is the same as Theorem 1.
The condition number of the linear system of  ̃︀  ℎ (•; •) also has a uniform upper bound.The proof is the same to Theorem 2. Theorem 4.There exists  such that where  ̃︀  ℎ is the linear system with respect to  ̃︀  ℎ (•; •).From Theorem 3, it can be seen that the scheme is also robust since the constants in the error bounds (59) are independent of  and the relative location between Γ and the mesh.Compared with the  2 norm least squares finite element method, the convergence rate is nearly optimal without a higher regularity assumption, but the linear system is more complicated.
Ultimately, let us give some details about solving the linear system  ̃︀  ℎ .Since  ̃︀  ℎ (•; •) involves the discrete minus inner product (•, •) −1/2,ℎ , we are required to solve the elliptic system (45) in the assembling of  ̃︀  ℎ .Let  be the sparse matrix of the bilinear form (45), which is invertible and satisfies that () ≤ ℎ −2 , see Remark 2. Let  be the sparse matrix corresponding to the inner product ([ ℎ,0 .We note that this definition ensures the matrix  has the same column size as the matrix  ̃︀  ℎ .Then,  ̃︀  ℎ has the form that where   2 corresponds to all  2 inner products in  ̃︀  ℎ (•, •), i.e.   2 is the matrix of the bilinear form   2 (•; •) that Suppose that we have a fast algorithm to solve the linear system y = z, then we can easily compute the matrix-vector products for the matrix  ̃︀  ℎ by (61).Consequently, the Krylov iterative method (e.g.GMRES) can be used as the solver for the linear system.Although we have shown that the condition number of  ̃︀  ℎ is (ℎ −2 ), an effective preconditioner is still expected especially when ℎ tends to zero.But  ̃︀  ℎ involves the inverse matrix  −1 , it is costly to form it and it is also not convenient to even extract its diagonal.Traditional diagonal or block diagonal preconditioners are hard to use for  ̃︀  ℎ .One method is to use the matrix-free preconditioning technique to  ̃︀  ℎ .For example, one can use the hierarchically semiseparable (HSS) approximation for the given matrix to accelerate the convergence of iterative methods [15,54,55].The construction of the HSS approximation may potentially only use matrix-vector products instead of the original matrix itself, and we refer to [44,54] for fully matrix-free techniques to the construction of the HSS approximation.This idea has been used in the immersed finite element method solving the elliptic interface problem [56] .Once the HSS approximation  for  ̃︀  ℎ is obtained in a structured form, it can be quickly factorized and the factors can be used as a preconditioner.Alternatively, we present a matrix-explicit preconditioning method for  ̃︀  ℎ .As (61), the matrix  ̃︀  ℎ can be split into two parts.From (62), it can be easily seen that   2 ( ℎ ,  ℎ ;  ℎ ,  ℎ ) = 0 implies  ℎ = 0,  ℎ = 0 for ( ℎ ,  ℎ ) ∈ Σ  ℎ × V  ℎ .Hence,   2 is symmetric positive definite.Any preconditioning technique can be applied to   2 because   2 is entirely explicit.As a numerical observation, the preconditioner from   2 can also significantly accelerate the iterative methods.We show the results in Example 1 of Section 5 by constructing the HSS approximation from   2 to obtain the preconditioner.For solving the linear system y = z, we also use the HSS approximation as the preconditioner.The codes of the construction and the factorization of the HSS approximation are freely available in STRUMPACK [50].A comprehensive analysis and a more appropriate method to solve the linear system are now left in the future study.

Numerical results
In this section, a series of numerical tests are presented to show the numerical performance of the proposed method.For all tests, the source function  and the jump condition ,  are taken from the exact solution accordingly.In Examples 1-5, we consider the elasticity interface problems in two dimensions on the squared domain Ω = (−1, 1) 2 .We adopt a family of triangular meshes with the mesh size ℎ = 1/5, 1/10, . . ., 1/40 to solve all tests.In Examples 1-3, the interface Γ is of class  2 and is described by a level set function, see Figure 1.In Examples 4 and 5, the interface is taken to be the boundary of the L-shaped domain.In Example 6, we solve a three-dimensional interface problem to illustrate the numerical performance.Example 1.We first consider a linear elasticity interface problem with a circular interface centered at the origin with radius  = 0.7, see Figure 1.We solve this problem to study the convergence rates.The exact solution is taken as with the discontinuous parameters | Ω0 =  0 = 5, | Ω1 =  1 = 1, | Ω0 =  0 = 2, | Ω1 =  1 = 1.The numerical results are gathered in Table 1 by the  2 norm least squares finite element method and the least squares finite element method with the discrete minus norm.For both methods, the energy norms ‖ • ‖ e ℎ and ‖ • ‖ ̃︀ e ℎ are stronger than the norm ‖ • ‖ e .We calculate ‖( −  ℎ ,  −  ℎ )‖ e to represent the errors under the energy norm.
From Table 1, it can be seen that the errors under the energy norm approach zero at the optimal rates.For the  2 errors, ‖ −  ℎ ‖  2 (Ω) and ‖ −  ℎ ‖  2 (Ω) have optimal/suboptimal convergence speeds.All numerically observed convergence orders are in agreement with the theoretical analysis in Theorems 1 and 3.
In this example, we demonstrate the numerical performance of the iterative method for solving the resulting linear system.For all tests, we use GMRES as the iterative solver and the iteration stops when the relative error ‖x  − b‖  2 /‖b‖  2 < 10 −8 at the stage  is smaller than the tolerance 10 −8 .For the  2 norm least squares finite element method, the resulting matrix is entirely explicit and any preconditioning technique can be used.Here, we try to use the HSS approximation and its diagonal to construct the preconditioners.The convergence steps are collected in Table 2.It can be seen that the convergence is apparently accelerated by the HSS preconditioning technique.
For the method with the discrete minus norm, we shall solve the linear system  ̃︀  ℎ x = b.As stated in the end of Section 4.2, we construct the HSS approximation   2 from   2 to obtain the preconditioner, and we also directly construct the HSS approximation  from  to obtain the preconditioner as a comparison.The  convergence steps are collected in Table 3.The convergence speeds of both preconditioned GMRES methods are significantly faster than the standard GMRES method.The convergence histories for the iterative methods on the mesh ℎ = 1/20 are depicted in Figure 2. The numerical performances of both preconditioned methods are very close, and it is more convenient to construct the HSS approximation from   2 in a standard procedure.In every step of Krylov iteration, we are required to solve the linear system y = z.We still use GMRES with the preconditioner from the HSS approximation as the solver for this system.The average iterative steps for different meshes are collected in Table 3.The iterative method also has a fast convergence speed.In the rest examples, we use the factors of   2 as the preconditioner to solve the linear system.Developing a more appropriate method to solve the linear system is now left in the future study.
Example 2. In this example, we consider the same interface problem as Example 1 but the Lamé parameters have a large jump across the interface.We solve this problem to demonstrate the robustness of the proposed method when  → ∞.The parameters are selected as  0 = 10000/100,  1 = 1,  0 =  1 = 1.The numerical results for both methods are displayed in Tables 4 and 5, respectively.We observe that the discrete solutions from both methods converge uniformly as  → ∞.The results illustrates the robustness when the parameter  → ∞.Example 3. In this example, we consider an elasticity interface problem with a star-shaped interface consisting of both concave and convex curve segments, see Figure 1.The interface Γ is governed by the polar angle  that The analytic solution  is taken in a piecewise manner as where  0 (, ) = [cos() cos(), cos()]  ,  1 (, ) = [sin() sin(), (1 − ) sin()]  .
The parameters are chosen as  0 =  1 =  0 =  1 = 1.We present the numerical results in Table 6.The predicted convergence rates under the energy norms and the  2 norms for both unknowns are verified from the convergence histories.The exact solution and the parameters are taken the same as Example 1.The numerical errors are displayed in Table 7.For the case of the polygonal interface, all numerically detected convergence speeds still agree with the theoretical results.The source term  = 0 on the domain Ω 0 ∪Ω 1 .The exact solution has the regularity that (, ) ∈ H − (div; Ω 0 ∪ Ω 1 ) × H 1+− (Ω 0 ∪ Ω 1 ) for ∀ > 0. Hence, we solve this problem by the least squares finite element with the discrete minus norm.We consider the linear accuracy  = 1 for this test, and the results are shown in Table 8.It can be observed that the convergence rate for the energy norm is (ℎ 0.52 ), which is consistent with the regularity of the exact solution and the theoretical analysis.The  2 errors for both variables are less than the optimal convergence rates, i.e. (ℎ 1.3 ) and (ℎ 0.51 ) for  and , respectively.The reason may be traced back to the singularity of the exact solution.
Example 6.In this test, we consider the interface problem in the cubic domain Ω = (0, 1) 3 .We take Γ as a spherical interface with the radius  = 0.35 centered at the point (0.5, 0.5, 0. The parameters ,  are discontinuous across the interface with  0 = 10,  1 = 1,  0 = 10,  1 = 1.We adopt a series of tetrahedral meshes with the mesh ℎ = 1/4, 1/8, 1/16, 1/32 to solve this problem, see Figure 3.The numerical errors under all error measurements are reported in Table 9 for both methods.The numerical results illustrate the accuracy of the methods in three dimensions.

Figure 1 .
Figure 1.The unfitted mesh and the interfaces in two dimensions.

Figure 3 .
Figure 3.The spherical domain and the tetrahedral mesh of Example 3.

Table 1 .
The numerical results for Example 1 by the  2 norm least squares finite element method (left)/the least squares finite element method with the discrete minus norm (right).

Table 2 .
Convergence steps for the linear system of the  2 norm least squares finite element method in Example 1.

Table 3 .
Convergence steps for the linear system of the least squares finite element method with the discrete minus norm in Example 1.

Table 4 .
Numerical results for Example 2 by the  2 norm least squares finite element method with  = 100, 10000.

Table 5 .
Numerical results for Example 2 by the least squares finite element method of the discrete minus norm with  = 100, 10000.

Table 6 .
Numerical results for Example 3 by the  2 norm least squares finite element method (left)/the least squares finite element method with the discrete minus norm (right).

Table 9 .
Numerical results for Example 6 by the  2 norm least squares finite element method (left)/the least squares finite element method with the discrete minus norm (right).