1 Introduction
Interface problems arise from diverse physical and engineering applications in which the coefficients of the governing partial differential equations are discontinuous across material interfaces that separate the physical domains. The bodyfitted finite element methods resolve the geometry of the interface by requiring the vertices of the finite element mesh located on the interfaces
[3, 23, 19]. For domains with complex geometry, the construction of bodyfitted shape regular finite element meshes may be difficult and timeconsuming, which is the main driving force of the study of unfitted finite element methods. In this paper we will show that the shape regular bodyfitted mesh can indeed be constructed for any shaped smooth interface based on our new merging cell algorithm (see remarks below Theorem 3.1). We emphasize, however, that even when the bodyfitted shape regular mesh is available, the construction of highorder finite element methods still requires substantial new ideas including, for example, the isoparametric finite element method [24, 34] or unfitted finite element methods which are the focus of this paper. We remark that the shape regularity assumption of the finite element mesh is not only fundamental in the mathematical theory of finite element methods (see, e.g., [24]) but also essential in controlling the condition number of the finite element stiffness matrix for elliptic equations (see, e.g. [10]).Let be a bounded Lipschitz domain which is divided by a smooth interface into two nonintersecting subdomains , , see Fig.1.1. We consider the following elliptic interface problem
(1.1)  
(1.2) 
where , , is the unit outer normal to , and stands for the jump of a function across the interface . We assume that the coefficient is positive and piecewise constant, namely, , , where
denotes the characteristic function of
, .Unfitted finite element methods in the discontinuous Galerkin (DG) framework have attracted considerable interests in the literature in the last twenty years starting from the seminal work [31] in which an unfitted finite element method is proposed for elliptic interface problems. The method is defined on a fixed background mesh and uses different finite element functions in different cut cells which is the intersection of the mesh elements and physical domains. The jump condition on the interface is enforced by penalties which extends an earlier idea of Nitsche [41]. This unfitted finite element method can also be viewed as the interior penalty discontinuous Galerkin method (see, e.g., [2]) defined on meshes allowing curveshaped elements. The main difficulty in using the unfitted finite element methods is the socalled small cut cell problem: the cut cells can be arbitrarily small and anisotropic, which can make the stiffness matrix extremely illconditioned, especially for highorder finite element methods [44, 8]. For other approaches to design unfitted discretization methods by constructing special finite element bases on interface elements or finite difference stencils along the interface, we refer to the immersed boundary method [43], the immersed interface method [35, 36], or the immersed finite element method [37, 22].
There are two approaches in the literature to attack the small cut cell problem. One is by appropriate techniques of stabilization [16, 17, 47, 38, 48, 30]. Among them, for example, the method of ghost penalty [16, 17, 30] adds additional penalties on the jumps of derivatives across sides or facets of interface elements. The other approach is by merging the small cut cells with neighboring large elements [33, 32, 9, 21, 15] so that the merged macroelements have enough support. While the DG formulation is still used in [33, 32, 21], the aggregated unfitted finite element method in [9] relies on the construction of stable extension operators so that the finite element space is still . We refer to recent works [13, 18, 7, 8] for further information about ghost penalty and the aggregated unfitted finite element method.
In [21] an adaptive highorder unfitted finite element method is proposed for elliptic interface problems in which the
a priori and a posteriori error estimates are derived based on novel
domain inverse estimates and the concept of interface deviation. The interface deviation is a measure to quantify the mesh resolution of the geometry of the interface. The macroelements, which are the union of small interface elements and their surrounding elements, are assumed to be rectangular in [21]. This assumption is different from those in [33, 32, 9], see Fig.1.2. The macroelements in [33, 32, 9] need not to be of rectangular shape, which makes the implementation simpler but the crucial inverse estimates on extended elements in [33, 32] or the stability of the extension operators [9] are shown without considering the dependence on the finite element approximation order . The assumption that the macroelements should be rectangular in [21] raises the question of how to construct the merging algorithm in practical applications.The first objective of this paper is to propose a reliable algorithm to merge small interface elements with their surrounding elements to generate the macroelements. The algorithm is based on the concept of admissible chain of interface elements, the classification of patterns for merging elements, and appropriate ordering in generating macroelements from the patterns so that the reliability of the algorithm in the sense that it terminates in finite number of steps can be proved. This algorithm also leads to a reliable algorithm of automatically generating 2D shape regular bodyfitted finite element meshes for arbitrarily shaped smooth interfaces. To the authors’ best knowledge, this algorithm introduces a new way to generate bodyfitted finite element meshes and may be of independent interest.
The second objective of the paper is to study the condition number of the stiffness matrix of highorder unfitted finite element methods which are known to be of the order in the literature [16, 33, 32, 9, 6] on quasiuniform meshes with the mesh size . For high order methods, it is known [44] that the condition number of the stiffness matrix may grow exponentially with the finite element approximation order in terms of the measure of cut cells. This indicates that the geometry of the cut cells is essential in controlling the condition number of the stiffness matrix.
In this paper, we will take the basis functions of the spectral element, that is, the Lagrangian interpolation functions at the GaussLobatto points on elements not intersecting with the interface. For the interface elements, extra care must be taken as the basis of the spectral element on
is illconditioned on the subsets , , which is similar to the observation in [27, P.346] for Legendre polynomials. Here we choose the orthogonal functions on some special polygons inside , , as the basis functions for the interface elements . We show that the condition number of the stiffness matrix is bounded by up to a logarithmic factor, where is the number of total elements, is the number of interface elements, and depends on the interface deviation and . This bound is optimal and indicates that the mesh has to sufficiently resolve the geometry of the interface to control the condition number of the stiffness matrix.The results of this paper allow for extensions in several directions. Firstly, for the ease of exposition, we consider in this paper the case when the domain is a union of rectangles and the interface is smooth. The extension to the general domains with smooth boundary is straightforward. Secondly, the case when the interface is piecewise smooth will be pursued in our forthcoming work by combining the ideas in [21] on large elements and interface deviation for interfaces with singularities with the merging algorithm developed in this paper. Thirdly, the theoretical results in this paper and in [21] including the domain inverse estimates and the concept of the interface deviation can be extended to study threedimensional interface problems. The merging algorithm in the threedimensional case is more challenging. Nevertheless, we believe that with the new insights gained in this paper for the twodimensional case, reliable algorithms for constructing cubic macroelements can be achieved in future. Finally, we remark that our argument to analyze and control the condition number of the stiffness matrix is fairly general, it can be used in other unfitted finite element methods including threedimensional cases.
The layout of the paper is as follows. In section 2 we introduce our unfitted finite element method. In section 3 we construct the merging algorithm to generate the induced mesh. In section 4 we prove the discrete Poincaré inequality and the estimate for the condition number of the stiffness matrix. In section 5 we present several numerical examples to confirm our theoretical results.
2 The unfitted finite element method
Let be a domain which is a union of rectangles and a Cartesian finite element mesh of with possible hanging nodes. This allows us to locally refine the mesh near the interface to resolve the geometry to save the computational costs away from the interface. The elements of the mesh are (open) rectangles whose sides are parallel to the coordinate axes. We assume that the interface intersects the boundary of twice at different sides (including the end points).
For any element , let stand for its diameter. Denote the set of interface elements. We recall the definition of large element in Chen et al [21, Definition 2.1].
Definiton 2.1.
(Large element) For , an element is called a large element with respect to if or for which there exists a constant such that for each side of having nonempty intersection with . Specially, is called a large element if is large with respect to both and . Otherwise, is called a small element.
Note that it is possible that may not be a large element. The following assumption in [21] is inspired by Johansson and Larson [33] in which a fictitious boundary method is considered.
Assumption (H1) For each , there exists a rectangular macroelement which is a union of and its surrounding element (or elements) such that is a large element. We assume for some constant .
In section 3 we will construct a merging algorithm to find the macroelement for each small element in an admissible chain of interface elements. This indicates that the assumption (H1) can always be satisfied by using the algorithm. In the following, we will always set if is a large element. Then, the induced mesh of is defined as
We will write . Note that is also a Cartesian mesh of in the sense that either or for any two different elements . All elements in are large elements.
For any , denote , , , and the open line segment connecting the two intersection points of and . divides the element into two polygons and which are the polygonal approximation of and , respectively. An important property of being a large element is that , , is a strongly shape regular polygon in the sense that it is the union of shape regular triangles in the sense of Ciarlet [24]. We remark that there are different definitions of shape regular polygons in the literature, see, e.g., Ming and Shi [40] and Brenner and Sung [14].
The following concept of interface deviation is introduced in [21].
Definiton 2.2.
For any , the interface deviation is defined as , where for , if is the vertex of which has the maximum distance to among all vertices of in ,
Here and .
The interface deviation is a measure on how well the mesh resolves the geometry of the interface. We will show in section 4 that this concept also links to the control of the condition number of the stiffness matrix.
It is known that if is smooth, (see, e.g., Feistauer [29, §3.3.2]) and thus for some constant independent of . Therefore, the interface deviation can be made arbitrarily small by locally refining the mesh near the interface. When the interface is Lipschitz and piecewise smooth, the definition of the large element and interface deviation has to be modified in the elements containing the singular points of the interface, see [21] for the details.
For any integer and , denote the set of polynomials in which is of degree in each variable. The following domain inverse estimate is proved in [21, Lemma 2.4].
Lemma 2.1.
Let be a triangle with vertices , , , where . Let and . Then we have
where .
The proof of this lemma makes use of the following onedimensional domain inverse estimate in [21, Lemma 2.3]
(2.1) 
where , , and is the set of polynomials of order in . We remark that the growing factor is sharp which is attained by the Chebyshev polynomials whose explicit expression is , , see DeVore and Lorentz [26, P.76].
Let , We also define two polygons , , as follows. Let be the line segment which is parallel to and its distance to is . Let be the polygon bounded by sides of and .
Lemma 2.2.
Let and , Then for , we have
(2.2) 
where the constant is independent of , , and .
Proof.
The numerical results in Example 1 in section 5 indicate that the bound in Lemma 2.2 is sharp. Now for any , we denote
Based on the concept of interface deviation, the following inverse estimates on curved domains are proved in [21, Lemma 2.8, (2.12)].
Lemma 2.3.
Let and , Then for , we have
where the constant is independent of and .
We remark that inverse estimates on starshaped curve elements are studied in Massjung [38], Wu and Xiao [47], and Cangiani et al [20] which can be viewed as different forms of assumption on the mesh to resolve the geometry. Lemma 2.3 does not require the locally starshaped assumption on the interface and is robust with respect to small variations of the interface as long as the interface deviation is the same.
Notice that if , for , where , we have with . Thus for some constant independent of and . This motivates us to make the following assumption in the remainder of this paper which can be easily satisfied for smooth interfaces if the mesh is locally refined near the interface.
Assumption (H2) For any , .
Now we introduce some notation for DG methods. Let , where , , and . For , denote by . Then . We denote the set of all sides of interior to , that is, not on the boundary . Finally, we set .
For any
, we fix a unit normal vector
of with the convention that is the unit outer normal to if and is the unit outer normal to if . For any , we define the jump of across aswhere is the trace of on in the direction. We define the normal vector function by .
For any subset and , we use the notation
where is the inner product of and is the inner product of .
The unfitted finite element method is based on the idea of “doubling of unknowns” in Hansbo and Hansbo [31]. We define the unfitted finite element space as
For any , we denote , where is the characteristic function of , . For any , we define the liftings , such that
(2.3) 
Our unfitted finite element method is to find such that
(2.4) 
where the bilinear form , and the functional are given by
(2.5)  
(2.6) 
where is the surface gradient on . For any ,
The interface penalty function is
(2.7) 
where is a fixed constant, , , and the mesh function if and if or .
We remark that our unfitted finite element method (2.4) is the socalled local discontinuous Galerkin (LDG) method in Cockburn and Shu [25] which is different from the interior penalty discontinuous Galerkin (IPDG) method used in [31]. We choose the LDG method because the penalty constant in (2.7) can be any fixed constant, while the corresponding penalty constant in the IPDG method has to be sufficiently large to ensure the stability. We refer to Arnold et al [2] for a review of different DG methods for elliptic equations.
Notice that the last term in the bilinear form (2.5) is not present in [21]. It is included in this paper in order to show the discrete Poincaré inequality for unfitted finite element functions in Lemma 4.2 which is crucial for us to study the condition number of the stiffness matrix. We also remark that penalizes the tangential gradient of the finite element solution, not the normal flux of the solution as in Burman and Hansbo [16], Xiao and Wu [47].
For any , we introduce the DG norm
where and . By Lemma 2.3, it is easy to show that
Moreover, by [21, Theorem 2.1] we know that
where is a constant independent of the mesh sizes, , and the interface deviations for all .
Theorem 2.1.
Proof.
For the sake of completeness, we sketch a proof by using the argument in e.g., Perugia and Schötzau [42], Wu and Xiao [47]. For , let be the Stein extension (cf., e.g., Adams and Fournier [1, Theorem 5.14]) of , which is available for any Lipschitz domains, such that . Let , where is the interpolation operator defined in Babuška and Suri [5, Lemma 4.5]. For any , it satisfies that for any ,
(2.8) 
where the constant is independent of , but may depend on . By the multiplicative trace inequality, we have
For any , by Xiao et al [48, Lemma 3.1], [21, Lemma 2.6], we have that for ,
(2.9) 
Thus we obtain by using (2.8) that for any , ,
This implies easily that
(2.10) 
On the other hand, since , we use (2.4) to conclude that
This completes the proof by (2.10) and the triangle inequality. ∎
3 The merging algorithm
In this section, we construct a merging algorithm for the admissible chain of interface elements so that each small interface element in the chain is included in some macroelement which is a large element. We first introduce the concept of admissible chain in §3.1 and five types of patterns of merging small interface elements with their surrounding elements in §3.2. We propose our merging algorithm and prove its reliability in §3.3.
3.1 The admissible chain of interface elements
A chain of interface elements orderly consists of interface elements , , such that is a continuous curve, . We call the length of and denote , .
For any element , we call a neighboring element of if and share a common side, and a diagonal element of if and only share one common vertex. Set , and for , denote , that is, is the set of all th layer elements surrounding , . Obviously, for any .
Definiton 3.1.
A chain of interface elements is called admissible if the following rules are satisfied.

For any , all elements in have the same size as that of .

If has a side such that , then must be a side of some neighboring element , .

Any elements can be neighboring at most two elements in .

For any , the interface elements in must be connected in the sense that the interior of the closed set is a connected domain.
We remark that the four rules of the admissible chains can be easily satisfied if the mesh is well refined near the interface. The purpose of Rules and is to exclude the situations illustrated in Fig.3.1, in which refinements are required to resolve the geometry of the interface. By Rule , the two cases illustrated in Fig.3.2 are not allowed since the interface elements in are not connected, where is the dark element.
3.2 The patterns
Since the interface intersects the boundary of twice at different sides (including the end points), the interface intersects any element only in four possible ways as shown in Fig.3.3. We denote the set of interface elements shown in Fig.3.3(a), the set of interface elements shown in Fig.3.3(b) and (c), and the set of interface elements shown in Fig.3.3(d). By Definition 2.1, each element in is a large element. Thus we only need to consider the merging of type and elements.
A pattern is a set of interface elements and their neighboring and diagonal elements whose union consists of a macroelement. We introduce five types of patterns according to the combination of different types of interface elements, which will be used in our merging algorithm for the admissible chain of interface elements. In the following, for any , stands for its length of the side of which is parallel to the axis, .
Pattern 1: has two neighboring elements , see Fig.3.4. and are respectively the thick part of the sides of and in the figure. We use Algorithm 1 to obtain the macroelements , , and . Here for any closed set , stands for the interior of .
Algorithm 1: Pattern 1
Input:
Output:
if , , and are large elements then
, , ;
else
if and then
let ;
else if and then
let ;
else if and then
let ;