# A hierarchical NXFEM for fictitious domain simulations

Roland Becker, Erik Burman, Peter Hansbo

## ► To cite this version:

Roland Becker, Erik Burman, Peter Hansbo. A hierarchical NXFEM for fictitious domain simulations. International Journal for Numerical Methods in Engineering, 2011, 86 (4-5), pp.549-559.
10.1002/nme.3093 . hal-01352895

HAL Id: hal-01352895

https://hal.science/hal-01352895

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
A hierarchical NXFEM for fictitious domain simulations

Roland Becker¹, Erik Burman² and Peter Hansbo³

¹Laboratoire de Mathématiques Appliques, Université de Pau et des Pays de l'Adour, BP 1155,
64013 PAU Cedex, France

²Department of Mathematics, University of Sussex Falmer, BN1 9RF UK, U.K.

³Department of Mathematical Sciences, Chalmers University of Technology and University of Gothenburg,
SE-41296 Göteborg, Sweden

SUMMARY

We suggest a fictitious domain method, based on the Nitsche XFEM method of (Comput. Meth. Appl. Mech. Engrg 2002; **191**:5537–5552), that employs a band of elements adjacent to the boundary. In contrast, the classical fictitious domain method uses Lagrange multipliers on a line (surface) where the boundary condition is to be enforced. The idea can be seen as an extension of the Chimera method of (ESAIM: Math. Model Numer. Anal. 2003; **37**:495–514), but with a hierarchical representation of the discontinuous solution field. The hierarchical formulation is better suited for moving fictitious boundaries since the stiffness matrix of the underlying structured mesh can be retained during the computations.

Our technique allows for optimal convergence properties irrespective of the order of the underlying finite element method.

KEY WORDS: Nitsche's method; fictitious domain; extended finite element method

# 1. INTRODUCTION

Fictitious domain methods were introduced in order to be able to use Cartesian meshes also for solving problems on domains with complex boundaries. The idea is to enforce Dirichlet boundary conditions on a given curve (surface) that is discretized independently of the mesh, cf. Glowinski *et al.* [1]. On this curve, the boundary condition is enforced, typically by the use of Lagrange multipliers. The system of equations can then be set up on a Cartesian mesh and the degrees of freedom falling outside of the boundary are discarded. The problem with this approach is that the derivatives of the finite element solution normal to the curve cannot accommodate the jump necessary to achieve optimal order convergence, cf. [2]. Another problem is how to choose the discretization of the curve relative to the elements it crosses in order for the problem to be well posed. Guidelines for this are given in [2] but they are by necessity rather vague.

In this paper we introduce an alternative method based on the use of Nitsche's method in the vein of Hansbo *et al.* [3] (building on [4, 5]), where overlapping meshes were considered. We shall also employ overlapping meshes in the form of (Figure 1)

1. the (structured) mesh on which the problem is set up and

2. a narrow band of elements that overlays the first mesh.
<figure><img src="image_2.png" /><figcaption>Figure 1. Underlying and overlying meshes.</figcaption></figure>

This allows for the direct use of the method proposed in [3], where the elements on the underlying mesh were cut by the overlying and the solution pasted together by the use of Nitsche's method [6]. The outer boundary of the band can then be used as the Dirichlet boundary. We remark at this point that another strength of the approach in [3] is that any boundary condition can be applied at the outer boundary of the band. This is not so straightforward in a classical fictitious domain method.

However, in [3] this was achieved by modifying the elements of the underlying mesh, which does not allow for the system matrix on the underlying mesh to remain unchanged. If the boundary was to move continuously, the entries in the system matrix would also have to be changed continuously. It is desirable to have a fixed matrix for the underlying problem and to see the imposition of the boundary condition as a set of additional degrees of freedom (as in the original fictitious domain approach where the additional degrees of freedom are the Lagrange multipliers).

As has been noted by Areias and Belytschko [7], the method of [3] can be interpreted as an extended finite element method (XFEM) method (which instead adds degrees of freedom in a hierarchical fashion) by a reordering of the degrees of freedom. In this paper, we seize on the hierarchical concept to construct a fictitious domain method which is formally identical to that of [3] (thus benefitting from the optimal order error analysis therein) while still keeping the underlying system matrix unchanged by the location of the boundary.

Our method shares some characteristics with the fat boundary method proposed by Maury [8], which also makes use of an auxiliary band-like domain. However, in [8], the original Poisson problem is recast on the continuous level and requires a fixed point iteration scheme to solve, like in classical domain decomposition methods. In contrast, our method is defined on the discrete level and can be solved monolithically.

## 2. MODEL PROBLEM AND FINITE ELEMENT SPACES

### 2.1. The continuous problem

As a model problem, we consider the Poisson equation

$$ -\Delta u = f \quad \text{in } \Omega, \quad u=0 \quad \text{on } \partial\Omega, \qquad (1) $$

where $\Omega$ is a domain in $\mathbb{R}^2$ with polygonal boundary $\partial\Omega$ and $f$ is a given forcing term. We embed $\Omega$ in a larger rectangular domain $\hat{\Omega}$ so that $\Omega$ is completely contained in the interior on $\hat{\Omega}$. Finally, we introduce a third domain $\Omega_1$ consisting of a band (the width of which may be mesh dependent) whose outer boundary coincides with $\partial\Omega$ and whose inner boundary forms a line, denoted $\Gamma$, in the
<figure><img src="image_3.png" /><figcaption>Figure 2. Boundary and domain definitions.</figcaption></figure>

interior of $\Omega$. The remainder of $\Omega$ is denoted by $\Omega_2 := \Omega \setminus \Omega_1$, cf. Figure 2). Clearly, the extension to three dimensions is straightforward.

We now rephrase the problem (1) as an interface problem. For $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_i = u|_{\Omega_i}$ is the restriction of $u$ to $\Omega_i$. 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$. Then we may formulate (1) as:

$$ 
\begin{align} 
-\Delta u &= f \quad \text{in } \Omega_1 \cup \Omega_2, \notag \\ 
u &= 0 \quad \text{on } \partial\Omega, \notag \\ 
[u] &= 0 \quad \text{on } \Gamma, \notag \\ 
[\partial_n u] &= 0 \quad \text{on } \Gamma. \tag{2}
\end{align} 
$$

Here $n$ is the outward pointing unit normal to $\Omega_1$ and $\partial_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}$. The inner products in $H^0(D) = L_2(D)$ are denoted as $(\cdot, \cdot)_D$.

## 2.2. Finite element spaces

Denote by $T_1^h$ the triangulation of $\Omega_1$ and by $T_2^h$ the triangulation of $\hat{\Omega}$. We shall make a discretization on the whole of $\hat{\Omega}$ even though the solution has no physical significance outside $\partial\Omega$. This is in line with classical fictitious domain methods and means that the stiffness matrix assembled from $T_2^h$ remains fixed even if the domain should change, as it must do in many applications.

We will use the following notation for mesh-related quantities. Let $h_K$ be the diameter of an element $K \in T_1^h$ and $h = \max_{K \in T_1^h, i=1,2} h_K$. To distinguish elements from the two meshes, we will sometimes use indexed element notation $K_i \in T_1^h$ for clarity. Furthermore, we introduce $G^h$ as the set of elements in $T_2^h$ intersected by $\Gamma$,

$$ G^h = \{K \in T_2^h : K \cap \Gamma \neq \emptyset\}, $$

and the corresponding mesh-dependent domain

$$ \Omega_{G_h} = \bigcup_{G_h} K. $$

We shall also need the mesh-dependent boundary $\partial\Omega_{G_h}$, which consists of the edges on elements in $G^h$ that form the boundary of $\Omega_{G_h}$. This boundary is also split in the part exterior to $\Gamma$, $\partial\Omega_{G_h}^{\text{ext}}$ and interior to $\Gamma$, $\partial\Omega_{G_h}^{\text{int}}$.
The nodes on Γ of the elements in $T_1^h$, together with the points of intersection between elements in $T_2^h$ and Γ, define a partition of Γ, $\Gamma = \bigcup_{j \in J_h} \bar{\Gamma}^j$. Note that each part $\Gamma^j$ belongs to two elements, one from each mesh. We denote these elements by $K_1^j$ and $K_2^j$, respectively. A local meshsize on Γ is defined by:

$$h(x) = h_{K_1^j}, \quad x \in \Gamma^j. \qquad (3)$$

For any element $K \in T_i^h$, let $P_K = K \cap \Omega_i$ denote the part of K in $\Omega_i$.

We make the following assumptions regarding the meshes:

A1: The triangulations are non-degenerate, i.e.

$$h_K / \rho_K \le C \quad \forall K \in T_i^h, \quad i=1,2,$$

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

A2: The meshes have locally compatible meshsize over Γ. More precisely, let $K_1^j \in T_1^h$ and $K_2^j \in T_2^h$ be the elements which contain a specific part $\Gamma^j$ of Γ. We assume that

$$ch_{K_1^j} \le h_{K_2^j} \le Ch_{K_1^j} \quad \forall j \in J_h.$$

Here and below, C and c denote generic constants.

We shall seek a discrete solution $U = (U_1, U_2, U_3, U_4)$ in the space $V^h = V_1^h \times V_2^h \times V_3^h \times V_4^h$, where

$$V_1^h = \{\phi \in H^1(\Omega_1) : \phi|_K \text{ is a polynomial of degree } p \quad \forall K \in T_1^h, \phi|_{\partial\Omega} = 0\},$$

$$V_2^h = \{\phi \in H^1(\hat{\Omega}) : \phi|_K \text{ is a polynomial of degree } p \quad \forall K \in T_2^h, \phi|_{\partial\hat{\Omega}} = 0\},$$

$$V_3^h = \{\phi \in H^1(\Omega_{G_h}) : \phi|_K \text{ is a polynomial of degree } p \quad \forall K \in G_h, \phi|_{\partial\Omega_{G_h}^{\text{int}}} = 0\},$$

$$V_4^h = \{\phi \in H^1(\Omega_{G_h}) : \phi|_K \text{ is a polynomial of degree } p \quad \forall K \in G_h, \phi|_{\partial\Omega_{G_h}^{\text{ext}}} = 0\}.$$

Note that functions in $V^h$ are, in general, discontinuous across Γ; the discontinuity is represented by the hierarchical space $V_3^h \times V_4^h$. In Figure 3 we illustrate the concept: a discontinuous function on a one-dimensional element occupying the domain (0, 1) (solid line) can be written as the sum of a continuous function (dashed line) from $V_h^1$ and piecewise continuous functions which are zero in the nodes of the element (dash-dotted line) from $V_3^h$ and $V_4^h$.

<figure><img src="image_4.png" /><figcaption>Figure 3. A discontinuous trial function and its split into a continuous and a discontinuous part.</figcaption></figure>
### 2.3. The finite element method (FEM)

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

$$a^h(U, \phi) = l(\phi) \quad \forall \phi \in V^h, \quad (4)$$

where, if we denote $\Omega_{G_h}^{\text{int}} := \Omega_{G_h} \cap \Omega_2$ and $\Omega_{G_h}^{\text{ext}} := \Omega_{G_h} \setminus \Omega_{G_h}^{\text{int}}$,

$$\begin{aligned} a^h(U, \phi) = & (\nabla U_1, \nabla \phi_1)_{\Omega_1} + (\nabla U_2, \nabla \phi_2)_{\hat{\Omega}} + (\nabla U_2, \nabla \phi_3)_{\Omega_{G_h}^{\text{int}}} + (\nabla U_2, \nabla \phi_4)_{\Omega_{G_h}^{\text{ext}}} \\ & + (\nabla U_3, \nabla \phi_3)_{\Omega_{G_h}^{\text{int}}} + (\nabla U_4, \nabla \phi_4)_{\Omega_{G_h}^{\text{ext}}} + (\nabla U_3, \nabla \phi_2)_{\Omega_{G_h}^{\text{int}}} + (\nabla U_4, \nabla \phi_2)_{\Omega_{G_h}^{\text{ext}}} \\ & - (\partial_n U_1, [\phi])_{\Gamma} - ([U], \partial_n \phi_1)_{\Gamma} + (\lambda h^{-1}[U], [\phi])_{\Gamma}, \\ l(\phi) = & (f, \phi_1)_{\Omega_1} + (f, \phi_2)_{\hat{\Omega}} + (f, \phi_3)_{\Omega_{G_h}^{\text{int}}} + (f, \phi_4)_{\Omega_{G_h}^{\text{ext}}}, \end{aligned}$$

with $f$ extended, e.g. by zero outside $\Omega$, and where $h$ is the local meshsize (3). Here, the jump $[U]$ is interpreted as $U_1 - (U_2 + U_3)$. The continuity conditions of $u$ and $\partial_n u$ at $\Gamma$ are satisfied weakly by means of a variant of Nitsche's method [6] for consistent weak enforcement of Dirichlet boundary conditions. To ensure stability, the parameter $\lambda$ has to be taken sufficiently large, cf. [3].

In order to show consistency, we first have to define in what sense the method is consistent since the bilinear form $a^h(\cdot, \cdot)$ is defined also outside of $\Omega$. Thus, we introduce a second bilinear form

$$a^h_*(U, \phi) = (\nabla U, \nabla \phi)_{\Omega_1 \cup \Omega_2} - (\partial_n U_1, [\phi])_{\Gamma} - ([U], \partial_n \phi_1)_{\Gamma} + (\lambda h^{-1}[U], [\phi])_{\Gamma},$$

and right-hand side

$$l_*(\phi) = (f, \phi)_{\Omega_1 \cup \Omega_2}.$$

It is straightforward to show that the method analyzed in [3] would then read: find $U^* \in V^h$ such that

$$a^h_*(U^*, \phi) = l_*(\phi) \quad \forall \phi \in V^h, \quad (5)$$

where $U^*|_{\Omega_1} = U_1$, $U^*|_{\Omega_2} = U_2 + U_3$, $\phi|_{\Omega_1} = \phi_1$ and $\phi|_{\Omega_2} = \phi_2 + \phi_3$, since we may ignore the solution outside $\Omega$ we set $\phi_2|_{\hat{\Omega} \setminus (\Omega_2 \cup \Omega_{G_h})} = 0$. Note that the integrals here are taken over only the domain $\Omega$, unlike the proposed method. However, since the solution is completely decoupled by the cut, it does not matter what we do outside of the domain (e.g. how we extend $f$ and how the boundary conditions on $\partial \hat{\Omega}$ are specified). Explaining the method of Areias and Belytschko [7] without the use of function spaces could be done as follows: only the standard FEM is used on the mesh $T_1^h$ on the band, and on the part of the mesh $T_2^h$ on $\hat{\Omega}$ which is inside $\Gamma$. Cut the elements of $T_2^h$ intersecting $\Gamma$ and discard all elements completely outside $\Gamma$. The solution on these two meshes is then joined by Nitsche's method.

In the next section we prove that indeed $U^* = U|_{\Omega}$ and hence the analysis of [3] carries over to our formulation.

## 3. A PRIORI ERROR ESTIMATES

Consider the following mesh-dependent norms:

$$\|v\|_{1/2, h, \Gamma}^2 := \|h(x)^{-1/2}v\|_{0, \Gamma}^2 = \sum_{j \in J_h} h_{K_1^j}^{-1} \|v\|_{0, \Gamma^j}^2,$$

$$\|v\|_{-1/2, h, \Gamma}^2 := \|h(x)^{1/2}v\|_{0, \Gamma}^2 = \sum_{j \in J_h} h_{K_1^j}^{1/2} \|v\|_{0, \Gamma^j}^2,$$
and

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

**Proposition 3.1**

Let $U$ denote the solution of (4) and $U^*$ the solution of (5). Then $U|_{\Omega}=U^*$.

**Proof**

The proof proceeds in three steps

1. Show the existence of a subspace $\tilde{V}_h \subset V_h$ such that

$$ a^h(U, \tilde{\phi}) = a_*^h(U, \tilde{\phi}) \quad \text{and} \quad l(\tilde{\phi}) = l_*(\tilde{\phi}) \quad \forall \tilde{\phi} \in \tilde{V}_h. $$

2. Show that $\exists \xi_U \in \tilde{V}_h$ such that $\xi_U|_{\Omega}=U$ and that $\exists \xi_{U^*} \in \tilde{V}_h$ such that $\xi_{U^*}|_{\Omega}=U^*$.

3. Apply coercivity and Galerkin orthogonality to the discrete error using the results of [3].

First, let

$$ \tilde{V}_h := \{ \tilde{\phi} \in V_h : \tilde{\phi}_2|_{\hat{\Omega}\setminus(\Omega_2 \cup \Omega_{G_h}^{\text{ext}})} = 0 \text{ and } \tilde{\phi}_4|_{\Omega_{G_h}^{\text{ext}}} = -\tilde{\phi}_2|_{\Omega_{G_h}^{\text{ext}}} \}. $$

Since the integrals on $\Gamma$ are the same in $a^h(\cdot, \cdot)$ and $a_*^h(\cdot, \cdot)$ we only need to prove the equivalence of the volume integrals. By using the decomposition $\hat{\Omega} = \hat{\Omega} \setminus (\Omega_2 \cup \Omega_{G_h}^{\text{ext}}) \oplus (\Omega_2 \setminus \Omega_{G_h}^{\text{int}}) \oplus \Omega_{G_h}^{\text{int}} \oplus \Omega_{G_h}^{\text{ext}}$, we may write for all $\phi \in V_h$

$$ \begin{aligned} a^h(U, \phi) ={}& (\nabla U_1, \nabla \phi_1)|_{\Omega_1} + (\nabla U_2, \nabla \phi_2)|_{\hat{\Omega} \setminus (\Omega_2 \cup \Omega_{G_h}^{\text{ext}})} + (\nabla U_2, \nabla \phi_2)|_{\Omega_2 \setminus \Omega_{G_h}^{\text{int}}} \\ & + (\nabla (U_2+U_4), \nabla (\phi_2+\phi_4))|_{\Omega_{G_h}^{\text{ext}}} + (\nabla (U_2+U_3), \nabla (\phi_2+\phi_3))|_{\Omega_{G_h}^{\text{int}}} + B_{\Gamma}(U, \phi), \quad (6) \end{aligned} $$

where $B_{\Gamma}(U, \phi)$ denotes the integrals over $\Gamma$. It then follows that for all $\tilde{\phi} \in \tilde{V}_h$

$$ \begin{aligned} a^h(U, \tilde{\phi}) ={}& (\nabla U_1, \nabla \tilde{\phi}_1)|_{\Omega_1} + (\nabla U_2, \nabla \tilde{\phi}_2)|_{\Omega_2 \setminus \Omega_{G_h}^{\text{int}}} + (\nabla (U_2+U_3), \nabla (\tilde{\phi}_2+\tilde{\phi}_3))|_{\Omega_{G_h}^{\text{int}}} + B_{\Gamma}(U, \tilde{\phi}) \\ ={}& (\nabla U_1, \nabla \tilde{\phi}_1)|_{\Omega_1} + (\nabla (U_2+U_3), \nabla (\tilde{\phi}_2+\tilde{\phi}_3))|_{\Omega_2} + B_{\Gamma}(U, \tilde{\phi}) = a_*^h(U, \tilde{\phi}). \end{aligned} \quad (7) $$

The equality $l(\tilde{\phi})=l_*(\tilde{\phi})$ is shown in a similar fashion.

Second, observe that $\tilde{V}_h$ only imposes constraints on components of the function outside $\Omega$, since the constraint on $\phi_2$ in $\Omega_{G_h}^{\text{int}}$ is compensated for by the freedom of $\phi_3$. Hence, the existence of the sought $\xi_U$ and $\xi_{U^*}$.

Finally we recall the following coercivity result from [3], for some $c>0$ and for all $v \in V_h$ there holds

$$ c |||v|||^2 \le a_*^h(v, v). \quad (8) $$

In particular this holds for $v=U-U^*$ and hence, since $\xi_U - \xi_{U^*} \in \tilde{V}_h$

$$ \begin{aligned} c |||U - U^*|||^2 \le {}& a_*^h(U - U^*, U - U^*) = a_*^h(U - U^*, \xi_U - \xi_{U^*}) \\ = {}& a^h(U, \xi_U - \xi_{U^*}) - a_*^h(U^*, \xi_U - \xi_{U^*}) = l(\xi_U - \xi_{U^*}) - l_*(\xi_U - \xi_{U^*}) = 0. \end{aligned} $$

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

We have the following consistency relation.

**Proposition 3.2**

The discrete problem (4) is consistent in the sense that, for $u$ solving (2) there holds

$$ a_*^h(u, \phi) = l_*(\phi) \quad \forall \phi \in V^h, $$
or, equivalently,

$$a_{*}^{h}(u - U, \phi) = 0 \quad \forall \phi \in V^{h}. \qquad (9)$$

The proof is given in [3].

We can now directly take advantage of the theory developed in [3] which shows that we have optimal error estimates for any polynomial degree of the underlying FEM:

**Theorem 3.3**

For $U$ solving (4) and $u$ solving (2), the following a priori error estimates hold

$$\|u - U\| \leq C h^p |u|_{p+1, \Omega}, \qquad (10)$$

and

$$\|u - U\|_{0, \Omega} \leq C h^{p+1} |u|_{p+1, \Omega}. \qquad (11)$$

We refer to [3] for the proof. Here we shall only point out one of the crucial points in the analysis. Accuracy of the method is expressed by the orthogonality relation (9), but to show convergence we also need stability of the discrete problem as expressed by (8). In order to show that our system matrix is positive definite, we rely on the following inverse inequality, see, e.g. Warburton and Hesthaven [9].

**Lemma 3.4**

For $\phi \in V^h$, the following inverse inequality holds:

$$\|\partial_n \phi_1\|_{-1/2, h, \Gamma} \leq C_I p \|\nabla \phi\|_{0, \Omega_1}.$$

The size of the constant $C_I$ can be found by solving a small local eigenvalue problem: explicit bounds are discussed in [9]. The constant $C_I$ determines the size of $\lambda$, we are obliged to take $\lambda > C_I^2 p^2$ in order to ensure coercivity. Here is the point of the band: if we just cut the mesh with $\partial\Omega$ and apply Nitsche's method on the cut elements, sliver elements would be generated that would require extremely large $C_I$ in Lemma 3.4, leading to severe ill conditioning of the discrete system. We illustrate the problem in the case $p=1$. Thus, consider using only $V_1^h$ and cutting the mesh with $\partial\Omega$. Then $C_I$ can be found, for a given element $K$ as the largest eigenvalue $\lambda_{max}$ in the eigenproblem of finding $U \in V_1^h$ and $\lambda \in \mathbb{R}$ such that

$$(h_K^{1/2} \partial_n U, \partial_n v)_{K \cap \partial \Omega} = \lambda (\nabla U, \nabla v)_{K \cap \Omega} \quad \forall v \in V_1^h.$$

For $p=1$, $\nabla v$ is constant on each element and thus we have

$$\|\partial_n v\|_{L_2(K \cap \partial \Omega)}^2 = \mathrm{meas}(K \cap \partial \Omega) |\partial_n v|^2, \qquad (12)$$

where meas($\cdot$) denotes the length, area, or volume of the object in question, and

$$\|\nabla v\|_{L_2(K \cap \partial \Omega)}^2 \geq \|\partial_n v\|_{L_2(K \cap \partial \Omega)}^2 = \mathrm{meas}(K \cap \partial \Omega) |\partial_n v|^2 \qquad (13)$$

and it follows that

$$\|h_K^{1/2} \partial_n v\|_{L_2(K \cap \partial \Omega)}^2 \leq \frac{h_K \mathrm{meas}(K \cap \partial \Omega)}{\mathrm{meas}(K \cap \partial \Omega)} \|\nabla v\|_{L_2(K \cap \partial \Omega)}^2, \qquad (14)$$

and thus we must choose

$$\lambda > C_I^2 \geq \frac{h_K \mathrm{meas}(K \cap \partial \Omega)}{\mathrm{meas}(K \cap \partial \Omega)}.$$
<figure><img src="image_6.png" /><figcaption>Figure 4. An element cut by the boundary.</figcaption></figure>

The situation is illustrated in Figure 4 where $\partial\Omega$ is represented by the dashed line. We note that as $\partial\Omega$ moves to the left in Figure 4, the area meas$(K \cap \partial\Omega)$ will approach zero, whereas the length meas$(K \cap \partial\Omega)$ remains bounded from below, which means that $\lambda$ must grow without bound. This situation is remedied by inserting the band of elements between the cut and the boundary.

In Burman and Hansbo [10], where no band was used, this conditioning problem was instead eliminated by the use of additional stabilizing terms, limiting the analysis to linear polynomials. In the present formulation, we avoid the use of additional terms since the mesh on the band $\Omega_1$ is always shape regular. We also remark in this context that the use of higher-order polynomials does not require a curved interface since the interface is artificial. Thus, the interface can always be assumed to be piecewise affine and no further complexity is introduced by higher-order approximations, cf. also [3].

## 4. IMPLEMENTATION ASPECTS

The difference between the suggested approach and that of [3] is in the new hierarchical formulation. Conceptually, this means that in the method of [3], there are only two finite element spaces: the one on the band and the one inside the band consisting of elements up to the interface $\Gamma$. This gives two sets of unknowns, $u_1$ and $u_2$ say, and the system of equations becomes

$$ \begin{bmatrix} S_1 & B \\ B^T & S_2 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \end{bmatrix} $$

where $S_1$ is the stiffness matrix from the discretization on the band, $S_2$ from the cut mesh, $B$ represents the coupling terms, and $f_1, f_2$ represents load terms. The degrees of freedom in $u_2$ outside the cut can then be discarded already at the outset. One problem with this approach is that all matrices $B, S_1$, and $S_2$ change if we want to move the interface. This could be in a time-dependent problem, or if we want to use the scheme for the purpose of shape optimization. In the present approach we instead obtain the system

$$ \begin{bmatrix} S_1 & B_{12} & B_{13} & B_{14} \\ B_{12}^T & S_2 & B_{23} & B_{24} \\ B_{13}^T & B_{23}^T & S_3 & B_{34} \\ B_{14}^T & B_{24}^T & B_{34}^T & S_4 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \\ f_3 \\ f_4 \end{bmatrix} $$
where $S_2$ now denotes the stiffness matrix on the underlying mesh which does not change. Since this matrix is by far the largest of the involved matrices, this means that we only have to recompute small matrices (corresponding in a sense to the multiplier matrices of the original fictitious domain method).

In our implementation, two different approximations of the interface have been used: the first defining the approximation of $\Gamma$ given by the piecewise affine (internal) boundary of $T_1^h$, the second defined by the approximation of the intersection of $\Gamma$ with the elements on $T_2^h$ (which defines a different affine approximation of $\Gamma$). The first has been used as a master surface for the integration of jump and consistency terms and should also, for exactness, be used for the computation of the bilinear form $a^h(\cdot, \cdot)$. We have however, for simplicity, used the second affine approximation of the interface for computing $a^h(\cdot, \cdot)$ which removes some technical difficulties in the implementation but introduces a conformity error. This error does not, however, affect the convergence rates in our example below. For higher-order methods, the intersection of the meshes must be handled more carefully, although, as noted above, the interface can remain piecewise affine since it is given by the user-defined internal boundary of $T_1^h$.

<figure><img src="image_7.png" /><figcaption>Figure 5. Convergence in broken $H^1$ and $L_2$ norms.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Elevation of the computed solution on $\Omega$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7. Elevation of the solution on $\hat{\Omega}$.</figcaption></figure>

## 5. NUMERICAL EXAMPLE

We consider a problem posed on a disc of radius $r_0=0.95$. With $r$ the length of the radius vector, we use $f=r$ to obtain the exact solution $u=u=(r_0^3-r^3)/9$. The stabilization parameter was set to $\gamma=10$. In Figure 5 we show the obtained convergence rates in $L_2(\Omega)$- and $H^1(\Omega)$-norms, which are optimal. An elevation of the solution is given in Figure 6, and an elevation of the solution on the whole of $\hat{\Omega}$ is shown in Figure 7. Note that we have extended $f=r$ to hold on the whole of $\hat{\Omega}$ and imposed zero boundary conditions on $\partial\hat{\Omega}$. This is of no consequence since the solution is decoupled at $\Gamma$.

## 6. CONCLUDING REMARKS

In this contribution we have shown that the NXFEM method is well suited for fictitious domain-type simulations. It has optimal convergence for arbitrary polynomial order and does not require Lagrange multipliers to enforce Dirichlet boundary conditions. Indeed, since the boundary conditions are prescribed on a regular mesh, we can handle all types of boundary conditions in the usual way.

Applications for the proposed method include shape optimization, where the boundary of the domain has to be moved in order to calculate sensitivities, and for computations involving objects moving across a background mesh in general.

## REFERENCES

1. Glowinski R, Pan TS, Periaux J. A fictitious domain method for Dirichlet problem and applications. *Computer Methods in Applied Mechanics and Engineering* 1994; **111**:283–303.

2. Girault V, Glowinski R. Error analysis of a fictitious domain method applied to a Dirichlet problem. *Japan Journal of Industrial and Applied Mathematics* 1995; **12**:487–514.

3. Hansbo A, Hansbo P, Larson MG. A finite element method on composite grids based on Nitsche's method. *ESAIM: Mathematical Modelling and Numerical Analysis* 2003; **37**:495–514.

4. Becker R, Hansbo P, Stenberg R. A finite element method for domain decomposition with non-matching grids. *ESAIM: Mathematical Modelling and Numerical Analysis* 2003; **37**:209–225.

5. Hansbo A, Hansbo P. An unfitted finite element method, based on Nitsche's method, for elliptic interface problems. *Computer Methods in Applied Mechanics and Engineering* 2002; **191**:5537–5552.
6. Nitsche J. Über ein Variationsprinzip zur Lösung Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. *Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg* 1971; **36**:9–15.

7. Areias PMA, Belytschko T. A comment on the article: a finite element method for simulation of strong and weak discontinuities in solid mechanics. *Computer Methods in Applied Mechanics and Engineering* 2006; **195**:1275–1276.

8. Maury B. A fat boundary method for the Poisson problem in a domain with holes. *Journal of Scientific Computing* 2001; **16**:319–339.

9. Warburton T, Hesthaven JS. On the constants in hp-finite element trace inverse inequalities. *Computer Methods in Applied Mechanics and Engineering* 2003; **192**:2765–2773.

10. Burman E, Hansbo P. Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method. Preprint, Department of Mathematical Sciences, Chalmers University of Technology, 2010.