Analysis of optimal preconditioners for CutFEM

In this paper we consider a class of unfitted finite element methods for scalar elliptic problems. These so-called CutFEM methods use standard finite element spaces on a fixed unfitted triangulation combined with the Nitsche technique and a ghost penalty stabilization. As a model problem we consider the application of such a method to the Poisson interface problem. We introduce and analyze a new class of preconditioners that is based on a subspace decomposition approach. The unfitted finite element space is split into two subspaces, where one subspace is the standard finite element space associated to the background mesh and the second subspace is spanned by all cut basis functions corresponding to nodes on the cut elements. We will show that this splitting is stable, uniformly in the discretization parameter and in the location of the interface in the triangulation. Based on this we introduce an efficient preconditioner that is uniformly spectrally equivalent to the stiffness matrix. Using a similar splitting, it is shown that the same preconditioning approach can also be applied to a fictitious domain CutFEM discretization of the Poisson equation. Results of numerical experiments are included that illustrate optimality of such preconditioners for the Poisson interface problem and the Poisson fictitious domain problem.

1. Introduction. In recent years many papers appeared in which the so-called CutFEM paradigm is developed and analyzed, cf. the overview references [5,3]. In this approach, for discretization of a partial differential equation a fixed unfitted mesh is used that is not aligned with a (moving) interface and/or a complex domain boundary. On this mesh standard finite element spaces are used. For treating the boundary and/or interface conditions, either a Lagrange multiplier technique or Nitsche's method is applied. In the setting of the present paper we restrict to Nitsche's method. Furthermore, to avoid extreme ill-conditioning of the resulting discrete systems (due to "small cuts") a stabilization technique is used. The most often used approach is the ghost-penalty stabilization [4]. In the literature the different components of this general CutFEM are studied, error analyses are presented and different fields of applications are studied [5,3]. Related unfitted finite element methods are popular in fracture mechanics [14]; in that community these are often called extended finite element methods (XFEM).
Almost all papers on CutFEM (or XFEM) either treat applications of this methodology or present discretization error analyses. In relatively very few papers efficient solvers for the resulting discrete problems are studied. In [7,32,19], for the resulting stiffness matrix condition number, bounds of the form ch −2 , with h a mesh size parameter and c a constant that is independent of how an interface or boundary intersects the triangulation, have been derived. In [7] a fictitious domain variant of CutFEM is introduced and it is shown that discretization of a Poisson equation using this method yields a stiffness matrix with such a condition number bound. In [32] a similar result is derived for CutFEM applied to a Poisson interface problem. In [19] a condition number bound is derived for CutFEM applied to a Stokes interface problem. These papers do not treat efficient preconditioners for the stiffness matrix.
There are few papers in which (multigrid type) efficient preconditioners for Cut-FEM or closely related discretizations (e.g., XFEM) are treated, e.g., [2,1,11,21,10,26,25]. In none of these papers a rigorous analysis of the spectral quality of the preconditioner is presented. The only paper that we know of that contains such a rigorous analysis is [24]. In that paper a CutFEM without stabilization is analyzed for a two-dimensional Poisson interface problem.
The main topic of the present paper is an analysis of a (new) subspace decomposition preconditioning technique for a CutFEM discretization of elliptic interface problems and for a CutFEM fictitious domain method. These discretization methods are known in the literature and are typical representatives of the CutFEM methodology [6,7,27]. This preconditioning technique leads to very natural and optimal preconditioners, in a sense as explained in section 6. We expect that similar preconditioners can be developed and rigorously analyzed for other CutFEM applications such as a Stokes fictitious domain method and Stokes interface problems.
We explain the key idea of the preconditioner for the interface problem. In the CutFEM applied to such an elliptic interface problem one uses a standard H 1conforming finite element space on a triangulation that is not fitted to the interface. For treating the interface conditions a Nitsche technique is used, leading to additional bilinear forms (consistency and penalty terms) in the variational formulation of the discrete problem. To damp the instabilities due to "small cuts" a ghost-penalty stabilization term is also added in the discrete variational formulation. The finite element space used in the CutFEM has a natural splitting into two subspaces, a "global" and a "local" one. The global subspace is spanned by all standard nodal basis functions on the whole triangulation, and the local space is spanned by nodal cut basis functions "close to" the interface. We will show that this splitting is stable, uniformly in the discretization parameter h and in the location of the interface in the triangulation. We also prove that the Galerkin discretization in the local subspace leads (after diagonal scaling) to a uniformly well-conditioned matrix and that the Galerkin discretization in the global subspace is uniformly equivalent to the standard finite element discretization of the Poisson interface problem on the global domain. Using the latter property it follows that a multigrid method yields an optimal preconditioner for the Galerkin discretization in the global subspace. An additive Schwarz subspace correction method (or, equivalently, block Jacobi) thus yields an optimal preconditioner for the CutFEM discretization of the interface problem. The same approach applies, with minor modifications, to a CutFEM fictitious domain discretization of scalar elliptic problems.
We briefly address relations between the results in this paper and in [24]. In the latter a CutFEM variant without stabilization is studied and the preconditioner is based on a subspace splitting that is similar to the one studied in this paper. The rather technical analysis in [24] is restricted to linear finite elements and twodimensional problems. In this paper we consider the CutFEM with stabilization. It turns out that this allows an elegant, rather simple and much more general analysis. In particular, the analysis covers two-and three-dimensional problems, arbitrary polynomial degree finite elements and triangulations that are shape regular but not necessarily quasi-uniform. Furthermore, the analysis of this paper can also be applied to related CutFEM discretizations such as, for example, the CutFEM fictitious domain method. A preliminary preprint version of this paper, in which only the preconditioner for the CutFEM fictitious domain method is treated, is [16].
The paper is organized as follows. In Section 2 we describe a CutFEM discretiza-2 tion of elliptic interface problems known from the literature. In Section 3 two related matrix-vector representations of the discrete problem are introduced. In Section 4 several uniform norm equivalences are derived that are used in Section 5 to prove a stable splitting property. Based on this stable splitting we propose (optimal) preconditioners in Section 6. In Section 7 results of numerical experiments with these preconditioners are presented.

2.
CutFEM for interface problems. We recall a class of CutFEM methods known from the literature [18,5,19]. On a bounded connected polygonal domain Ω ⊂ R d , d = 2, 3, we consider the following standard model problem for scalar elliptic interface problems: (2.1) Here, f ∈ L 2 (Ω) is a given source term, Ω 1 ∪ Ω 2 = Ω a non-overlapping partitioning of the domain, Γ = Ω 1 ∩ Ω 2 is the interface, [[·]] denotes the usual jump operator across Γ and n Γ denotes the unit normal at Γ pointing from Ω 1 into Ω 2 . The weak formulation of the problem (2.1) is as follows: determine u ∈ H 1 0 (Ω) such that Here and in the remainder, (·, ·) Ω denotes the L 2 scalar product on Ω. We assume that for discretization a family of simplicial triangulations {T h } h>0 of Ω is used which are not fitted to Γ. Let T h denote a simplicial triangulation of Ω and V k h the corresponding standard finite element space of continuous piecewise polynomials up to degree k that have zero values on ∂Ω. The set of all simplices that are cut by the interface Γ is denoted by T Γ h and the domain formed by these simplices is denoted by Ω Γ h . The domain formed by all simplices with nonzero intersection with Ω i ("extended subdomain") is denoted by Ω ex i,h , i = 1, 2. Note that Ω Γ h ⊂ Ω ex i,h holds. In the CutFEM one uses pairs of finite element functions u h : Based on this space we formulate a discretization of (2.1) using the Nitsche technique: Here  [18,9,29] or the overview in [25]. For the case of linear finite elements optimal discretization error bounds for this method are derived in [18]. For the higher order case, but without the ghost-penalty term G h (·, ·), optimal discretization error bounds are derived in [29,30]. These analyses can be extended to the case with the ghost-penalty stabilization.
Since we do not assume quasi-uniformity of the triangulation, the scalings with h and with h −1 are element-wise, e.g., The parameters γ > 0, β > 0 are fixed. The bilinear form G h (·, ·) is the ghost penalty stabilization. Different equivalent variants of this stabilization are known in the literature, cf. [4,28,23]. The choice of a particular variant of this stabilization is not relevant for the analysis in this paper.
Remark 1. In practice the method above is not feasible because integrals over cut simplices T ∩ Ω i and over the interface segments T ∩ Γ are difficult to compute. For linear finite elements (k = 1) one usually replaces Γ by a suitable piecewise linear approximation Γ h . For higher order finite elements the isoparametric approach introduced in [22] can be used. In that approach one assumes that the interface is represented as the zero level of a level set function. The fundamental idea is the introduction of a (level set function based) parametric mapping Θ h of the underlying mesh from a geometrical reference configuration to a final configuration, cf. Fig.  2.1. We refer to [22] for the definition of Θ h . The discretization approach consists  [22]: The piecewise linear approximation Γ lin is mapped to a higher order approximation Γ h using a mesh transformation Θ h . of two steps. First, a (higher order) finite element space is considered with respect to the reference configuration. Then the transformation Θ h is applied to this space and to the geometries in the variational formulation, resulting in a new unfitted finite element discretization with an accurate treatment of the geometry. The mapping renders the finite element spaces into isoparametric finite element spaces. The mapping Θ h and corresponding quadrature rules are implemented in the add-on library ngsxfem to Netgen/NGSolve. The isoparametric Nitsche unfitted FEM is a transformed version of the original Nitsche unfitted FE discretization [18] with respect to the interface approximation Γ h = Θ h (Γ lin ), where Γ lin is the zero level of a piecewise linear interpolation of a sufficiently accurate higher order finite element approximation of the level set function φ. In the isoparametric approach one uses the spaces For further explanation of this method and its discretization error analysis we refer to [29,30]. We will not consider this "perturbation" due to the isoparametric transformation because it makes the presentation of the analysis below less transparent. We restrict to the method with exact geometry approximation as defined in (2.3) since this "geometric error" does not play an essential role with respect to the spectral accuracy of the preconditioner introduced in this paper.
It turns out that the preconditioner that we treat in this paper can easily be modified for application to a CutFEM applied in a fictitious domain approach. To explain this more precisely, we describe in the remark below a Nitsche fictitious domain discretization known from the literature. The corresponding preconditioner for this problem is discussed in Remark 7.
Remark 2. Instead of the interface problem (2.1) we consider the Poisson equation For discretization we apply a fictitious domain method known from the literature [7,27]: where the bilinear form is defined by Here the Nitsche method is used to satisfy (approximately) the boundary condition u = g on Γ, whereas for the interface problem the Nitsche method is used to enforce the interface condition [[u]] = 0 on Γ. The discretization of the interface problem can be seen as a fictitious domain discretization "from both sides" Ω i , i = 1, 2, with a coupling condition u 1 | Γ = u 2 | Γ on Γ.
3. Discrete problems in matrix vector formulation. In this section we briefly recall two matrix vector formulations of the discretization (2.3). We introduce the (closed) subdomains formed by all simplices that are completely contained in The finite element nodal basis functions of V h := V k h are denoted by φ j , j ∈ I 0 , for a suitable index set I 0 . The finite element nodes that are in Ω ex i,h are labeled by j ∈ I i ⊂ I 0 , i = 1, 2. Let I Γ ⊂ I 0 be the subset of labels corresponding to finite element nodes in Ω Γ h and To simplify the presentation, we assume that there are no nodes on ∂Ω i . Note that Note that for k ≥ 2 and interior nodes (i.e., nodes strictly inside an element Using this basis one obtains a matrix-vector representation of (2.3) that is denoted by Ax = b. For preconditioning it is convenient to use another representation of the discrete solution, namely as a suitable global finite element function in the space V h that is corrected using a finite element function with support only on Ω Γ h . This representation is more in the spirit of the extended finite element method. Fig. 3.1: Sketch of interface Γ and (part of) triangulation T h of Ω h with interface nodes I Γ 1 (blue circles) and I Γ 2 (red rectangles). All rectangular nodes form the set More precisely, we introduce the local spaces The bases used in (3.1) and (3.2) are the same, hence 6 We introduce the projections With the compact notationû h : cf. Figure 3.2, or in basis notation The discrete problem can be reformulated in the space The corresponding matrix vector problem is denoted bŷ In the remainder we introduce and analyze a preconditioner for this discrete problem. The matrix representation L of the isomorphism (3.4) is simple, as for the part corresponding to φ j , j ∈ I 0 \ I Γ and φ Γ j , j ∈ I Γ it is only a permutation matrix, and for the remaining part (φ j , j ∈ I Γ ) there are exactly two non-zero entries per column. To understand the latter, consider an index j ∈ I Γ 2 . Then also j ∈ I 1 \ I Γ 1 and, hence, L(φ j , 0) T = (φ j , φ Γ j ) T . This implies that given A, the matrixÂ = L T AL is easy to obtain by permutation and summation of corresponding pairs of rows and columns. Also L −1 is easy to determine. A spectrally equivalent preconditioner PÂ ofÂ induces a corresponding spectrally equivalent preconditioner L −T PÂL −1 of A.
Remark 3. In case of the fictitious domain discretization (2.5) of the Poisson problem (2.4) on Ω 1 , an analogous splitting of the fictious finite element space V FD h is given by have Note that elements from the global space V − 1 are zero on ∂Ω ex 1,h . 4. Fundamental norm equivalences. In this section several norm equivalences are presented that will be used to derive a new spectral equivalence result for the bilinear formÂ h (·, ·) in the main theorem 5.2 below. We start with a norm equivalence, which is known from the literature. We use the notation ∼ to denote estimates in both directions with constants that are independent of h and of the location of the interface Γ in the triangulation. We recall the notation introduced above: for From the literature on discretization error analyses of CutFEM, e.g. [27], the following fundamental norm equivalence is known: For this uniform norm equivalence to hold it is essential that a ghost penalty type stabilization is added. We derive preliminaries in the following lemmas. We will use the trace inequality [18]: For a subdomain ω ⊂ Ω we use the notation The result in the next lemma gives a useful uniform norm equivalence for finite element functions restricted to the local interface strip Ω Γ h . Lemma 4.1. The following uniform norm equivalence holds: Proof. Using (4.2) we get Combining this with a standard finite element inverse inequality yields i.e., a uniform estimate in one direction in (4.3). We now derive the estimate in the other direction. We introduce, for T ∈ T Γ h , the subdomain consisting of all simplices in T Γ h that have at least a common vertex with T , i.e., ω T := {T ∈ T Γ h |T ∩ T = ∅ }. Note that due to shape regularity we have hT ∼ h T forT ∈ ω T and diam(ω T ) ∼ h T . Take The area |T ∩ Γ| can be arbitrary small ("small cuts"), but it follows from [12,Proposition 4.2] that there is an elementT ∈ ω T such that |T ∩ Γ| ≥ c 0 h d−1 T , with a constant c 0 > 0 that depends only on shape regularity of T Γ h and on smoothness of Γ. Take such anT ∈ ω T . Take a fixed ξ ∈ Γ ∩T such that |v h (ξ)| = max x∈T ∩Γ |v h (x)| =: v h ∞,T ∩Γ . Take x ∈ T and let S be a smooth shortest curve in ω T that connects x and ξ. Due to shape regularity we have |S| h T , independent of x. This yields with s the arclength parametrization of S. Hence, Using integration over T , |T | ∼ h d T and the standard FE norm estimate ∇v h and combining this with the result (4.4) and hT ∼ h T yields Summing over T ∈ T Γ h completes the proof. Remark 4. Results similar to (4.3) are known in the literature. For example, in the papers [8,15], for the case of a quasi-uniform triangulation the following uniform estimate is derived: Note that due to the quasi-uniformity assumption we have a simpler scaling with the global mesh parameter h and that in (4.5) we have the normal derivative term n · ∇v h Ω Γ h , with n the normal on Γ (constantly extended in the neighborhood Ω Γ h ) instead of the full derivative term ∇v h Ω Γ h . The proofs of (4.5) in [8,15] are much more involved than the simple proof of Lemma 4.1 above. This is due to the fact that in the bound in (4.5) only the normal derivative occurs.
A second norm equivalence is derived in the following lemma.
The estimate in the one direction directly follows from a standard finite element inverse inequality.
and thus which is this estimate in the other direction.
Thus we obtain the following corollary.
Corollary 4.3. The following uniform norm equivalence holds Besides this norm equivalence result for finite element functions from the local correction spaces v Γ i ∈ V Γ i , i = 1, 2, there also holds a strengthened Cauchy-Schwarz inequality for the two spaces V Γ i , i = 1, 2. This is shown in Lemma 4.5. For the proof of that lemma it is convenient to use the following elementary estimate.
which proves the result. Using this we obtain the following uniform strengthened Cauchy-Schwarz inequality and a corresponding norm equivalence.
Note that β, γ = 0 holds. Let M ∈ R m×m , M i,j = (φ i , φ j ) T , be the element mass matrix. Note that κ(M ) = κ(M ) holds. Thus we obtain, using Lemma 4.4: which yields the result (4.8). Multiplying by h −2 T and summing over T ∈ T Γ h we get which yields the estimate (4.9) in one direction with constant κ(M ) − 1 2 . The estimate in the other direction follows from the triangle inequality.

Stable subspace splitting.
Based on results from the previous section we now derive a stable splitting result which essentially states that the angles (in the energy scalar product) between the subspaces V h , V Γ in V h × V Γ are uniformly bounded away from zero. Based on classical theory cf. [17,31] this then immediately leads to optimal block-Jacobi type preconditioners. We recall three norm equivalences from the previous section that we need to derive the stable splitting property, namely the ones in (4.3), (4.7) and (4.9): with notation as in (3.3). Forû h = (u 0 , u Γ ) ∈ V h × V Γ , projections P i on the two subspaces are defined by P 0ûh := (u 0 , 0), P 1ûh := (0, u Γ ).
Lemma 5.1. The following uniform norm equivalence holds Hence the result (5.4) holds. Theorem 5.2. The following uniform norm equivalences hold: Proof. The result in (5.6) is a direct consequence of (5.5) and (4.1). We prove the result (5.5) as follows: where in the last step we used Lemma 5.1. From this and P 0ûh = (u 0 , 0), P 1ûh = (0, u Γ ) the result (5.5) follows. Remark 5. With similar arguments as in the proof of Theorem 5.2 one can show that the norm equivalence and V Γ 2 is stable. However, concerning preconditioning this does not yield significant advantages compared to the stable splitting of V h × V Γ in the subspaces V h and V Γ .
Remark 6. The constants in ∼ in (5.5)-(5.6) will depend on the jump in the diffusion coefficient α across the interface. Therefore, the preconditioners proposed in the next section are not expected to be robust with respect to large jumps in this coefficient. We expect that robustness can be obtained using suitable scalings in (5.5)-(5.6) that depend on the diffusion coefficient. This will be analyzed in future work.
A stable subspace splitting result similar to (5.5) also holds for the fictitious domain bilinear form with subspaces V − 1 and V Γ 1 , cf. Remarks 2 and 3. On we define the energy norms with notation as in (3.7). From the literature [7,27] we have (for γ sufficiently large) the uniform norm equivalence Along the same lines as in the proof of (5.5) with u 0 replaced by u − 1 , Q 1 u Γ replaced by u Γ 1 and u 2,h = Q 2 u Γ = 0 one obtains forû 1,h Thus we get the uniform norm equivalence which yields the stable subspace splitting result for the fictitious domain method.
6. Optimal preconditioners. We return to the linear systemÂx =b in (3.6). We introduce some notation to represent the subspace splitting in matrix-vector format. The coefficient vector x that represents the unknown finite element function u h = (u 0 , u Γ ) is split into the parts corresponding to u 0 and u Γ , i.e., x = (x 0 , x 1 ) with We define corresponding projections P i by P 0 x = (x 0 , 0), P 1 x = (0, x 1 ). The Galerkin projections on the subspaces are denoted byÂ i , i.e., we have the relations Let D A := blockdiag(Â 0 ,Â 1 ) be the blockdiagonal matrix corresponding to the Galerkin projections on the subspaces. The result (5.6) in matrix formulation yields that D A is spectrally equivalent toÂ: Hence D A is an optimal preconditioner forÂ in the sense that the spectral condition number λ max (D −1 AÂ )/λ min (D −1 AÂ ) is uniformly bounded both with respect to the mesh size h and the location of Γ in the triangulation. Note that this condition number may depend on the size of the jumps in the diffusion coefficient α, cf. Remark 6.
Clearly the preconditioner D A , which we call the exact preconditioner, is not computationally efficient. We now explain how the diagonal blocksÂ i , i = 0, 1, can be replaced by computationally efficient spectrally equivalent approximations, which then yields a computationally efficient optimal preconditioner forÂ.
We first consider the blockÂ 0 that corresponds to the Galerkin projection onto the global H 1 0 (Ω)-conforming finite element space V h . We have It is natural to consider a spectrally equivalent preconditioner, denoted by B 0 , for the interface problem (2.2) discretized in the standard conforming finite element space An option for such a B 0 is a multigrid preconditioner. From (6.1) it follows that B 0 is then also uniformly spectrally equivalent toÂ 0 , i.e., B 0 ∼Â 0 .
We finally consider computationally efficient optimal preconditioners for the block A 1 , which corresponds to the local correction space V Γ . Lemma 6.1. For D 1 := diag(Â 1 ) the uniform spectral equivalence

holds.
Proof. From Lemma 5.1 it follows thatÂ 1 is spectrally equivalent to a mass matrix and it is well-known that the diagonally scaled mass matrix has a uniformly bounded spectral condition number. For completeness we give the details. Recall the relation between x 1 = (x 1,j ) j∈I Γ and u Γ that is given by u Γ = j∈I Γ x 1,j φ Γ j . Using Lemma 5.1 we get For T ∈ T Γ h we denote by N (T ) ⊂ I Γ the subset of indices with corresponding nodes in T . Standard arguments yield that u Γ 2 T ∼ |T | j∈N (T ) x 2 1,j holds. Using this we get For j ∈ I Γ defineû j := (0, φ Γ j ). Hence, (D 1 ) j,j = (Â 1 ) j,j = P 1ûj 2 a and T |T | and thus we get Comparing (6.2) and (6.3) we obtain the spectral equivalence. Remark 7. Based on the stable splitting result (5.8) the same approach can be applied to derive optimal block Jacobi preconditioners for the fictitious domain discretization. In that case theÂ 0 "global" block corresponds to a finite element discretization of the Laplace problem in V − 1 := span{φ j } j∈I1\I Γ 1 with homogeneous Dirichlet boundary condition on the boundary of the domain formed by these basis functions. As spectrally equivalent preconditioner B 0 for this block one can again use a multigrid solver. The other diagonal blockÂ 1 corresponds to Galerkin discretization in V Γ 1 and the result in Lemma 6.1 implies that the diagonally scaled version of this matrix has a uniformly bounded condition number.
7. Numerical experiments. The analysis above leads to the following preconditioners for the linear system in (3.6). Preconditioners ofÂ i are denoted by B i , i = 0, 1. We define the block Jacobi preconditioners Here, we choose for B −1 0 a few (3) multigrid sweeps (V-cycle) with symmetric Gauss-Seidel smoothing applied toÂ 0 and for B −1 1 the symmetric Gauss-Seidel preconditioner (one iteration) applied toÂ 1 . In the following, we apply a preconditioned conjugate gradient (PCG) method to the linear system (3.6) and examine different choices of preconditioners P. Starting with x 0 = 0, the PCG iteration is stopped when the preconditioned residual is reduced by a factor tol = 10 −6 , i.e.
with · 2 the Euclidean norm. In the following, results for the Poisson interface problem and Poisson fictitious domain problem are presented. All numerical experiments have been performed with the DROPS package [13].
We present results for the symmetric Gauss-Seidel preconditioner P SGS and the block Jacobi preconditioners P A , P D , P B defined in (7.1). and PCG iteration numbers for different refinement levels are reported in Table 7.3.
For finer grid levels ≥ 4 the condition number κ 2 (Â) behaves like ∼ h −2 as for stiffness matrices of standard conforming finite element discretizations of a Poisson problem. For the symmetric Gauss-Seidel preconditioner P SGS , on the finer grid levels the iteration numbers grow approximately like h −1 . For the block preconditioners P A , P D , P B , we observe almost constant iteration numbers for increasing level . For each grid level, the iteration numbers of the block preconditioners are very similar (and even the same for P D and P B ). Note the very small increase in iteration numbers when we change from the exact block preconditioner P A to the inexact ones P D and P B . The third preconditioner, P B , is the only one with computational costs O(N ), N := N 0 + N 1 , with a constant independent of .
We now fix the grid refinement level = 3 and vary the midpoint x 0 = (δ, 2δ, 3δ) penalty stabilization, the condition number κ 2 (Â) has the same order of magnitude. The PCG iteration numbers for all considered preconditioners are (almost) constant for varying δ.
For our analysis to be applicable it is essential that we consider the Nitsche method with stabilization, i.e., β > 0 in (2.6). We also performed numerical experiments with β = 0, where, to ensure positive definiteness of the systems, the Nitsche parameter γ = 100 was chosen larger. The results show that for β = 0 the condition numbers κ 2 (Â) can be extremely large (due to "bad cuts"). However, this does not significantly affect the PCG iteration numbers, which show a similar behavior as for the case with β > 0. These results are consistent with the ones presented in [24]. For higher order finite elements the use of the (ghost penalty) stabilization will be essential.

Poisson fictitious domain problem.
We now consider the Poisson fictitious domain problem in (2.5). Let Ω and Ω 1 be defined as in section 7.1. For the function u : Ω → R, u(x) := (3x 2 1x2 −x 3 2 ) exp(1 − x 2 2 ), the right-hand side f (x) = u(x)(−4 x 2 2 + 18) and boundary data g = u are chosen such that u is a solution of (2.4) on Ω 1 . For discretization the same initial triangulation T 0 of Ω as in section 7.1 is chosen. Applying an adaptive refinement algorithm, where all tetrahedra T ∈ T 0 with meas 3 (T ∩ Ω 1 ) > 0 are marked for regular refinement, we obtain the refined grid T 1 . Repeating this refinement process yields the grids T with refinement levels = 2, . . . , 6 and corresponding grid sizes h = 2 − · 3 4 . We use linear finite elements (k = 1) and construct finite element spaces V h on the respective grids T , = 0, 1, . . . , 6. Table 7.5 reports the numbers N 0 = dim V − 1 (the number of grid points inside the fictitious domain) and N 1 = dim V Γ 1 (the number of grid points on ∂Ω ex 1,h ) for different grid levels. We observe that N 0 and N 1 grow with the expected factors of approximately 8 and 4, respectively. Choosing γ = 10 and β = 0.1, we obtain numerical solutions u h ∈ V h of the discrete problem (2.5), with discretization errors w.r.t. the L 2 and H 1 norm as in Table 7.6. Optimal convergence rates in the L 2 and in the H 1 norm are observed.
We present results for the symmetric Gauss-Seidel preconditioner P SGS and the block Jacobi preconditioners defined in (7.1), where this time B −1 0 denotes one iteration of an algebraic multigrid solver (HYPRE BoomerAMG [20]) applied toÂ 0 and B −1 1 denotes three symmetric Gauss-Seidel iterations. The condition numbers κ 2 (Â) and PCG iteration numbers for different refinement levels are reported in Table 7 Table 7.7: Condition numbers and PCG iteration numbers for different preconditioners and varying grid refinement levels for the fictitious domain problem.
As seen for the interface Poisson problem before, for ≥ 4 the condition number κ 2 (Â) behaves like ∼ h −2 and the iteration numbers for the symmetric Gauss-Seidel preconditioner P SGS grow approximately like h −1 . For the block preconditioners P A , P D , P B , we observe almost constant iteration numbers for increasing level . For all three block preconditioners the number of iterations roughly doubles when going from the coarsest level = 0 to the finest one = 6. The influence of the interface position on condition numbers and PCG iteration numbers shows a similar behavior as for the Poisson interface problem in Table 7.4. We therefore do not report the numbers here.