An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems

In this paper we propose a method for the ﬁnite element solution of elliptic interface problem, using an approach due to Nitsche. The method allows for discontinuities, internal to the elements, in the approximation across the interface. We show that optimal order of convergence holds without restrictions on the location of the interface relative to the mesh. Further, we derive a posteriori error estimates for the purpose of controlling functionals of the error and present some numerical examples.


Introduction
As a model elliptic interface problem, we consider a stationary heat conduction problem in two dimensions with a conduction coefficient which is discontinuous across a smooth internal interface. When solving such problems numerically using the standard finite element method, one usually takes the discontinuity of the data into account by enforcing mesh lines along the interface. If this is not done, suboptimal convergence behaviour will occur, cf. [1,8].
As a motivation for this work, we also have in mind more complicated, time dependent or non-linear, problems where the interface moves with time or during iteration. In that case, it may be advantageous to use the same mesh on the domain for different, nearby, locations of the interface, since repeated remeshing of the domain to obtain fitted meshes is very costly. We are thus led to study unfitted finite element methods where the interface is allowed to cross the elements.
In this paper, we propose an unfitted finite element method, based on a variant of NitscheÕs method [9], allowing for discontinuities, internal to the elements, in the approximation across the interface. This method is of optimal order; in particular we show second order convergence in L 2 for appropriately modified piecewise linears on a non-degenerate triangulation. We also consider a posteriori error estimates 1 for functionals of the solution, in the spirit of Becker and Rannacher [3], and use these estimates as a basis for adaptively refining the mesh.
Fitted mesh FE methods for elliptic problems with discontinuous coefficients and homogeneous interface conditions are analysed in [1,6,10]. In [2,4,5], problems with inhomogeneous interface conditions are considered.
As for unfitted mesh methods for interface problems, Barrett and Elliott [2] show first order of convergence in energy-norm and interior second order L 2 error estimates for a piecewise linear method based on boundary penalty and numerical integration over approximate domains. This approach is close in spirit to the method to be presented here, but contains a consistency error that we avoid.
An alternative approach is to construct a FEM-basis where the basis functions fulfill homogeneous interface conditions exactly, and to use the continuous bilinear form (without penalty) to define the method. MacKinnon and Carey [8] use linear basis functions with this property in two dimensions and numerical examples of optimal order of convergence are presented. Li et al. [7] analyse the max-norm interpolation error on cartesian grids of a non-conforming piecewise linear method where the basis functions fulfill homogeneous interface conditions exactly, as well as a conforming method based on further subdivision of the triangles which intersect the interface. The latter method is of optimal order in energy norm and numerical examples indicate second order convergence in the max-norm for this method.
An outline of the paper is as follows. In Section 2 we formulate the continuous problem that we aim to solve, in Section 3 we define the numerical method used for the approximation, and in Section 4 we prove the approximation properties of the corresponding finite element spaces. In Section 5 we prove optimal a priori error estimates and in Section 6 we give corresponding a posteriori error estimates that serve as a basis for adaptive mesh refinement. Finally, in Section 7, we give some implementation details and numerical examples.

Problem formulation and preliminaries
Let X be a bounded domain in R 2 , with convex polygonal boundary oX and an internal smooth boundary C dividing X into two open sets X 1 and X 2 . For any sufficiently regular function u in X 1 [ X 2 we define the jump of u on C by ½u :¼ u 1 j C À u 2 j C , where u i ¼ uj Xi is the restriction of u to X i . Conversely, for u i defined in X i we identify the pair fu 1 ; u 2 g with the function u which equals u i on X i . We consider the following stationary heat conduction problem with a discontinuity in the conductivity across C and an inhomogeneous conormal derivative condition on the interface: Here n is the outward pointing unit normal to X 1 and r n v ¼ n Á rv.
For a bounded open connected domain D we shall use standard Sobolev spaces H r ðDÞ with norm k Á k r;D and spaces H r 0 ðDÞ with zero trace on oD.
We assume that f 2 L 2 ðXÞ, g 2 H 1=2 ðCÞ and, for simplicity, that a is constant in X i with a i > 0. The weak form of (1) is as follows: find u 2 H 1 0 ðXÞ such that aðu; vÞ ¼ ðf ; vÞ X þ ðg; vÞ C ; 8v 2 H 1 0 ðXÞ: ð2Þ Here aðu; vÞ ¼ ðaru; rvÞ X is the bilinear form corresponding to the elliptic operator. It is known that this problem has a unique solution which is in H 2 on each subdomain. The following a priori estimate is valid, see [5]: Here and below, C and c denote generic constants.

The approximation
In a standard finite element method, the jump in normal derivative resulting from the continuity of the flux, when a 1 6 ¼ a 2 , can be taken into account by letting C coincide with mesh lines. We will take an alternative approach and solve (1) approximately using piecewise linear finite elements on a family of conforming triangulations T h of X which are independent of the location of the interface C. Instead, we shall allow the approximation to be discontinuous inside elements which intersect the interface.
We will use the following notation for mesh related quantities. Let h K be the diameter of K and h ¼ max K2T h h K . For any element K, let K i ¼ K \ X i denote the part of K in X i . By G h :¼ fK 2 T h : K \ C 6 ¼ ;g we denote the set of elements that are intersected by the interface. For an element K 2 G h , let C K :¼ C \ K be the part of C in K.
We make the following assumptions regarding the mesh and the interface.
A1: We assume that the triangulation is non-degenerate, i.e., where h K is the diameter of K and q K is the diameter of the largest ball contained in K. A2: We assume that C intersects each element boundary oK exactly twice, and each (open) edge at most once. A3: Let C K;h be the straight line segment connecting the points of intersection between C and oK. We assume that C K is a function of length on C K;h ; in local coordinates C K;h ¼ fðn; gÞ : 0 < n < jC K;h j; g ¼ 0g and C K ¼ fðn; gÞ : 0 < n < jC K;h j; g ¼ dðnÞg: Since the curvature of C is bounded, the assumptions A2 and A3 are always fulfilled on sufficiently fine meshes. Thus the assumptions are natural and not very restrictive; they ensure that the curvature of the interface is well resolved by the mesh.
We shall seek a discrete solution U ¼ ðU 1 Note that functions in V h may be discontinuous across C. Since C may intersect two edges of a triangle arbitrarily, the size of the parts K i are not fully characterized by the meshsize parameters. To define the method, we will therefore use the function j ¼ ðj 1 ; j 2 Þ defined on each element by where jKj :¼ meas K. Clearly, 0 6 j i 6 1 and j 1 þ j 2 ¼ 1 so that The method is defined by the variational problem of finding U 2 V h such that where a h ðU ; /Þ :¼ ða i rU i ; r/ i Þ X 1 [X 2 À ð½U ; far n /gÞ C À ðfar n U g; ½/Þ C þ ðk½U ; ½/Þ C with k sufficiently large (see Lemma 5 below), and In this method, the conditions at C are satisfied weakly by means of a variant of NitscheÕs method. With these definitions, we have the following consistency relation. (4) is consistent in the sense that, for u solving (1), Proof. We first note that, for u solving (1), g À far n ug ¼ ðj 1 þ j 2 Þg À far n ug À j 1 ðg À ½ar n uÞ ¼ j 2 g À j 1 a 1 r n u 1 À j 2 a 2 r n u 2 þ j 1 a 1 r n u 1 À j 1 a 2 r n u 2 ¼ j 2 g À a 2 r n u 2 ; and, similarly, g À far n ug ¼ ðj 1 þ j 2 Þg À far n ug þ j 2 ðg À ½ar n uÞ Since ½u ¼ 0, we may use (5) and GreenÕs formula to obtain which is the statement of the lemma. Ã An immediate consequence of Lemma 1 is the condition which we will refer to as Galerkin orthogonality. A FE basis for V h is easily obtained from a standard FE basis on the mesh by the introduction of new basis functions for the elements that intersect C. Thus, we replace each standard basis function living on an element that intersects the interface by two new basis functions, namely its restrictions to X 1 and X 2 , respectively. The collection of basis functions with support in X i is then clearly a basis for V h i , and hence we obtain a basis for V h by the identification w ¼ ðwj X 1 ; wj X 2 Þ. If the interface coincides exactly with an element edge, no new basis functions are introduced on these elements but the approximating functions may still be discontinuous over such an edge. As a consequence, there are six non-zero basis functions on each element that properly intersects C. Further implementation details are considered in Section 7.

Approximation property of V h
Recall that G h denotes the set of elements that are intersected by the interface. We will use the following mesh dependent norms: and ;h;C : We note for future reference that ðu; vÞ C 6 kvk 1=2;h;C kvk À1=2;h;C : ð7Þ To show that functions in V h approximates functions v 2 H 1 0 ðXÞ \ H 2 ðX 1 [ X 2 Þ to the order h in the norm jjj Á jjj, we construct an interpolant of v by nodal interpolants of H 2 -extensions of v 1 and v 2 as follows. Choose extension operators E i : H 2 ðX i Þ ! H 2 ðXÞ such that ðE i wÞj Xi ¼ w and Let I h be the standard nodal interpolation operator and define The following theorem is valid.
Theorem 2. Let I Ã h be an interpolation operator defined as in (9). Then In the proof of this result, we need to estimate the interpolation error at the interface. To that end, we shall use the following variant of a well known trace inequality on a reference element. The crucial fact is that the constant in this inequality is independent of the location of the interface relative to the mesh. Lemma 3. Map a triangle K 2 G h onto the unit reference triangle e K K by an affine map and denote by e C CK K the corresponding image of C K . Under assumptions A1-A3 of Section 3 there exist a constant C, depending on C but independent of the mesh, such that Proof. We start by showing that Recall that C K;h is the straight line connecting the points of intersection between C and the element K and Assume first that dðnÞ > 0. Since the curvature of the interface is bounded, jd 0 ðnÞj 6 CjC K;h j. As the mesh is non-degenerate this implies that on the reference element we may write, using again ðn; gÞ as local coordinates, where jd d 0 ðnÞj 6 C. We now let D denote the domain bounded by e C CK K and e C CK K;h and note that by the divergence theorem, Asd d 0 is bounded, whence (11) follows from (12) using Cauchy-SchwarzÕ inequality. In a general case where d may switch sign, the same argument may be applied for each part between the intersections of e C CK K and e C CK K;h . It remains to show that the first term on the right in (11) is appropriately bounded. To that end we shall map the triangular part e K K t and the quadrilateral part e K K q of e K K onto new reference domains. We may assume that e C CK K;h intersects e K K in ða; 0Þ and in ð0; bÞ, and, by symmetry, that 0 6 a 6 b 6 1. For a ¼ b ¼ 1, the desired trace inequality is valid. For 1=2 < a 6 b < 1, we may map the triangular part e K K t onto the unit reference triangle by a linear map. By the bound from below on a and b, this map is bounded, uniformly in a and b, with uniformly bounded inverse, and hence (13) is valid also in this case. For 1=2 < a < b ¼ 1 the same argument holds, choosing this time e K K t as the triangular part which contains the origin. Assume now that a 6 1=2. Let ðx x;ŷ yÞ ¼ Mðx x;ỹ yÞ ¼ ðỹ y; ð1 À aÞ À1 ðx x þỹ y À 1ÞÞ: Then the image b K K q ¼ Mð e K K q Þ has its corners in ð0; 0Þ, ð1; 0Þ, ð0; 1Þ, b P P ¼ ðb; ð1 À bÞ=ð1 À aÞÞ, and there holds An additional argument is needed to show uniformity in b P P . Since 0 6 a 6 1=2 and a 6 b 6 1, b P P varies in the domain b D D :¼ f0 6x x 6 1=2; 1 Àx x 6ŷ y 6 1g [ f1=2 6x x 6 1; 1 Àx x 6ŷ y 6 2ð1 Àx xÞg as a and b vary. Let We will show that F ð b P P Þ ¼ sup w2H 1 ðK K q Þ F ð b P P ;ŵ wÞ is uniformly bounded. For points b R R and b S S in b D D, assuming without restriction that F ð b R RÞ P F ð b S S Þ, we have for any w that Given > 0 we may chooseŵ w such that I 6 =2. Note that F ð b R R;v vÞ is continuous for fixedv v since the only dependence of b R R lies in the domains of integration. We may thus take j b R R À b S Sj small enough so that II 6 =2. Hence F ð b P P Þ is continuous on the compact set b D D, and thus (14) holds uniformly in b P P . Finally, since M is bounded, uniformly in a and b, with uniformly bounded inverse, (13) follows and the proof is complete. Ã

Proof of Theorem 2. Recall that
By a standard interpolation estimate we obtain Summing over all triangles that intersect X i , it follows by (8) that Next we consider the jumps on the interface. Since the mesh is non-degenerate, it follows from Lemma 3, scaled by the map from the reference triangle, that h À1 K kwk 2 0;C K 6 Cðh À2 K kwk 2 0;K þ kwk 2 1;K Þ; 8w 2 H 1 ðKÞ: Hence it follows, using again a standard interpolation estimate, that Summing the contributions from K 2 G h , we get from (8) that Finally, Lemma 3 applied to r n w and scaling gives h K kr n wk 2 0;C K 6 Cðkwk 2 1;K þ h 2 K kwk 2 2;K Þ; 8w 2 H 2 ðKÞ; whence similar arguments as above yield kr n ðv i À I Ã h;i v i Þk À1=2;h;C 6 Chkv i k 2;Xi : ð17Þ Since j i < 1, the theorem now follows from (15)-(17). Ã 7

A priori error estimates
We will first show coercivity of the discrete form, for which purpose we will need the following inverse inequality.

A posteriori error estimates
In this Section, we prove a posteriori error estimates and formulate adaptive algorithms for the finite element method (4), following Becker and Rannacher [3].
We will consider control of linear functionals jðeÞ of the error, and define the local and global estimators as and We then have the following a posteriori error estimate.

Implementation and numerical examples
Recall that a FE basis for V h is obtained from a standard FE basis on the mesh by replacing each standard basis function living on an element that intersects the interface by two new basis functions, namely its restrictions to X 1 and X 2 , respectively. Both these new basis functions are, however, represented in the implementation using the same nodes from the original triangulation. The points of intersection between the element edges and the interface are not used to represent the new basis function, and the geometry of the interface and the element parts does not come into play until integrating the terms in the bilinear form. One may thus alternatively think of the approximation as being defined on two overlapping meshes formed by triangles covering the subdomains.
In order to implement the discontinuous approximation, we first determined the set G h of triangles intersected by C. For each K 2 G h , we assigned two identical copies K 0 and K 00 . We assumed that K, with nodes fi; j; kg was split by C into a triangle and a quadrilateral; the case of a split into two triangles was handled by creating a quadrilateral with one side of zero length. We assigned to K 0 the triangular part and to K 00 the quadrilateral part. Thus, on K 00 we created a new node i 0 on the far side of C, and on K 0 we created two new nodes, k 0 and j 0 , on the other side of C, see Fig. 1. To ensure continuity across the edges of K 0 and K 00 (away from C), we also checked if the new nodes had already been created by the same process on the neighboring elements. After having completed this process, we thus had two independent meshes, one completely covering X 1 , and the other completely covering X 2 . The elements crossed by C had been doubled, but coincided geometrically.
To numerically evaluate a h ðÁ; ÁÞ and LðÁÞ, we used the following strategy. The triangles that were not crossed by C were handled in the usual way. On the elements crossed by C, we used centroid quadrature to evaluate all terms, both on the quadrilateral side and on the triangular side. On the interface, we used twopoint Gaussian quadrature. All contributions were then assembled using the old and new nodes defined by the splitting process. We emphasize that the new nodes i 0 , j 0 , and k 0 are to be considered convenient support points for the definition of a continuous, piecewise linear, approximation rather than nodes in the standard finite element sense. The solution at these points is computed but is outside the domain of interest.

Example 1
We considered solutions to the ordinary differential equation The domain is ð0; 1Þ, with an interface at x ¼ 1=2. While this is a one-dimensional problem, we solved it numerically in 2D on the domain ð0; 1Þ Â ð0; 1Þ, with zero Neumann boundary conditions at y ¼ 0 and y ¼ 1. The equation has a closed-form solution, given by We chose a 1 ¼ 1=2, a 2 ¼ 3 and performed a numerical convergence test for different approaches: a standard unfitted FE method, a fitted standard FE method, and the proposed unfitted method. The convergence results are given in Fig. 2, and show the suboptimal behaviour of a standard unfitted method, and the small difference in computational error between the proposed method and a fitted method. We remark that an exact correspondence cannot be obtained, since the meshes will not be exactly the same for the fitted and unfitted methods. In Fig. 3

Example 2
Here we considered a less trivial two-dimensional example (from [7]). The exact solution is given by where r :¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi , and, on the domain ð0; 1Þ Â ð0; 1Þ, we chose r 0 ¼ 3=4, a 1 ¼ 1, and a 2 ¼ 1000. The boundary conditions were symmetry boundaries at x ¼ 0 and y ¼ 0 and Dirichlet boundary conditions corresponding to the exact solution at x ¼ 1 and y ¼ 1. The right-hand side yielding this exact solution is f ¼ À4. In Fig. 4, we give the elevation of the approximate solution on the last mesh in a sequence. The corresponding L 2 -and maximum norm convergence is given in Fig. 5. We note in particular that the order of convergence in maximum norm seems to be approximately equal to the (second order) L 2 -norm convergence, though we have no theoretical results to back this up.

Example 3
The third example was solved on the domain ð0; 1Þ Â ð0; 1Þ, with zero Dirichlet boundary conditions. Centred in the domain is an ellipsoidal inclusion with conduction parameter a min ¼ 1; outside of the ellipse we set a max ¼ 6. In Fig. 6 we show the first mesh and associated solution. We used an adaptive algorithm corresponding to control of the L 2 ðXÞ-error, but made no attempt to tune the interpolation constants or solve a dual problem. Thus, the example only gives an indication of how the adapted meshes appear in such a case. For each successive mesh we refined the third of the elements containing the highest element contribution to the total error as defined by (29). In Fig. 7 we give the last adapted mesh in a sequence, together with an elevation of the corresponding solution.

Concluding remarks
In this paper, we have introduced and analysed a new method for elliptic interface problems on unfitted meshes, i.e., meshes which are independent of the location of the interface. Unlike the standard unfitted finite element method, the proposed approach leads to optimal convergence rates. This have been shown in the model situation of piecewise linear approximations in two dimensions. In future work, we will address the general situation of higher order polynomial approximations in three dimensions. Other extensions under consideration include Stefan problems and fictitious domain type simulations.