DeepAI

# An Arbitrarily High Order Unfitted Finite Element Method for Elliptic Interface Problems with Automatic Mesh Generation

We consider the reliable implementation of high-order unfitted finite element methods on Cartesian meshes with hanging nodes for elliptic interface problems. We construct a reliable algorithm to merge small interface elements with their surrounding elements to automatically generate the finite element mesh whose elements are large with respect to both domains. We propose new basis functions for the interface elements to control the growth of the condition number of the stiffness matrix in terms of the finite element approximation order, the number of elements of the mesh, and the interface deviation which quantifies the mesh resolution of the geometry of the interface. Numerical examples are presented to illustrate the competitive performance of the method.

• 6 publications
• 175 publications
09/03/2020

### An adaptive high-order unfitted finite element method for elliptic interface problems

We design an adaptive unfitted finite element method on the Cartesian me...
05/12/2020

### Finite Element Methods For Interface Problems On Local Anisotropic Fitting Mixed Meshes

A simple and efficient interface-fitted mesh generation algorithm is dev...
08/14/2019

### Stability of explicit Runge-Kutta methods for high order finite element approximation of linear parabolic equations

We study the stability of explicit Runge-Kutta methods for high order La...
05/24/2022

### Accelerating High-Order Mesh Optimization Using Finite Element Partial Assembly on GPUs

In this paper we present a new GPU-oriented mesh optimization method bas...
09/10/2019

### Recurrent Neural Networks as Optimal Mesh Refinement Strategies

We show that an optimal finite element mesh refinement algorithm for a p...
03/14/2022

### Extended Finite Elements for 3D-1D coupled problems via a PDE-constrained optimization approach

In this work we propose the application of the eXtended Finite Element M...
11/14/2017

### Generation of unstructured meshes in 2-D, 3-D, and spherical geometries with embedded high resolution sub-regions

We present 2-D, 3-D, and spherical mesh generators for the Finite Elemen...

## 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 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., [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

 −div(a∇u)=f  in Ω1∪Ω2, (1.1) [[u]]Γ=0,[[a∇u⋅n]]Γ=0  on Γ,  u=g  on ∂Ω, (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 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 [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 macro-elements 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 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 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 macro-elements, 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 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 [9] are shown without considering the dependence on the finite element approximation order . The assumption that the macro-elements 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 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 [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 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 on

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

###### 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 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

 M={M(K):K∈TΓ}∪{K∈T:K⊄M(K′) for some K′∈TΓ}.

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 ,

 ηiK=distH(ΓK,ΓhK)dist(AiK,ΓhK).

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

 ∥v∥L2(Δ)≤T(1+δa−121−δa−12)2p+3/2∥v∥L2(Δδ)  ∀v∈Qp(Δ),

where .

The proof of this lemma makes use of the following one-dimensional domain inverse estimate in [21, Lemma 2.3]

 ∥g∥2L2(Iλ∖¯I)≤12[(λ+√λ2−1)2p+1−1]∥g∥2L2(I)  ∀g∈Qp(Iλ), (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

 ∥vi∥L2(Kh−δKi)≤∥vi∥L2(Ki)≤CT(1+3ηK1−ηK)2p+3/2∥vi∥L2(Kh−δKi)  ∀v∈Qp(K), (2.2)

where the constant is independent of , , and .

###### Proof.

The left inequality (2.2) is trivial since . Here we prove the right inequality in (2.2) when intersects at neighboring sides. The other cases can be proved similarly.

We use the notation in Fig.2.1 in which , are parallel to and the distances of to are . Then and is the polygon bounded by sides of and . Let , . By definition, the interface deviation , . By Lemma 2.1, for any ,

 ∥v∥L2(K1)≤∥v∥L2(ΔA1KB′′C′′) ≤ T(1+2δ/(d1+δ)1−2δ/(d1+δ))2p+3/2∥v∥L2(ΔA1KB′C′) ≤ T(1+3ηK1−ηK)2p+3/2∥v∥L2(Kh−δK1).

The case for can be proved similarly. This completes the 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

 aK={a1+a22if K∈MΓ,aiif K∈Ωi.,ΘK=⎧⎨⎩T(1+3ηK1−ηK)4p+3if K∈MΓ,1otherwise.

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

 ∥∇v∥L2(Ki)≤Cp2h−1KΘ1/2K∥v∥L2(Ki)  ∀v∈Qp(K), ∥v∥L2(∂Ki)≤Cph−1/2KΘ1/2K∥v∥L2(Ki)  ∀v∈Qp(K),

where the constant is independent of and .

We remark that inverse estimates on star-shaped 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 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 .

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 as

 [[v]]e:=v−−v+  ∀e∈Eside∪EΓ,    [[v]]e:=v−  ∀e∈Ebdy,

where is the trace of on in the direction. We define the normal vector function by .

For any subset and , we use the notation

 (u,v)ˆM:=∑K∈ˆM(u,v)K,  ⟨u,v⟩^E:=∑e⊂^E⟨u,v⟩e,

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

 Xp(M)={v1χΩ1+v2χΩ2:vi|K∈Qp(K),K∈M,i=1,2}.

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

 ah(U,v)=Fh(v)    ∀v∈Xp(M), (2.4)

where the bilinear form , and the functional are given by

 ah(v,w)= (a(∇hv−L(v)),∇hw−L(w))M+⟨α[[v]],[[w]]⟩¯E+⟨p−2h∇T[[v]],∇T[[w]]⟩EΓ, (2.5) Fh(v)= (f,v)M−(aL1(g),∇hv−L(v))M+⟨αg,v⟩Ebdy, (2.6)

where is the surface gradient on . For any ,

The interface penalty function is

 α|e=α0aeΘeh−1ep2  ∀e∈E, (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 so-called 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

 ∥v∥2DG:=∥a1/2∇v∥2M+∥α1/2[[v]]∥2¯E+∥p−1h1/2∇T[[v]]∥2EΓ,

where and . By Lemma 2.3, it is easy to show that

 ah(v,v)≤C∥v∥2DG   ∀v∈Xp(M).

Moreover, by [21, Theorem 2.1] we know that

 ah(v,v)≥cstab∥v∥2DG   ∀v∈Xp(M),

where is a constant independent of the mesh sizes, , and the interface deviations for all .

###### Theorem 2.1.

Let the solution of the problem (1.1)-(1.2) , , and be the solution of the problem (2.4). Then we have

 ∥u−U∥DG≤CΘ1/2hmin(p+1,k)−1pk−3/2∥u∥Hk(Ω1∪Ω2),

where , , and the constant is independent of the mesh sizes, , and the interface deviations for all .

###### 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 ,

 ∥w−Ihp(w)∥Hj(K)≤Chmin(p+1,k)−jKpk−j∥v∥Hk(K)  ∀v∈Hk(K), (2.8)

where the constant is independent of , but may depend on . By the multiplicative trace inequality, we have

 ∥w∥L2(∂K)≤Ch−1/2K∥w∥L2(K)+C∥w∥1/2L2(K)∥∇w∥1/2L2(K)  ∀w∈H1(K).

For any , by Xiao et al [48, Lemma 3.1], [21, Lemma 2.6], we have that for ,

 ∥w∥L2(ΓK)≤C∥w∥1/2L2(Ki)∥∇w∥1/2L2(Ki)+∥w∥L2(∂Ki∖¯ΓK)  ∀w∈H1(K). (2.9)

Thus we obtain by using (2.8) that for any , ,

 ∥w−Ihp(w)∥Hj(∂Ki)≤Chmin(p+1,k)−j−1/2pk−j−1/2∥w∥Hk(K)  ∀w∈Hk(K).

This implies easily that

 ∥u−uI∥DG≤CΘ1/2hmin(p+1,k)−1pk−3/2∥u∥Hk(Ω1∪Ω2). (2.10)

On the other hand, since , we use (2.4) to conclude that

 ∥uI−U∥2DG≤c−1stabah(uI−U,uI−U) = c−1stabah(uI−u,uI−U) ≤ C∥uI−u∥DG∥uI−U∥DG.

This completes the proof by (2.10) and the triangle inequality. ∎

To conclude this section, we remark that the same a posteriori error estimate in [21, Theorem 3.1] also holds for the solution in (2.4). Here we omit the details.

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

###### 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 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

Input:

Output:

if , , and are large elements then

, , ;

else

if and then

let ;

else if and then

let ;

else if and then

let ;