High-Order Finite Element Methods for Interface Problems: Theory and Implementations

We propose arbitrary order extended finite element (XFE) methods for solving both two and three dimensional elliptic interface problems. The meshes in our methods do not need to fit the interface and discontinuous Galerkin (DG) schemes are used in the vicinity of the interface to incorporate the transmission conditions between subdomains. Optimal error estimates in the piecewise H1-norm and in the L2-norm are rigorously proved for the scheme. Two implementation aspects are addressed in this paper: (1) An optimal multigrid solver is introduced for the generated linear system, which converges uniformly with respect to the mesh size. (2) A high order algorithm is presented to compute the integral over elements intersecting the interface. Numerical experiments are reported to support the theoretical results.


Introduction
The interface problems which involve partial differential equations having discontinuous coefficients across certain interfaces are often encountered in fluid dynamics, electromagnetics and materials science. Because of the low global regularity and the irregular geometry of the interface, the standard numerical methods which are efficient for smooth solutions usually lead to loss in accuracy across the interface.
For arbitrarily shaped interface , it is known that optimal or nearly optimal convergence rate can be recovered if body-fitted finite element meshes are used, see e.g. [6,8,20,29]. Here, by "body-fitted meshes" we mean an element of the underlying mesh is required to intersect with the interface only through its boundaries (Fig. 1). Unfortunately, when the geometry is complex, this usually leads to a nontrivial interface meshing problem. Therefore, numerous modified finite difference methods based only on simple Cartesian grids have been proposed in the literature. We refer to the immersed boundary method [24], the immersed interface method [17,18], the ghost fluid method [21], and the references therein. In the Fig. 1 A body-fitted, shape regular mesh finite element setting, we refer to the work of the immersed finite element method [7,11,19], the multiscale finite element method [9], the penalty finite element method [1].
In the past decade, a combination of the extended finite element method (XFEM) with the Nitsche scheme has become a popular discretization method. As the first attempt, an unfitted finite element method was proposed in [13] which can be viewed as a linear and consistent modification of [1]. This approach has motivated a number of works, e.g., the unfitted finite element method [4,5,12], the Ghost penalty method [2,3], the unfitted discontinuous Galerkin methods [22]. Although significant progresses in the error analyses of some methods have been made, the development of high-order accurate unfitted FEMs with rigorous error analysis is still challenging. We refer to the work of [14-16, 22, 27, 28] which claim high order approximations. In [22], an hp-unfitted discontinuous Galerkin method for Problem (1) was considered, and optimal h-convergence for arbitrary p was shown for the two-dimensional case in the energy norm and in the L 2 -norm. With an extra flux penalty term applied on the interface, [27] gave better hp a priori error estimates in both two and three dimensions. In [15,16], an isoparametric finite element method with a high order geometrical approximation of level set domains was presented. The analysis reveals optimal order error bounds with respect to h for the geometry approximation and for the finite element approximation. In [14,28], various issues related to unfitted methods was addressed, including the dependence of error estimates on the diffusion coefficients, the condition number of the discrete system, and the choice of stabilization parameters.
The Nitsche-XFEM can be interpreted as applying interior penalty (IP) methods on the interface, and our method falls into this category. The major step in our variant is an appropriate choice of the mesh and geometry dependent weights in the average (see (6)), which lead to trace and inverse inequalities for possibly degenerated subelements (see (9)). We note that in our approach, the penalization is applied only to the jump of the solution values across the interface (compared with the bilinear form in [27]). The optimal h-convergence rate for arbitrary high-order discretization in the energy and L 2 -norm are proved regardless of the dimension. We refer to [14][15][16] for the similar estimates with respect to h and [27] for a refined version with respect to both h and p.
Efficient implementations of this method are then discussed in two aspects. We first consider an optimal multigrid solver for the generated linear system. We use the continuous FE space as a "background" subspace, with some smoothing operations added near the interface, to formulate a nested geometrical multigrid method. We prove the optimality of this special multiplicative multigrid method, which means the method converges uniformly with respect to the mesh size, and is independent of the location of the interface relative to the meshes. Since the assembling of the stiffness matrix will require integration over curved surfaces and volumes, we then implement a robust and arbitrarily high order numerical quadrature algorithm by transforming surface and volume integrals into multiple 1-D integrals. The code for the algorithm is freely available in the open source finite element toolbox Parallel Hierarchical Grid (PHG) [26]. We also refer to [23,25] for different approaches to compute integrals on curved sub-elements and their curved boundaries.
The layout of this paper is as follows. In Sect. 2 we introduce the XFE spaces and reformulate the interface problem (1) in DG schemes. The H 1 -and L 2 -error estimates of both schemes-which attain the optimal order of the convergence rate in respect to mesh size h-are given. In Sect. 3, we give an optimal multigrid method for the aforementioned DG-XFE schemes. Numerical examples for both two and three dimensions are reported in Sect. 4, to illustrate the high accuracy of the algorithm.

XFE and DG Schemes for Interface Problems
We consider the following elliptic interface problem for u: Let = 1 ∪ ∪ 2 be a bounded and convex polygonal or polyhedral domain in R d , d = 2 or 3, where 1 and 2 are two subdomains of and are separated by a C 2 -smooth interface (see Fig. 2 for an illustration of a unit square that contains a circle as an interface), Here α(x) = α i , i = 1, 2, is a piecewise constant function on the partition 1 ∪ 2 . Denote by {T h }, a family of conforming, quasi-uniform, and regular partitions of into triangles and parallelograms/tetrahedrons and parallelepipeds. As K is of regular shape, there is a constant γ 0 such that We define the set of all elements intersected by as T h = {K ∈ T h : |K ∩ | = 0}. Each T h induces a partition of interface , which we denote by and n i be the unit outward normal vector on ∂K i with i = 1, 2. As is of class C 2 , it is easy to prove that (cf. [6,31]) each interface segment/patch e K is contained in a strip of width δ and satisfies We define the weighted average {·} and the jump [·] on e ∈ E h by For the stability analysis of our schemes, we define (κ 1 , κ 2 ) on each element as follows: Clearly, 0 ≤ κ i ≤ 1 and κ 1 + κ 2 = 1 so that {·} is a convex combination along . Roughly speaking, we adopt the weight κ i = |K i | |K| suggested in [13] for general subelements and we set κ i = 0 for |K i | < ch d+1 K . Here, the user-defined constant c 0 ≥ 2γ 0 γ 1 and γ 0 , γ 1 are constants defined in (2) and (3), respectively. The dependence of c 0 on these generic constants is elaborated in Lemma 1.
Let χ i be the characteristic function on i with i = 1, 2. Given a mesh T h , let V h be the continuous piecewise polynomial function space of degree p ≥ 1 on the mesh.
Then, the DG-XFE method for the interface problem is: where For η β sufficiently large, the norm corresponding to the bilinear form B h (·, ·) is uniformly equivalent to · B h , which is defined by The crucial component in regard to establishing this equivalence result and also the stability of bilinear forms is the control on the weighted normal derivatives, which is stated as a trace and inverse inequality in Lemma 1.

Lemma 1 ([27, 28])
Let γ 0 and γ 1 be constants defined in (2) and (3), respectively. If we choose c 0 ≥ 2γ 0 γ 1 in the definition (6) of κ, there exists a positive constant h 0 such that for all h ∈ (0, h 0 ] and any interface segment/patch e K = K ∩ ∈ E h , the following estimates hold on both sub-elements of K: The coercivity and boundedness of B h (·, ·) in its norm · 2 B h is then a direct consequence of the Cauchy-Schwarz inequality. and provided the penalty parameter η β is chosen sufficiently large.
The XFE space has optimal approximation quality for piecewise smooth functions in H p ( 1 ∪ 2 ). The following theorem is proved in [28] as an analogue of Cea's lemma. (1) satisfies u ∈ H s ( 1 ∪ 2 ), where s ≥ 2 is an integer. Let μ = min{p + 1, s}. The following error estimates hold for any h ∈ (0, h 0 ]: If η β is chosen sufficiently large (see (11)) and u h is the solution to the first scheme of (7), then

Theorem 1 Assume that the interface is C 2 smooth and that the solution of the elliptic interface problem
The hidden constants in the above estimates are dependent on the angle condition of the mesh T h , the degree of the polynomials, the parameter in the scheme, and α(x), but are independent of the location of the interface relative to the mesh. Here, the constant h 0 is from Lemma 1.

An Optimal Multigrid Method for (7)
In this section, we propose a two-level geometric multigrid solver of the finite element problem (7). It is well known that the element K with a "small" cut (i.e. |K ∩ i |/|K| 1) would have adverse effect on the conditioning of the resulting stiffness matrices (see e.g. [3]). Our approach is based on the general theory of the successive subspace correction (SSC) method of solving on a linear vector spacẽ V = J i=0 V i with inner product (·, ·) the equation (Au, v) = (f, v), where A :Ṽ →Ṽ is a symmetric positive definite operator.
We apply SSC for a relatively simple case of two subspaces (i.e. J = 2), that is, is the space of nodal basis functions that vanish on N h := {x j : |supp(ψ j ) ∩ | = 0}. With a slight abuse of notation, the DG-XFE scheme induces a symmetric positive definite operator B h for β = 1. LetB h andB h be the restrictions of B h on V 0 h and V h , respectively. LetR h : V 0 h → V 0 h be approximately an inverse ofB h . We have this two-level successive subspace correction method (Algorithm 1). The similar idea has been employed in a special linear case in [32] and analyzed using the framework given in [30,33].

Algorithm 1 The multigrid method for (7)
Implement this iterative procedure until converge: 1. do subspace correction on V 0 h with an inexact solverR h ; 2. do subspace correction on V h with an exact solver (B h ) −1 .
Obviously, Algorithm 1 defines an iterative method for solving B h u h = f h . Denote by A h the iterator of the method, then the error contract property is summarized as the following theorem.

Theorem 2 ([28]) Assume that I −R hBh B h ρ < 1. Then Algorithm 1 is uniformly convergent with respect to the mesh size with
where is a constant independent of h.
When T h is a shape-regular grid with a geometrical multilevel structure, then a geometric multigrid process can be implemented on V 0 h , and the approximate inverseR h ofB h can be chosen to be the iterator of V -cycle multigrid method.

Numerical Tests
In this section, we present some initial results to demonstrate the high-order accuracy and robustness of our method.

High-Order Numerical Quadratures on "cut" Elements
Assembling the local stiffness matrix and the corresponding RHS for K ∈ T h requires integration over irregularly shaped manifolds: where is defined by the zero level set of a piecewise smooth function. Our implementation of (13) relies on a general-purpose and arbitrarily high order numerical quadrature algorithm proposed in [10]. The basic idea is to choose a local coordinate system with three orthogonal directions, decompose integrals in (13) into multiple 1-D integrals along these directions, and use 1-D Gaussian quadratures to compute these integrals. For 1-D Gaussian quadratures to work, the local coordinate system should be suitably chosen according to properties of K and to prevent essential singularities from appearing in the 1-D integrands, and the integration intervals are divided into subintervals at the non essential singularities of the integrands. We note that the proposed algorithm only requires finding roots of univariate nonlinear functions in given intervals and evaluating the integrand, the level set function, and the gradient of the level set function at given points. It can achieve arbitrarily high order by increasing the orders of Gaussian quadratures, and does not need extra a priori knowledge about the integrand and the level set function.
This algorithm has been implemented in the file src/quad-interface.c and include/phg/quad-interface.h in PHG [26]. Extensive h− and p−convergence tests have been performed in [10] and included in a sample code test/quad_test2.c.

2-D Numerical Examples
Let domain be the unit square (0, 1) 2 and interface be the zero level set of the function ϕ(x) = (x 1 −0.5) 2 +(x 2 −0.5) 2 −1/7. The subdomain 1 is characterized by ϕ(x) < 0 and 2 by ϕ(x) > 0. The domain is partitioned into grids of squares with the same size h. The exact solution is chosen as The right-hand side can be computed accordingly.
We implement Algorithm 1, with V -cycle geometric multigrid based on the unfitted grid T h playing as the coarse grid corrector. In each pre-and post-smoothing stage of V -cycle iterator, we perform Gauss-Seidel for two times. We record the numerical results in Table 1. In these examples, the initial guess is 0, and the stopping criterion is From Table 1, we can see that the multigrid method converges uniformly with respect to the mesh size, which confirms our theoretical results.

3-D Numerical Examples
The settings of this numerical experiment are as follows. The domain = (0, 1) 3 . The interfaces are two touched spheres of radius 0.1 centered at (0.4, 0.5, 0.5) and (0.6, 0.5, 0.5). The exact solution is given by The discontinuous coefficient function is defined such that α 1 = 1 and α 2 = 100. A convergence study is performed on a series of meshes generated by uniform refinements of an initial mesh consisting of 6 congruent tetrahedra. Relative errors and convergence rates of numerical solutions for P p elements for p = 1, 2, 3 and 4 are listed in Table 2, with the quadrature order q = 2p + 3. The convergence rates are optimal for both H 1 ( )-errors (order p) and L 2 ( )-errors (order p + 1). For the time being, however, the design of multigrid solver for 3-D case is still on-going. The computations for P 1 , P 2 and P 3 elements were done using the 64-bit double precision and the linear systems were solved using MUMPS, but for P 4 element, to eliminate influences of roundoff errors, the computations were done using the 80bit extended double precision and the linear systems were solved using the GMRES method with MUMPS in double precision as its preconditioner. The performance of Algorithm 1 will be reported in a future work.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made. The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.