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 body-fitted 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 body-fitted shape regular finite element meshes may be difficult and time-consuming, which is the main driving force of the study of unfitted finite element methods. In this paper we will show that the shape regular body-fitted 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 body-fitted shape regular mesh is available, the construction of high-order 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., ) but also essential in controlling the condition number of the finite element stiffness matrix for elliptic equations (see, e.g. ).
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
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  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 . This unfitted finite element method can also be viewed as the interior penalty discontinuous Galerkin method (see, e.g., ) defined on meshes allowing curve-shaped elements. The main difficulty in using the unfitted finite element methods is the so-called small cut cell problem: the cut cells can be arbitrarily small and anisotropic, which can make the stiffness matrix extremely ill-conditioned, especially for high-order 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 , 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 macro-elements have enough support. While the DG formulation is still used in [33, 32, 21], the aggregated unfitted finite element method in  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  an adaptive high-order unfitted finite element method is proposed for elliptic interface problems in which the
a priori and a posteriori error estimates are derived based on noveldomain 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 macro-elements, which are the union of small interface elements and their surrounding elements, are assumed to be rectangular in . This assumption is different from those in [33, 32, 9], see Fig.1.2. The macro-elements 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  are shown without considering the dependence on the finite element approximation order . The assumption that the macro-elements should be rectangular in  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 macro-elements. 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 macro-elements 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 body-fitted finite element meshes for arbitrarily shaped smooth interfaces. To the authors’ best knowledge, this algorithm introduces a new way to generate body-fitted 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 high-order unfitted finite element methods which are known to be of the order in the literature [16, 33, 32, 9, 6] on quasi-uniform meshes with the mesh size . For high order methods, it is known  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 Gauss-Lobatto points on elements not intersecting with the interface. For the interface elements, extra care must be taken as the basis of the spectral element onis ill-conditioned 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  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  including the domain inverse estimates and the concept of the interface deviation can be extended to study three-dimensional interface problems. The merging algorithm in the three-dimensional case is more challenging. Nevertheless, we believe that with the new insights gained in this paper for the two-dimensional case, reliable algorithms for constructing cubic macro-elements 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 three-dimensional 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].
(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.
Assumption (H1) For each , there exists a rectangular macro-element 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 macro-element 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 . We remark that there are different definitions of shape regular polygons in the literature, see, e.g., Ming and Shi  and Brenner and Sung .
The following concept of interface deviation is introduced in .
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  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].
Let be a triangle with vertices , , , where . Let and . Then we have
The proof of this lemma makes use of the following one-dimensional domain inverse estimate in [21, Lemma 2.3]
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 .
Let and , Then for , we have
where the constant is independent of , , and .
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)].
Let and , Then for , we have
where the constant is independent of and .
We remark that inverse estimates on star-shaped curve elements are studied in Massjung , Wu and Xiao , and Cangiani et al  which can be viewed as different forms of assumption on the mesh to resolve the geometry. Lemma 2.3 does not require the locally star-shaped 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 .
, we fix a unit normal vectorof 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 as
where 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 . 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
Our unfitted finite element method is to find such that
where the bilinear form , and the functional are given by
where is the surface gradient on . For any ,
The interface penalty function is
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 so-called local discontinuous Galerkin (LDG) method in Cockburn and Shu  which is different from the interior penalty discontinuous Galerkin (IPDG) method used in . 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  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 . 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 , Xiao and Wu .
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 .
For the sake of completeness, we sketch a proof by using the argument in e.g., Perugia and Schötzau , Wu and Xiao . 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 ,
where the constant is independent of , but may depend on . By the multiplicative trace inequality, we have
Thus we obtain by using (2.8) that for any , ,
This implies easily that
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 macro-element 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 .
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 macro-element. 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 macro-elements , , and . Here for any closed set , stands for the interior of .
Algorithm 1: Pattern 1
if , , and are large elements then
, , ;
if and then
else if and then
else if and then