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

Anita Hansbo, Peter Hansbo

► To cite this version:

Anita Hansbo, Peter Hansbo. An unfitted finite element method, based on Nitsche's method, for elliptic interface problems. Computer Methods in Applied Mechanics and Engineering, 2002, 191 (47-48), pp. 5537-5552. 10.1016/S0045-7825(02)00524-8. hal-01352903

HAL Id: hal-01352903

https://hal.science/hal-01352903

Submitted on 10 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An unfitted finite element method, based on Nitsche's method, for elliptic interface problems

Anita Hansbo<sup>a</sup>, Peter Hansbo<sup>b</sup>

<sup>a</sup> Department of Informatics and Mathematics, University of Trollhättan-Uddevalla, Box 957, S-461 39 Trollhättan, Sweden

<sup>b</sup> Department of Applied Mechanics, Chalmers University of Technology, S-412 96 Göteborg, Sweden

In this paper we propose a method for the finite 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.

## 1. 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, sub-optimal 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 linear on a non-degenerate triangulation. We also consider a posteriori error estimates
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.

## 2. Problem formulation and preliminaries

Let $\Omega$ be a bounded domain in $\mathbb{R}^2$, with convex polygonal boundary $\partial\Omega$ and an internal smooth boundary $\Gamma$ dividing $\Omega$ into two open sets $\Omega_1$ and $\Omega_2$. For any sufficiently regular function $u$ in $\Omega_1 \cup \Omega_2$ we define the jump of $u$ on $\Gamma$ by $[u] := u_1|_{\Gamma} - u_2|_{\Gamma}$, where $u_1 = u|_{\partial\Omega_1}$ is the restriction of $u$ to $\Omega_1$. Conversely, for $u_i$ defined in $\Omega_i$ we identify the pair $\{u_1, u_2\}$ with the function $u$ which equals $u_i$ on $\Omega_i$. We consider the following stationary heat conduction problem with a discontinuity in the conductivity across $\Gamma$ and an inhomogeneous conormal derivative condition on the interface:

$$ 
\begin{array}{l}
-\nabla \cdot (\alpha \nabla u) = f \quad \text{in } \Omega_1 \cup \Omega_2, \\
u = 0 \quad \text{on } \partial\Omega, \\
[u] = 0 \quad \text{on } \Gamma, \\
[\alpha \nabla_n u] = g \quad \text{on } \Gamma.
\end{array}
\qquad (1) $$

Here $n$ is the outward pointing unit normal to $\Omega_1$ and $\nabla_n v = n \cdot \nabla v$.

For a bounded open connected domain $D$ we shall use standard Sobolev spaces $H^r(D)$ with norm $\|\cdot\|_{r,D}$ and spaces $H_0^k(D)$ with zero trace on $\partial D$. The inner products in $H^0(D) = L_2(D)$ is denoted $(\cdot, \cdot)_D$. For a bounded open set $G = \cup_{i=1}^2 D_i$, where $D_i$ are open mutually disjoint components of $G$, we let $H^k(D_1 \cup D_2)$ denote the Sobolev space of functions in $G$ such that $u|_{D_i} \in H^k(D_i)$ with norm

$$ \| \cdot \|_{k, D_1 \cup D_2} = \sum_{i=1}^{2} \| \cdot \|_{k, D_i}^2 \quad (1/2)^{1/2}. $$
We assume that $f \in L_2(\Omega)$, $g \in H^{1/2}(\Gamma)$ and, for simplicity, that $\alpha$ is constant in $\Omega_i$ with $\alpha_i > 0$. The weak form of (1) is as follows: find $u \in H_0^1(\Omega)$ such that

$$a(u, v) = (f, v)_{\Omega} + (g, v)_{\Gamma}, \quad \forall v \in H_0^1(\Omega). \qquad (2)$$

Here

$$a(u, v) = (\alpha \nabla u, \nabla v)_{\Omega}$$

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]:

$$\|u\|_{1,\Omega} + \|u\|_{2,\Omega_1 \cup \Omega_2} \le C(\|f\|_{0,\Omega} + \|g\|_{1/2,\Gamma}). \qquad (3)$$

Here and below, $C$ and $c$ denote generic constants.

### 3. The approximation

In a standard finite element method, the jump in normal derivative resulting from the continuity of the flux, when $\alpha_1 \ne \alpha_2$, can be taken into account by letting $\Gamma$ 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 $\Omega$ which are independent of the location of the interface $\Gamma$. 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_{K \in T_h} h_K$. For any element $K$, let $K_i = K \cap \Omega_i$ denote the part of $K$ in $\Omega_i$. By $G_h := \{K \in T_h : K \cap \Gamma \ne \emptyset\}$ we denote the set of elements that are intersected by the interface. For an element $K \in G_h$, let $\Gamma_K := \Gamma \cap K$ be the part of $\Gamma$ in $K$.

We make the following assumptions regarding the mesh and the interface.

A1: We assume that the triangulation is non-degenerate, i.e.,

$$h_K / \rho_K \le C \quad \forall K \in T_h,$$

where $h_K$ is the diameter of $K$ and $\rho_K$ is the diameter of the largest ball contained in $K$.

A2: We assume that $\Gamma$ intersects each element boundary $\partial K$ exactly twice, and each (open) edge at most once.

A3: Let $\Gamma_{K,h}$ be the straight line segment connecting the points of intersection between $\Gamma$ and $\partial K$. We assume that $\Gamma_K$ is a function of length on $\Gamma_{K,h}$; in local coordinates

$$\Gamma_{K,h} = \{(\xi, \eta) : 0 < \xi < |\Gamma_{K,h}|, \eta = 0\}$$

and

$$\Gamma_K = \{(\xi, \eta) : 0 < \xi < |\Gamma_{K,h}|, \eta = \delta(\xi)\}.$$

Since the curvature of $\Gamma$ 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, U_2)$ in the space $V^h = V_1^h \times V_2^h$, where

$$V_i^h = \{\phi_i \in H^1(\Omega_i) : \phi_i|_{K_i} \text{ is linear, } \phi_i|_{\partial\Omega} = 0\}.$$
Note that functions in $V_h$ may be discontinuous across $\Gamma$. Since $\Gamma$ 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 $\kappa = (\kappa_1, \kappa_2)$ defined on each element by

$$ \kappa_i|_{\kappa} = \frac{|K_i|}{|K|}, $$

where $|K| := \text{meas } K$. Clearly, $0 \le \kappa_i \le 1$ and $\kappa_1 + \kappa_2 = 1$ so that

$$ \{\phi\} := (\kappa_1\phi_1 + \kappa_2\phi_2)|_{\Gamma} $$

is a convex combination of $\phi = (\phi_1, \phi_2)$ along $\Gamma$.

The method is defined by the variational problem of finding $U \in V^h$ such that

$$ a_h(U, \phi) = L(\phi), \quad \forall \phi \in V^h, \qquad (4) $$

where

$$ a_h(U, \phi) := (\alpha_i \nabla U_i, \nabla \phi_i)_{\Omega_1 \cup \Omega_2} - ([U], \{\alpha \nabla_n \phi\})_\Gamma - (\{\alpha \nabla_n U\}, [\phi])_\Gamma + (\lambda[U], [\phi])_\Gamma $$

with $\lambda$ sufficiently large (see Lemma 5 below), and

$$ L(\phi) := (f, \phi)_{\Omega} + (\kappa_2 g, \phi_1)_{\Gamma} + (\kappa_1 g, \phi_2)_{\Gamma}. $$

In this method, the conditions at $\Gamma$ are satisfied weakly by means of a variant of Nitsche's method. With these definitions, we have the following consistency relation.

**Lemma 1.** The discrete problem (4) is consistent in the sense that, for $u$ solving (1),

$$ a_h(u, \phi) = L(\phi), \quad \forall \phi \in V^h. $$

**Proof.** We first note that, for $u$ solving (1),

$$ \begin{aligned} g - \{\alpha \nabla_n u\} &= (\kappa_1 + \kappa_2)g - \{\alpha \nabla_n u\} - \kappa_1(g - [\alpha \nabla_n u]) \\ &= \kappa_2 g - \kappa_1 \alpha_1 \nabla_n u_1 - \kappa_2 \alpha_2 \nabla_n u_2 + \kappa_1 \alpha_1 \nabla_n u_1 - \kappa_1 \alpha_2 \nabla_n u_2 = \kappa_2 g - \alpha_2 \nabla_n u_2, \end{aligned} $$

and, similarly,

$$ \begin{aligned} g - \{\alpha \nabla_n u\} &= (\kappa_1 + \kappa_2)g - \{\alpha \nabla_n u\} + \kappa_2(g - [\alpha \nabla_n u]) \\ &= \kappa_1 g - \kappa_1 \alpha_1 \nabla_n u_1 - \kappa_2 \alpha_2 \nabla_n u_2 - \kappa_2 \alpha_1 \nabla_n u_1 + \kappa_2 \alpha_2 \nabla_n u_2 = (1 + \kappa_2)g - \alpha_1 \nabla_n u_1, \end{aligned} $$

so that

$$ \{\alpha \nabla_n u\} = \alpha_1 \nabla_n u_1 - \kappa_2 g = \alpha_2 \nabla_n u_2 + \kappa_1 g. \qquad (5) $$

Since $[u] = 0$, we may use (5) and Green's formula to obtain

$$ \begin{aligned} a_h(u, \phi) &= (\alpha \nabla u, \nabla \phi)_{\Omega_1 \cup \Omega_2} - (\{\alpha \nabla_n u\}, \phi_1 - \phi_2)_\Gamma \\ &= (\alpha \nabla u, \nabla \phi)_{\Omega_1 \cup \Omega_2} - (\alpha_1 \nabla_n u_1 - \kappa_2 g, \phi_1)_\Gamma + (\alpha_2 \nabla_n u_2 + \kappa_1 g, \phi_2)_\Gamma \\ &= -(\nabla \cdot (\alpha \nabla u), \phi)_{\Omega_1 \cup \Omega_2} + (\kappa_2 g, \phi_1)_\Gamma + (\kappa_1 g, \phi_2)_\Gamma = (f, \phi)_\Omega + (\kappa_2 g, \phi_1)_\Gamma + (\kappa_1 g, \phi_2)_\Gamma = L(\phi), \end{aligned} $$

which is the statement of the lemma. $\square$
An immediate consequence of Lemma 1 is the condition

$$a_h(u - U, \phi) = 0, \quad \forall \phi \in V_h, \qquad (6)$$

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 $\Gamma$. Thus, we replace each standard basis function living on an element that intersects the interface by two new basis functions, namely its restrictions to $\Omega_1$ and $\Omega_2$, respectively. The collection of basis functions with support in $\Omega_i$ is then clearly a basis for $V_i^h$, and hence we obtain a basis for $V_h$ by the identification $\psi = (\psi|_{\Omega_1}, \psi|_{\Omega_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 $\Gamma$. Further implementation details are considered in Section 7.

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

$$\|v\|_{1/2,h,\Gamma}^2 := \sum_{K \in G_h} h_K^{-1} \|v\|_{0,\Gamma_K}^2,$$

$$\|v\|_{-1/2,h,\Gamma}^2 := \sum_{K \in G_h} h_K \|v\|_{0,\Gamma_K}^2,$$

and

$$\|v\|^2 := \|\nabla v\|_{0,\Omega_1 \cup \Omega_2}^2 + \|\{\nabla_n v\}\|_{-1/2,h,\Gamma}^2 + \|[v]\|_{1/2,h,\Gamma}^2.$$

We note for future reference that

$$(u, v)_{\Gamma} \le \|v\|_{1/2,h,\Gamma} \|v\|_{-1/2,h,\Gamma}. \qquad (7)$$

To show that functions in $V_h$ approximates functions $v \in H_0^1(\Omega) \cap H^2(\Omega_1 \cup \Omega_2)$ to the order $h$ in the norm $||| \cdot |||$, 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(\Omega_i) \to H^2(\Omega)$ such that $(E_i w)|_{\Omega_i} = w$ and

$$\|E_i w\|_{s,\Omega} \le C \|w\|_{s,\Omega_i} \quad \forall w \in H^s(\Omega_i), \quad s = 0, 1, 2. \qquad (8)$$

Let $I_h$ be the standard nodal interpolation operator and define

$$I_h^* v := (I_{h,1}^* v_1, I_{h,2}^* v_2) \quad \text{where} \quad I_{h,i}^* v_i := (I_h E_i v_i)|_{\Omega_i}. \qquad (9)$$

The following theorem is valid.

**Theorem 2.** Let $I_h^*$ be an interpolation operator defined as in (9). Then

$$\|v - I_h^* v\| \le C_A h \|v\|_{2,\Omega_1 \cup \Omega_2}, \quad \forall v \in H_0^1(\Omega) \cap H^2(\Omega_1 \cup \Omega_2).$$

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 \in G_h$ onto the unit reference triangle $\tilde{K}$ by an affine map and denote by $\tilde{\Gamma}_{\tilde{K}}$ the corresponding image of $\Gamma_K$. Under assumptions A1–A3 of Section 3 there exist a constant C, depending on $\Gamma$ but independent of the mesh, such that

$$ \|w\|_{0, \tilde{\Gamma}_{\tilde{K}}}^2 \le C \|w\|_{0, \tilde{K}} \|w\|_{1, \tilde{K}}, \quad \forall w \in H^1(\tilde{K}). \quad (10) $$

**Proof.** We start by showing that

$$ \|w\|_{0, \tilde{\Gamma}_{\tilde{K}}}^2 \le C (\|w\|_{0, \tilde{\Gamma}_{\tilde{K},h}}^2 + \|w\|_{0, \tilde{K}} \|w\|_{1, \tilde{K}}). \quad (11) $$

Recall that $\Gamma_{K,h}$ is the straight line connecting the points of intersection between $\Gamma$ and the element $K$ and

$$ \Gamma_K = \{(\xi, \eta) : 0 < \xi < |\Gamma_{K,h}|, \eta = \delta(\xi)\}. $$

Assume first that $\delta(\xi) > 0$. Since the curvature of the interface is bounded, $|\delta'(\xi)| \le C|\Gamma_{K,h}|$. As the mesh is non-degenerate this implies that on the reference element we may write, using again $(\xi, \eta)$ as local coordinates,

$$ \tilde{\Gamma}_{\tilde{K}} = \{(\xi, \eta) : 0 < \xi < |\tilde{\Gamma}_{\tilde{K},h}|, \eta = \tilde{\delta}(\xi)\}, $$

where $|\tilde{\delta}'(\xi)| \le C$. We now let $D$ denote the domain bounded by $\tilde{\Gamma}_{\tilde{K}}$ and $\tilde{\Gamma}_{\tilde{K},h}$ and note that by the divergence theorem,

$$ 2 \int_D w \frac{\partial w}{\partial \eta} d\xi d\eta = \int_D \operatorname{div}(0, w^2) d\xi d\eta = - \int_{\tilde{\Gamma}_{\tilde{K},h}} w^2 d\xi + \int_{\tilde{\Gamma}_{\tilde{K}}} w^2 (1 + (\tilde{\delta}')^2)^{-1/2} ds. \quad (12) $$

As $\tilde{\delta}'$ is bounded,

$$ \|w\|_{0, \tilde{\Gamma}_{\tilde{K}}}^2 \le C \int_{\tilde{\Gamma}_{\tilde{K}}} w^2 (1 + (\tilde{\delta}')^2)^{-1/2} ds, $$

whence (11) follows from (12) using Cauchy-Schwarz' inequality.

In a general case where $\delta$ may switch sign, the same argument may be applied for each part between the intersections of $\tilde{\Gamma}_{\tilde{K}}$ and $\tilde{\Gamma}_{\tilde{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 $\tilde{K}_t$ and the quadrilateral part $\tilde{K}_q$ of $\tilde{K}$ onto new reference domains. We may assume that $\tilde{\Gamma}_{\tilde{K},h}$ intersects $\tilde{K}$ in $(\alpha, 0)$ and in $(0, \beta)$, and, by symmetry, that $0 \le \alpha \le \beta \le 1$.

For $\alpha = \beta = 1$, the desired trace inequality

$$ \|w\|_{0, \tilde{\Gamma}_{\tilde{K},h}}^2 \le C \|w\|_{0, \tilde{K}} \|w\|_{1, \tilde{K}} \quad (13) $$

is valid. For $1/2 < \alpha \le \beta < 1$, we may map the triangular part $\tilde{K}_t$ onto the unit reference triangle by a linear map. By the bound from below on $\alpha$ and $\beta$, this map is bounded, uniformly in $\alpha$ and $\beta$, with uniformly bounded inverse, and hence (13) is valid also in this case. For $1/2 < \alpha < \beta = 1$ the same argument holds, choosing this time $\tilde{K}_t$ as the triangular part which contains the origin.

Assume now that $\alpha \le 1/2$. Let

$$ (\hat{x}, \hat{y}) = M(\tilde{x}, \tilde{y}) = (\tilde{y}, (1-\alpha)^{-1}(\tilde{x} + \tilde{y} - 1)). $$

Then the image $\tilde{K}_q = M(\tilde{K}_q)$ has its corners in $(0, 0)$, $(1, 0)$, $(0, 1)$, $\hat{P} = (\beta, (1-\beta)/(1-\alpha))$, and there holds

$$ \|w\|_{0, \tilde{\Gamma}_{\tilde{K},h}}^2 \le C(\hat{P}) \|w\|_{0, \tilde{K}_q} \|w\|_{1, \tilde{K}_q}. \quad (14) $$
An additional argument is needed to show uniformity in $\hat{P}$. Since $0 \le \alpha \le 1/2$ and $\alpha \le \beta \le 1$, $\hat{P}$ varies in the domain

$$ \hat{D} := \{0 \le \hat{x} \le 1/2, 1 - \hat{x} \le \hat{y} \le 1\} \cup \{1/2 \le \hat{x} \le 1, 1 - \hat{x} \le \hat{y} \le 2(1 - \hat{x})\} $$

as $\alpha$ and $\beta$ vary. Let

$$ F(\hat{P}, \hat{w}) = \frac{\|w\|_{0, \hat{K}_h}^2}{\|w\|_{0, \hat{K}_q} \|w\|_{1, \hat{K}_q}}. $$

We will show that $F(\hat{P}) = \sup_{w \in H^1(\hat{K}_q)} F(\hat{P}, \hat{w})$ is uniformly bounded. For points $\hat{R}$ and $\hat{S}$ in $\hat{D}$, assuming without restriction that $F(\hat{R}) \ge F(\hat{S})$, we have for any $w$ that

$$ F(\hat{R}) - F(\hat{S}) = \sup_{\hat{v}} F(\hat{R}, \hat{v}) - \sup_{\hat{v}} F(\hat{S}, \hat{v}) \le |\sup_{\hat{v}} F(\hat{R}, \hat{v}) - F(\hat{R}, \hat{w})| + |F(\hat{R}, \hat{w}) - F(\hat{S}, \hat{w})| = I + II. $$

Given $\epsilon > 0$ we may choose $\hat{w}$ such that $I \le \epsilon/2$. Note that $F(\hat{R}, \hat{v})$ is continuous for fixed $\hat{v}$ since the only dependence of $\hat{R}$ lies in the domains of integration. We may thus take $|\hat{R} - \hat{S}|$ small enough so that $II \le \epsilon/2$. Hence $F(\hat{P})$ is continuous on the compact set $\hat{D}$, and thus (14) holds uniformly in $\hat{P}$. Finally, since $M$ is bounded, uniformly in $\alpha$ and $\beta$, with uniformly bounded inverse, (13) follows and the proof is complete. $\square$

**Proof of Theorem 2.** Recall that $K_i = K \cap \Omega_i$ and let $v_i^* = E_i v_i$ denote the extension of $v_i$ to $\Omega$. By a standard interpolation estimate we obtain

$$ \|\nabla(v_i - I_{h,i}^* v_i)\|_{0, \hat{K}_i} = \|\nabla(v_i^* - I_h v_i^*)\|_{0, \hat{K}_i} \le \|\nabla(v_i^* - I_h v_i^*)\|_{0, \hat{K}} \le Ch \|v_i^*\|_{2, \hat{K}}. $$

Summing over all triangles that intersect $\Omega_i$, it follows by (8) that

$$ \|\nabla(v_i - I_{h,i}^* v_i)\|_{0, \Omega_i}^2 \le Ch^2 \sum_{K \cap \Omega_i \ne \emptyset} \|v_i^*\|_{2, \hat{K}}^2 \le Ch^2 \|v_i\|_{2, \Omega_i}^2. \quad (15) $$

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_K^{-1} \|w\|_{0, \Gamma_K}^2 \le C(h_K^{-2} \|w\|_{0, \Gamma_K}^2 + \|w\|_{1, \Gamma_K}^2), \quad \forall w \in H^1(K). $$

Hence it follows, using again a standard interpolation estimate, that

$$ \begin{aligned} h_K^{-1} \|[v - I_h^* v]\|_{0, \Gamma_K}^2 &\le Ch_K^{-1} \sum_i \|v_i - I_{h,i}^* v_i\|_{0, \Gamma_K}^2 = Ch_K^{-1} \sum_i \|v_i^* - I_h v_i^*\|_{0, \Gamma_K}^2 \\ &\le C \sum_i (h_K^{-2} \|v_i^* - I_h v_i^*\|_{0, \Gamma_K}^2 + \|v_i^* - I_h v_i^*\|_{1, \Gamma_K}^2) \le Ch_K^2 \sum_i \|v_i^*\|_{0, \Gamma_K}^2. \end{aligned} $$

Summing the contributions from $K \in G_h$, we get from (8) that

$$ \|[v - I_h^* v]\|_{1/2, h, \Gamma} \le Ch \sum_{i=1}^2 \|v_i^*\|_{2, \cup K \in G_h} \le Ch \|v\|_{2, \Omega_1 \cup \Omega_2}. \quad (16) $$

Finally, Lemma 3 applied to $\nabla_n w$ and scaling gives

$$ h_K \|\nabla_n w\|_{0, \Gamma_K}^2 \le C(\|w\|_{1, \Gamma_K}^2 + h_K^2 \|w\|_{2, \Gamma_K}^2), \quad \forall w \in H^2(K), $$

whence similar arguments as above yield

$$ \|\nabla_n (v_i - I_{h,i}^* v_i)\|_{-1/2, h, \Gamma} \le Ch \|v_i\|_{2, \Omega_i}. \quad (17) $$

Since $\kappa_i < 1$, the theorem now follows from (15)–(17). $\square$
## 5. A priori error estimates

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

**Lemma 4.** For $\phi \in V_h$, the following inverse inequality holds:

$$ \|\{\nabla_n \phi\}\|_{-1/2, h, \Gamma}^2 \le C_I \|\nabla \phi\|_{0, \Omega_1 \cup \Omega_2}^2. $$

**Proof.** Since $\phi \in V_h$ is linear on $K_i$, we have

$$ h_K \|\kappa_i \nabla_n \phi_i\|_{0, \Gamma_K}^2 \le h_K \kappa_i^2 |\Gamma_K| \|\nabla \phi_i\|_{0, K_i}^2 = h_K \frac{|\Gamma_K|}{|K_i|^2} \|\nabla \phi_i\|_{0, K_i}^2 \le C \|\nabla \phi_i\|_{0, K_i}^2. $$

In the last step above we have used that $|\Gamma_K| \le h_K$, $|K_i| \le h_K^2$, and, since the mesh is non-degenerate, $|K| \ge ch_K^2$. The result follows by summation over the elements. $\square$

**Lemma 5.** The discrete form $a_h(\cdot, \cdot)$ is coercive on $V^h$, i.e.,

$$ a_h(v, v) \ge C |||v|||^2 \quad \forall v \in V^h, $$

provided $\lambda$ is chosen sufficiently large. It is also continuous, i.e.,

$$ a_h(u, v) \le C |||u||| \quad |||v||| \quad \forall u \in V, \forall v \in V. $$

**Proof.** Continuity of the discrete form follows directly from the definitions. To prove coercivity, we use (7) to find that for any $\epsilon > 0$

$$ \begin{aligned} a_h(v, v) &= \|\alpha^{1/2} \nabla v\|_{0, \Omega_1 \cup \Omega_2}^2 - 2([v], \{\alpha \nabla_n v\})_\Gamma + \|\lambda^{1/2} [v]\|_{0, \Gamma}^2 \\ &\ge \|\alpha^{1/2} \nabla v\|_{0, \Omega_1 \cup \Omega_2}^2 - 2\|\{\alpha \nabla_n v\}\|_{-1/2, h, \Gamma} \|[v]\|_{1/2, h, \Gamma} + \|\lambda^{1/2} [v]\|_{0, \Gamma}^2 \\ &\ge \|\alpha^{1/2} \nabla v\|_{0, \Omega_1 \cup \Omega_2}^2 - \frac{1}{\epsilon} \|\{\alpha \nabla_n v\}\|_{-1/2, h, \Gamma}^2 + \sum_{K \in G_h} \left(\lambda - \frac{\epsilon}{h_K}\right) \|[v]\|_{0, \Gamma_K}^2. \end{aligned} $$

It then follows from Lemma 4 that

$$ \begin{aligned} a_h(v, v) &\ge \frac{1}{2} \|\alpha^{1/2} \nabla v\|_{0, \Omega_1 \cup \Omega_2}^2 + \left( \frac{1}{2} - \frac{2C_I \max_{\Omega} \alpha}{\epsilon} \right) \|\alpha^{1/2} \nabla v\|_{0, \Omega_1 \cup \Omega_2}^2 \\ &\quad + \frac{1}{\epsilon} \|\{\alpha \nabla_n v\}\|_{-1/2, h, \Gamma}^2 + \sum_{K \in G_h} \left( \lambda - \frac{\epsilon}{h_K} \right) \|[v]\|_{0, \Gamma_K}^2. \end{aligned} $$

Taking $\epsilon = 4C_I \max_{\Omega} \alpha$, coercivity follows if $\lambda|_K = \gamma h_K^{-1}$ where $\gamma > 4C_I \max_{\Omega} \alpha$. $\square$

**Theorem 6.** Under assumptions A1–A3 of Section 3, and for U solving (4) and u solving (1), the following a priori error estimates hold:

$$ |||u - U||| \le Ch \|u\|_{2, \Omega_1 \cup \Omega_2} \quad (18) $$

and

$$ \|u - U\|_{0, \Omega} \le Ch^2 \|u\|_{2, \Omega_1 \cup \Omega_2}. \quad (19) $$
**Proof.** For any $v \in V^h$, $|||u - U||| \le |||u - v||| + |||v - U|||$. Further, by Lemma 5 and orthogonality, we have that

$$ |||U - v|||^2 \le C a_h(U - v, U - v) = C a_h(u - v, U - v) \le C |||u - v||| |||U - v|||, $$

and it follows that

$$ |||u - U||| \le C |||u - v||| \quad \forall v \in V^h. $$

Taking $v = I_h^*u$ and invoking the interpolation result of Theorem 2, (18) follows.

For (19) we use a duality argument. Define $z = (z_1, z_2)$ by

$$ \begin{align} - \nabla \cdot (\alpha_i \nabla z_i) &= e_i \quad \text{in } \Omega_i, \quad i = 1, 2, \tag{20} \\ z_i &= 0 \quad \text{on } \partial\Omega \cap \partial\Omega_i, \notag \\ [z] &= 0 \quad \text{on } \Gamma, \notag \\ [\alpha \nabla_n z] &= 0 \quad \text{on } \Gamma, \notag \end{align} $$

where $e_i = u_i - U_i$. By Green's formula and (5) with $u = z, g = 0$, we have that

$$ \begin{align*} \|e\|_{0,\Omega}^2 &= -(\nabla \cdot (\alpha \nabla z), e)_{\Omega_1 \cup \Omega_2} = (\alpha \nabla z, \nabla e)_\Omega - (\alpha_1 \nabla_n z_1, e_1) + (\alpha_2 \nabla_n z_2, e_2) \\ &= (\alpha \nabla z, \nabla e)_\Omega - (\{\alpha \nabla_n z\}, [e])_\Gamma = a_h(z, e) \end{align*} $$

since $[z] = 0$. Thus, using the symmetry of $a_h(\cdot, \cdot)$ and applying the orthogonality relation (6) and Theorem 2, we find that

$$ \|e\|_{0,\Omega}^2 = a_h(z - I_h z, e) \le C |||z - I_h z||| \quad |||e||| \le C h \|z\|_{2,\Omega_1 \cup \Omega_2} \quad |||e|||. \quad (21) $$

Finally, by the elliptic regularity result (3), we have $\|z\|_{2,\Omega_1 \cup \Omega_2} \le C\|e\|_{0,\Omega}$, whence the estimate (19) follows from (21) and (18). $\square$

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

$$ E_K(U) = (h_K^2 \|f + \nabla \cdot (\alpha \nabla U)\|_{0,K_1 \cup K_2}^2 + h_K \|[\alpha_i \nabla U_i]\|_{0,\partial K}^2 + h_K \|g - [\alpha \nabla U]\|_{0,\Gamma_K}^2 + h_K^{-1} \|[U]\|_{0,\Gamma_K}^2)^{1/2} \quad (22) $$

and

$$ E(U) = \left( \sum_{K \in \mathcal{T}_h} h_K^2 E_K(U)^2 \right)^{1/2} . \quad (23) $$

We then have the following a posteriori error estimate.

**Theorem 7.** For a continuous linear functional $j(\cdot)$ on $L_2(\Omega)$, let $J \in L_2(\Omega)$ be defined by Riesz' representation theorem, i.e., $j(\cdot) := (J, \cdot)_\Omega$. Then there is a positive constant $C$ such that

$$ j(e) \le CE(U) \|J\|_{0,\Omega}. \quad (24) $$
**Proof.** Let $z$ be the solution to the problem

$$ 
\begin{array}{l}
-\nabla \cdot (\alpha_i \nabla z_i) = J \quad \text{in } \Omega_i, \quad i = 1, 2, \\
z_i = 0 \quad \text{on } \partial\Omega \cap \partial\Omega_i, \\
[z] = 0 \quad \text{on } \Gamma, \\
[\alpha \nabla_n z] = 0 \quad \text{on } \Gamma.
\end{array}
\qquad (25) $$

We first note that

$$ (e, J)_{\Omega} = (\alpha \nabla e, \nabla z)_{\Omega_1 \cup \Omega_2} - (\alpha_1 \nabla_n z_1, e_1)_{\Gamma} + (\alpha_2 \nabla_n z_2, e_2)_{\Gamma}. $$

Now, since $[u] = 0$, we see that, by (5),

$$ -(\alpha_1 \nabla_n z_1, e_1)_\Gamma + (\alpha_2 \nabla_n z_2, e_1)_\Gamma = (U_1, \alpha_1 \nabla_n z_1)_\Gamma - (U_2, \alpha_2 \nabla_n z_2)_\Gamma = ([U], \{\alpha \nabla_n z\}). $$

Thus,

$$ j(e) = (e, J)_{\Omega} = (\alpha \nabla e, \nabla z)_{\Omega_1 \cup \Omega_2} + ([U], \{\alpha \nabla_n z\}). \qquad (26) $$

Now, take $Z \in V_h$. From Galerkin orthogonality we then have $a_h(e, Z) = 0$, and since $[z] = 0$, $[u] = 0$ on the interface, we get

$$ 
\begin{aligned}
0 &= a_h(e, Z) = (\alpha \nabla e, \nabla Z)_{\Omega_1 \cup \Omega_2} - ([e], \{\alpha \nabla_n Z\})_\Gamma - (\{\alpha \nabla_n e\}, [Z])_\Gamma + (\lambda[e], [Z])_\Gamma \\
&= (\alpha \nabla e, \nabla Z)_{\Omega_1 \cup \Omega_2} + ([U], \{\alpha \nabla_n Z\})_\Gamma + (\{\alpha \nabla_n e\}, [z-Z])_\Gamma + (\lambda[U], [z-Z])_\Gamma.
\end{aligned}
\qquad (27) $$

Denote by $n_K$ the outward pointing unit normal to $K$. Subtracting (27) from (26) and integrating by parts we get

$$ 
\begin{aligned}
j(e) &= (\alpha \nabla e, \nabla (z-Z))_{\Omega_1 \cup \Omega_2} + ([U], \{\alpha \nabla_n (z-Z)\})_\Gamma - (\{\alpha \nabla_n e\}, [z-Z])_\Gamma - (\lambda[U], [z-Z])_\Gamma \\
&= \sum_{K \in T_h} (f - \nabla \cdot (\alpha \nabla U), z-Z)_{K_1 \cup K_2} + ([U], \{\alpha \nabla_n (z-Z)\})_\Gamma - (\{\alpha \nabla_n e\}, [z-Z])_\Gamma - (\lambda[U], [z-Z])_\Gamma \\
&\quad - \frac{1}{2} \sum_{K \in T_h} ([\alpha n_K \cdot \nabla U], z-Z)_{\partial K \setminus \Gamma} + (\alpha_1 \nabla_n e_1, z_1 - Z_1)_\Gamma - (\alpha_2 \nabla_n e_2, z_2 - Z_2)_\Gamma.
\end{aligned}
 $$

We now note that

$$ 
\begin{aligned}
& (\alpha_1 n \cdot \nabla e_1, z_1 - Z_1)_\Gamma - (\alpha_2 n \cdot \nabla e_2, z_2 - Z_2)_\Gamma - (\{\alpha \nabla_n e\}, [z-Z])_\Gamma \\
&= (\alpha_1 \nabla_n e_1, z_1 - Z_1)_\Gamma - (\alpha_2 \nabla_n e_2, z_2 - Z_2)_\Gamma - (\kappa_1 \alpha_1 \nabla_n e_1 + \kappa_2 \alpha_2 \nabla_n e_2, [z-Z])_\Gamma \\
&= (\kappa_2 [\alpha \nabla_n e], z_1 - Z_1)_\Gamma + (\kappa_1 [\alpha \nabla_n e], z_2 - Z_2)_\Gamma \\
&= (\kappa_2 (g - [\alpha \nabla_n U]), z_1 - Z_1)_\Gamma + (\kappa_1 (g - [\alpha \nabla_n U]), z_2 - Z_2)_\Gamma \\
&= \sum_{j=1}^2 ((1 - \kappa_j)(g - [\alpha \nabla_n U]), z_j - Z_j)_\Gamma
\end{aligned}
 $$

and thus we find that

$$ 
\begin{aligned}
j(e) = & \sum_{K \in T_h} \left( (f - \nabla \cdot (\alpha \nabla U), z-Z)_{K_1 \cup K_2} - \frac{1}{2} ([\alpha n_K \cdot \nabla U], z-Z)_{\partial K \setminus \Gamma} \right) \\
& + \sum_{K \in G_h} ([U], \{\alpha \nabla_n (z-Z)\} - \lambda[z-Z])_{\Gamma_K} \\
& + \sum_{K \in G_h} \sum_{j=1}^2 ((1 - \kappa_j)(g - [\alpha \nabla_n U]), z_j - Z_j)_{\Gamma_K}.
\end{aligned}
\qquad (28) $$
Further, by Cauchy–Schwarz' inequality, and choosing $\lambda|_K = \gamma h_K^{-1}$ with $\gamma > 4C_T \max_{\Omega} \alpha_i$,

$$j(e) \le \sum_{K \in T_h} \sum_{j=1}^{2} \rho_{K,j} \omega_{K,j} + \sum_{K \in G_h} \sum_{j=1}^{2} \rho_{S,j} \omega_{S,j}, \qquad (29)$$

where

$$\rho_{K,1} = h_K \|f + \nabla \cdot (\alpha \nabla U)\|_{0,K_1 \cup K_2}, \quad \omega_{K,1} = h_K^{-1} \|z - Z\|_{0,K_1 \cup K_2},$$

$$\rho_{K,2} = \frac{1}{2} h_K^{1/2} \|[\alpha \mathbf{n}_K \cdot \nabla U]\|_{0, \partial K \setminus \Gamma}, \quad \omega_{K,2} = h_K^{-1/2} \|z - Z\|_{0, \partial K \setminus \Gamma},$$

$$\rho_{S,1} = h_K^{-1/2} \|[u]\|_{0,\Gamma_K}, \quad \omega_{S,1} = h_K^{1/2} \|\{\alpha \nabla_n (z-Z)\} - \gamma[z-Z]/h_K\|_{0,\Gamma_K}$$

and

$$\rho_{S,2} = h_K^{1/2} \|g - [\alpha \nabla_n U]\|_{0,\Gamma_K}, \quad \omega_{S,2} = h_K^{-1/2} \max_j \|z_j - Z_j\|_{0,\Gamma_K}.$$

Note that the right hand side of (29) is a weighted sum of local residuals $\rho_{K,i}$, $\rho_{S,i}$. To conclude the proof we bound the weights in (29) by choosing $Z = I_h^* z$, applying the interpolation estimates (15)–(17) and finally the regularity estimate

$$\|z\|_{2,\Omega_1 \cup \Omega_2} \le C \|J\|_{0,\Omega},$$

whence (24) follows. $\square$

## 7. 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 $\Omega_1$ and $\Omega_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 $\Gamma$. For each $K \in G_h$, we assigned two identical copies $K'$ and $K''$. We assumed that $K$, with nodes $\{i, j, k\}$ was split by $\Gamma$ 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'$ the triangular part and to $K''$ the quadrilateral part. Thus, on $K''$ we created a new node $i'$ on the far side of $\Gamma$, and on $K'$ we created two new nodes, $k'$ and $j'$, on the other side of $\Gamma$, see Fig. 1. To ensure continuity across the edges of $K'$ and $K''$ (away from $\Gamma$), 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 $\Omega_1$, and the other completely covering $\Omega_2$. The elements crossed by $\Gamma$ had been doubled, but coincided geometrically.

To numerically evaluate $a_h(\cdot, \cdot)$ and $L(\cdot)$, we used the following strategy. The triangles that were not crossed by $\Gamma$ were handled in the usual way. On the elements crossed by $\Gamma$, we used centroid quadrature to
<figure><img src="image_2.png" /><figcaption>Fig. 1. The split of a triangle used in the implementation of the proposed method.</figcaption></figure>

evaluate all terms, both on the quadrilateral side and on the triangular side. On the interface, we used two-point 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', j', and k' 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.

## 7.1. Example 1

We considered solutions to the ordinary differential equation

$$ -\sum_{i} \frac{d}{dx} \left( \alpha_{i} \frac{du_{i}}{dx} \right) = 1; \quad [u(1/2)] = 0; \quad \alpha_{1} \frac{du_{1}}{dx}(1/2) = \alpha_{2} \frac{du_{2}}{dx}(1/2). $$

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) \times (0, 1)$, with zero Neumann boundary conditions at $y = 0$ and $y = 1$. The equation has a closed-form solution, given by

$$ u_1(x) = \frac{(3\alpha_1 + \alpha_2)x}{4\alpha_1^2 + 4\alpha_1\alpha_2} - \frac{x^2}{2\alpha_1}, \quad u_2(x) = \frac{\alpha_2 - \alpha_1 + (3\alpha_1 + \alpha_2)x}{4\alpha_2^2 + 4\alpha_1\alpha_2} - \frac{x^2}{2\alpha_2}. $$

We chose $\alpha_1 = 1/2$, $\alpha_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 elevations of the solutions obtained with the different unfitted methods are shown on the final mesh consisting of 2048 elements.
<figure><img src="image_3.png" /><figcaption>Fig. 2. L₂-norm convergence of different methods applied to the interface problem.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Elevation of the discontinuous approximation obtained from the proposed method (left) and the continuous standard FE approximation (right) on the final mesh.</figcaption></figure>

## 7.2. Example 2

Here we considered a less trivial two-dimensional example (from [7]). The exact solution is given by

$$u(x,y) = \begin{cases} \frac{r^2}{\alpha_1}, & \text{if } r \le r_0, \\ \frac{r^2}{\alpha_2} - \frac{r_0^2}{\alpha_2} + \frac{r_0^2}{\alpha_1}, & \text{if } r > r_0, \end{cases}$$

where $r := \sqrt{x^2 + y^2}$, and, on the domain $(0, 1) \times (0, 1)$, we chose $r_0 = 3/4$, $\alpha_1 = 1$, and $\alpha_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₂- 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₂-norm convergence, though we have no theoretical results to back this up.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Elevation of the approximate solution.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. $L_2$-norm and maximum norm convergence.</figcaption></figure>

### 7.3. Example 3

The third example was solved on the domain $(0, 1) \times (0, 1)$, with zero Dirichlet boundary conditions. Centred in the domain is an ellipsoidal inclusion with conduction parameter $\alpha_{\min} = 1$; outside of the ellipse we set $\alpha_{\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(\Omega)$-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.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. First mesh and elevation of the approximate solution.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Last adaptively refined mesh and elevation of the approximate solution.</figcaption></figure>

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

## References

[1] I. Babuška, The finite element method for elliptic equations with discontinuous coefficients, Computing 5 (1970) 207–213.

[2] J.W. Barrett, C.M. Elliott, Fitted and unfitted finite-element methods for elliptic equations with smooth interfaces, IMA J. Numer. Anal. 7 (1987) 283–300.

[3] R. Becker, R. Rannacher, A feed-back approach to error control in finite element methods: basic analysis and examples, East-West J. Numer. Math. 4 (1996) 237–264.

[4] J.H. Bramble, J.T. King, A finite element method for interface problems in domains with smooth boundaries and interfaces, Adv. Comput. Math. 6 (1996) 109–138.

[5] Z. Chen, J. Zhou, Finite element methods and their convergence for elliptic and parabolic interface problems, Numer. Math. 79 (1998) 175–202.

[6] M. Feistauer, V. Sobotíková, Finite element approximation of nonlinear problems with discontinuous coefficients, M2AN Math. Model. Numer. Anal. 24 (1990) 457–500.
[7] Z. Li, T. Lin, X. Wu, New Cartesian grid methods for interface problems using the finite element formulation, NCSU Technical Report: CRSC-TR99-12, Center for Research in Scientific Computation, North Carolina State University, 1999.

[8] R.J. MacKinnon, G.F. Carey, Treatment of material discontinuities in finite element computations, Int. J. Numer. Meth. Eng. 24 (1987) 393–417.

[9] J. Nitsche, Über ein Variationsprinzip zur Lösung von Dirichlet–Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind, Abh. Math. Univ. Hamburg 36 (1971) 9–15.

[10] A. Ženíšek, The finite element method for nonlinear elliptic equations with discontinuous coefficients, Numer. Math. 58 (1990) 51–77.