does not necessarily lead to a coercive bilinear form, as explained in Section
2.1. Therefore, in the results with standard finite element spaces, an
element-wise Nitsche parameter that follows from a local generalized
eigenvalue problem is employed. Note that, as mentioned previously, this local
value of the parameter is not bounded due to the lack of stability of the
method.

We consider hexahedral elements with continuous piecewise trilinear shape
functions. The geometry in the experiments is the complex body depicted in
Figure 9, which is often considered in the literature of immersed finite
element methods; see [74]. The bounding box used to define the background mesh
is the unit cube $[0,1]^{3}$, which is also depicted in Figure 9. We consider
background meshes of hexahedral elements generated and partitioned in parallel
with the $\mathtt{p4est}$ library. The test is performed for three different
local loads (referred to as “load 1”, “load 2” and “load 3”) with $20^{3}$,
$30^{3}$, and $40^{3}$ elements per parallel subdomain, respectively. We
generate the meshes by recursive subdivision of the unit cube. The finest mesh
considered in the examples has $1,073,741,824$ elements (corresponding to
refinement level 10), and it is partitioned into $16,777$ subdomains which are
mapped to the same number of MPI tasks.

Figure 9: View of the geometry and the bounding box considered in the
numerical example. Figure reproduced from [194].

First, we analyze the impact of using aggregated finite-element or standard
finite-element spaces on the number of solver iterations; see Figure 10(a).
This figure conveys that the use of aggregated finite element spaces is
beneficial (and generally indeed essential) to achieve good performance of the
linear solver. For the standard finite element spaces, some results are
missing in Figure 10(a). These correspond to cases in which the linear solver
was not able to provide a converged solution. This shows that standard finite
element spaces are not reliable in immersed discretization methods. In
contrast, aggregated finite element spaces allow one to effectively solve the
underlying linear systems. The results obtained with aggregated finite element
spaces (blue lines) are very close to the expected optimal performance of
multigrid methods (_i.e._ , number of linear solver iterations is
asymptotically independent of the problem size). Figure 10(b) reports the wall
clock time spent in the linear solver step, separated into time spent in the
solver setup and in the solver run. The solver wall clock time increases only
moderately as the problem size increases.

(a) Solver iterations

(b) Wall clock time

Figure 10: Influence of standard or aggregated finite element spaces on the
linear solver behavior. The figure on the top shows the conjugate gradient
(CG) iterations versus problem size for three different loads per processor.
At the bottom, we show the wall clock time for the two main linear solver
phases (setup and run) versus problem size for aggregated finite element
spaces and load 3. Figures reproduced from [194].

###  4.4 Ghost-penalty stabilization

In this section, we present common techniques for the weak stabilization of
immersed finite element methods, often referred to as the _ghost penalty_ and
generally associated with the _Cut Finite Element Method_ (CutFEM) [76, 73,
74]. In the ghost-penalty approach, additional terms are added to the bilinear
form which provide control over the solution on the elements that intersect
the boundary, while the original approximation space $V_{h}$ is unmodified.
These terms are denoted by $s_{h}(v_{h},w_{h})$ and provide a contribution to
the operator norm in the form of $\|v_{h}\|^{2}_{s_{h}}=s_{h}(v_{h},v_{h})$.
The underlying idea is to control the variation of the approximation function
across neighboring elements, such that the solution on small cut elements is
controlled by the solution on interior elements. A standard way to achieve
this is by controlling the jump in the normal derivatives across faces, by
adding appropriately scaled stiffnesses to these jumps. Another approach is to
add a volumetric penalization to the difference between a function itself and
the extension of the function on a neighboring element. In the subsequent
sections, both approaches are discussed.

Figure 11: Set of ghost-penalty faces $\mathcal{F}_{h}^{\mathrm{gh}}$ marked
in red. These faces are interior faces of the active mesh $\mathcal{T}_{h}$
(yellow and blue elements) and belong to at least one cut element in
$\mathcal{T}_{h}^{\textrm{cut}}$ (blue elements).

#### Face-based stabilization

The most common stabilization term for cut elements is the face-based ghost
penalty [73, 74], defined as

$\displaystyle
s_{h}(w_{h},v_{h})=\sum_{F\in\mathcal{F}_{h}^{\mathrm{gh}}}\sum_{j=1}^{p}\tau_{j}h^{2j-1}(\llbracket\partial_{n}^{j}w_{h}\rrbracket,\llbracket\partial_{n}^{j}v_{h}\rrbracket)_{F}$
(4.5)

where $\tau_{j}$ are positive constants, $\mathcal{F}_{h}^{\mathrm{gh}}$ is
the set of all interior faces that are shared by an element in
$\mathcal{T}_{h}^{\textrm{cut}}$ (see Figure 11), $(\cdot,\cdot)_{F}$
indicates the inner product over $F$, and

$\llbracket\partial_{n}v_{h}\rrbracket=\partial_{n_{1}}v_{h,1}+\partial_{n_{2}}v_{h,2}$
(4.6)

is the jump in the normal gradient across the face $F$. In this definition of
the jump on face $F$, $n_{1}$ and $v_{h,1}$, and $n_{2}$ and $v_{h,2}$,
correspond to the exterior unit normal and the function $v_{h}$ on the
elements $T_{1}$ and $T_{2}$ that share face $F$, respectively. Note that,
therefore, $n_{1}=-n_{2}$. With this penalty on the jumps in the normal
derivatives between two neighboring elements, the norm of a function $v_{h}\in
V_{h}$ on one of the elements can be controlled by the norm on the other.
Accordingly, gradients on cut elements are controlled by gradients on full
elements in the interior.

The parameters $\tau_{j}$ in the stabilization terms must be chosen in an
appropriate way. With too small parameters, the effect is negligible, while
with too large parameters, degrees of freedom can be constrained too strongly,
potentially degenerating the accuracy of the approximation [204]. On Dirichlet
boundaries, the required value to provide coercivity of the bilinear form
correlates inversely with the Nitsche parameter. On Neumann boundaries,
coercivity of the bilinear form does not depend on the stabilization term.
Nevertheless, stabilization can still be applied to the Neumann boundary in
order to repair the conditioning problems discussed in Section 3.1, for which
a scaling with $h^{2j+1}$ instead of $h^{2j-1}$ in (4.5) suffices [205].

To elaborate the ghost-penalty approach, let us for simplicity consider the
case of linear elements, _i.e._ , $p=1$. It then holds that

$\displaystyle\|\nabla v_{h}\|^{2}_{T_{1}}\lesssim\|\nabla
v_{h}\|^{2}_{T_{2}}+h\|\llbracket\partial_{n}v_{h}\rrbracket\|^{2}_{F},\qquad\|v_{h}\|^{2}_{T_{1}}\lesssim\|v_{h}\|^{2}_{T_{2}}+h^{3}\|\llbracket\partial_{n}v_{h}\rrbracket\|^{2}_{F},\qquad
v_{h}\in V_{h}$ (4.7)

To derive these estimates, we consider the element-wise restrictions
$v_{h,i}=v_{h}|_{T_{i}}$ ($i\in\\{1,2\\}$), both extended onto $T_{1}\cup
T_{2}$; see Figure 12. The second estimate in (4.7) follows from

$\displaystyle\|v_{h,1}\|^{2}_{T_{1}}$
$\displaystyle\lesssim\|v_{h,1}-v_{h,2}\|^{2}_{T_{1}}+\|v_{h,2}\|^{2}_{T_{1}}\lesssim\|(\boldsymbol{\xi}-\boldsymbol{\xi}_{F})\cdot\nabla(v_{h,1}-v_{h,2})\|^{2}_{T_{1}}+\|v_{h,2}\|^{2}_{T_{2}}$
(4.8) $\displaystyle\lesssim
h^{3}\|\partial_{n_{i}}(v_{h,1}-v_{h,2})\|^{2}_{F}+\|v_{h,2}\|^{2}_{T_{2}}=h^{3}\|\llbracket\partial_{n}v\rrbracket\|^{2}_{F}+\|v_{h,2}\|^{2}_{T_{2}}$

with $\partial_{n_{i}}$ indicating either of the normal derivatives. Here we
first added and subtracted $v_{h,2}$ and used the triangle inequality. We then
used the Taylor expansion
$v_{h,1}-v_{h,2}=(\boldsymbol{\xi}-\boldsymbol{\xi}_{F})\cdot\nabla(v_{h,1}-v_{h,2})$
at $\boldsymbol{\xi}$ with point $\boldsymbol{\xi}_{F}\in F$ the projection of
$\boldsymbol{\xi}$ along the normal direction, which holds since
$v_{h,1}-v_{h,2}=0$ on $F$ by virtue of the continuity of $v_{h}$.
Additionally, we used the estimate
$\|v_{h,2}\|_{T_{1}}\lesssim\|v_{h,2}\|_{T_{2}}$ which holds since $v_{h,2}$
is polynomial. The first estimate follows in the same way, but is simpler
since the gradients are piecewise constant, leading to an $h$ scaling factor
instead of the $h^{3}$ factor.

$T_{1}$$T_{2}$$v_{h,1}\,(=v)$$v_{h,1}-v_{h,2}$$v_{h,2}\,(\text{extension})$$v_{h,2}\,(=v)$$v_{h,1}\,(\text{extension})$$F$
Figure 12: Illustration of the extension and difference of $v_{h,1}$ and
$v_{h,2}$ in (4.8).

Using the pairwise bounds in (4.7) – either directly or through a chain of
neighboring elements – the degrees of freedom on cut elements in
$\mathcal{T}_{h}^{\textrm{cut}}$ can be controlled by the degrees of freedom
on interior elements $\mathcal{T}_{h}^{\textrm{in}}$ via the ghost penalty

$\|\nabla v_{h}\|^{2}_{\Omega_{h}}\lesssim\|\nabla
v_{h}\|^{2}_{\Omega}+s_{h}(v_{h},v_{h}),\qquad\|v_{h}\|^{2}_{\Omega_{h}}\lesssim\|v_{h}\|^{2}_{\Omega}+h^{2}s_{h}(v_{h},v_{h}),\qquad
v_{h}\in V_{h}$ (4.9)

To prove the first bound we proceed as

$\displaystyle\|\nabla v_{h}\|^{2}_{\Omega_{h}}$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}}\|\nabla v_{h}\|^{2}_{T}$ (4.10)
$\displaystyle=\sum_{T\in\mathcal{T}_{h}^{\textrm{in}}}\|\nabla
v_{h}\|^{2}_{T}+\sum_{T\in\mathcal{T}_{h}^{\textrm{cut}}}\|\nabla
v_{h}\|^{2}_{T}$
$\displaystyle\lesssim\sum_{T\in\mathcal{T}_{h}^{\textrm{in}}}\|\nabla
v_{h}\|^{2}_{T}+\Big{(}\sum_{T\in\mathcal{T}_{h}^{\textrm{in}}}\|\nabla
v_{h}\|^{2}_{T}+s_{h}(v_{h},v_{h})\Big{)}$ $\displaystyle\lesssim\|\nabla
v_{h}\|^{2}_{\Omega}+s_{h}(v_{h},v_{h})$

The second bound in (4.9) can be established similarly. Note that the
constants in these estimates will depend on the lengths of the paths
connecting cut elements to interior elements, which is related to the
properties and the resolution of the boundary, as well as the properties of
the mesh. For smooth boundaries and locally quasi-uniform meshes with
sufficiently small mesh sizes, one can show that the length of these paths is
uniformly bounded.

Conversely, the stabilization term can be bounded in terms of the $L^{2}$-norm
of the gradient of the approximation on the active domain $\Omega_{h}$

$\displaystyle s_{h}(v_{h},v_{h})\lesssim\|\nabla v_{h}\|_{\Omega_{h}}\qquad
v_{h}\in V_{h}$ (4.11)

which follows from

$\displaystyle
h^{2j-1}\|\llbracket\partial^{j}_{n}v_{h}\rrbracket\|^{2}_{F}\lesssim\sum_{i=1}^{2}h^{2j-1}\|\partial^{j}_{n}v_{h,i}\|^{2}_{F}\lesssim\sum_{i=1}^{2}h^{2j-2}\|\nabla^{j}v_{h,i}\|^{2}_{T_{i}}\lesssim\sum_{i=1}^{2}\|\nabla
v_{h,i}\|^{2}_{T_{i}}=\|\nabla v_{h}\|_{T_{1}\cup T_{2}}^{2}$ (4.12)

where we first employed the triangle inequality, then a trace inverse estimate
to pass from the face $F$ to the neighboring elements $T_{1}$ and $T_{2}$
[97], and finally another inverse estimate to arrive at a bound in terms of
$\nabla v_{h}$ [93]. Using the same inverse estimate it follows that

$h^{2}\|v_{h}\|^{2}_{s_{h}}\lesssim h^{2}\|\nabla
v_{h}\|^{2}_{\Omega_{h}}\lesssim\|v_{h}\|^{2}_{\Omega_{h}}$ (4.13)

The estimates (4.9), (4.11), and (4.13) are central to the coercivity and
continuity of the (stabilized) bilinear form with a (global) Nitsche parameter
that scales inversely with the mesh size of the _background_ grid. The
estimates are also essential for the condition number bounds derived in
Section 4.5. The principal notion of the ghost-penalty method is encapsulated
in the bounds in (4.9), which convey that the ghost penalty provides control
over the solution on the entire active domain including the fictitious part.
The bounds (4.11) and (4.13) in turn impart that the ghost-penalty term is
suitably well behaved.

#### Element- and patch-based stabilization forms

An alternative to the face-based penalty is a penalty that is based on
volumetric integrals. This is convenient for higher-order polynomials, since
it avoids the computation of higher-order derivatives at element interfaces,
which is a non-standard operation in many finite element codes. An overview
and analysis of different types of such ghost-penalty formulations is
presented in [206].

A stabilization term based on element integrals can be defined according to
any of the two forms

$\displaystyle
s_{h,0}^{\mathrm{elem}}(w_{h},v_{h})=\sum_{F\in\mathcal{F}_{h}^{\mathrm{gh}}}\tau_{0}h^{-2}([w_{h}],[v_{h}])_{\mathcal{T}_{h}(F)},\quad\textrm{or}\quad
s_{h,1}^{\mathrm{elem}}(w_{h},v_{h})=\sum_{F\in\mathcal{F}_{h}^{\mathrm{gh}}}\tau_{1}([\nabla
w_{h}],[\nabla v_{h}])_{\mathcal{T}_{h}(F)}$ (4.14)

where for each face $F\in\mathcal{F}_{h}^{\mathrm{gh}}$, $\mathcal{T}_{h}(F)$
is the pair of elements that share $F$ and

$\displaystyle[v_{h}]=v_{h,1}-v_{h,2},\qquad[\nabla v_{h}]=\nabla
v_{h,1}-\nabla v_{h,2}$ (4.15)

on $T\in\mathcal{T}_{h}(F)=\\{T_{1},T_{2}\\}$ with polynomials $v_{h,1}$ and
$v_{h,2}$ extended in the canonical way; see Figure 12. Note that, opposed to
the existing volumetric integrals in the weak formulation introduced in
Section 2.1, the integrals in the penalization are taken over the entire
background elements including the fictitious parts and are not restricted to
the physical domain $\Omega$. The form in (4.14) is presented in [207, 208]
and is obtained by example formulation 2 in [206, §2.4] in case the element
pairs are chosen based on the faces in $\mathcal{F}_{h}^{\mathrm{gh}}$.

We may also consider stabilization based on patches or macro elements, which
is the form that was presented with the introduction of the ghost penalty in
[76] and, for a certain choice of patches, is also considered in [208]. With
this form of stabilization, each cut element is associated with a patch of
neighboring elements that contains at least one element in the interior. We
may then penalize the difference between the finite element function and a
global polynomial on the patch. More precisely, let $M_{T}$ be a patch of
elements containing $T$, then we can define

$\displaystyle s_{h,0}^{\mathrm{patch}}(w_{h},v_{h})$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}^{\textrm{cut}}}\tau_{0}h^{-2}(w_{h}-P_{p}w_{h},v_{h}-P_{p}v_{h})_{M_{T}}$
(4.16) $\displaystyle s_{h,1}^{\mathrm{patch}}(w_{h},v_{h})$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}^{\textrm{cut}}}\tau_{1}(\nabla(w_{h}-P_{p}w_{h}),\nabla(v_{h}-P_{p}v_{h}))_{M_{T}}$

where $P_{p}:V_{h}|_{M_{T}}\rightarrow\mathcal{P}_{p}(M_{T})$ is the
$L^{2}$-projector or the Ritz projector, respectively. Note that by strictly
following this formulation, patches containing multiple cut elements will be
contained in the summation multiple times. In practice, every patch that
contains cut elements is included only once, or the boundary is simply covered
by a set of non-overlapping patches such as in [209].

Ghost-penalty stabilization can also be combined with the concept of
aggregation to design weakly aggregated schemes, as proposed in [206, 78]. In
these schemes, one makes use of the standard finite element space $V_{h}$ on
the active mesh, and instead of strongly imposing constraints via an extension
operator of interior degrees of freedom, the difference between $V_{h}$ and
$V_{h}^{\textrm{ag}}$ is weakly penalized. Recall the discrete extension
operator $\mathcal{E}_{h}^{\mathrm{ag}}:V_{h}^{\mathrm{in}}\rightarrow
V_{h}^{\mathrm{ag}}\in V_{h}$ defined in Section 4.3 to define a mapping
$P^{\textrm{ag}}_{h}:V_{h}\ni
v_{h}\mapsto\mathcal{E}_{h}^{\mathrm{ag}}(v_{h}|_{\mathcal{T}_{h}^{\textrm{in}}})\in
V_{h}^{\textrm{ag}}$. The difference $v_{h}-P^{\textrm{ag}}_{h}v_{h}$ can then
be penalized similarly as in (4.16)

$\displaystyle s_{h}^{\mathrm{ag}}(w_{h},v_{h})$
$\displaystyle=\sum_{\mathcal{T}_{h}}\tau_{0}h^{-2}\left(w_{h}-P^{\textrm{ag}}_{h}w_{h},v_{h}-P^{\textrm{ag}}_{h}v_{h}\right)_{T}$
(4.17) $\displaystyle s_{h}^{\mathrm{ag}}(w_{h},v_{h})$
$\displaystyle=\sum_{\mathcal{T}_{h}}\tau_{1}\left(\nabla(w_{h}-P^{\textrm{ag}}_{h}w_{h}),\nabla(v_{h}-P^{\textrm{ag}}_{h}v_{h})\right)_{T}$

These terms can be computed element-wise (in a similar manner as described in
Section 4.3) and are similar to local-projection stabilization techniques;
see, _e.g._ , [210, 211]. We note that the (strong) aggregated finite element
method of Section 4.3 is formally recovered in the limit $\tau_{j}\to\infty$.
Thus, the stabilization method is also accurate for large values of
$\tau_{j}$. This last statement does not hold for other forms of the ghost
penalty, in which a large parameter can constrain too many degrees of freedom,
which inhibits the accuracy [204] and is referred to as locking in [206].
Detailed analyses and experimentation with this type of stabilization can be
found in [78] and [206].

Stabilization techniques based on volumetric integrals have also been
developed for multimesh discretizations [212, 213, 214] and for problems on
manifolds [75]. Additionally, it is noteworthy that with the $H^{1}$-product
and with the projector $P_{p}=0$ or $P^{\textrm{ag}}_{h}=0$, the fictitious
domain stiffness that is commonly employed in the finite cell method is
recovered [215]. Fictitious domain stiffness requires $\tau_{1}\ll 1$,
however, as this stabilization term is not consistent with the original
problem. Therefore, with a fictitious domain stiffness, stability and
conditioning need to be balanced with consistency.

#### Implementation aspects

To implement any of the aforementioned stabilization forms, one must be able
to find face neighbors of given elements and assemble contributions from pairs
of elements sharing a face. For higher-order elements, the face-based ghost
penalty requires computation of higher-order derivatives. The element-based
formulation may then be preferable. All the stabilization terms lead to
additional coupling and, as a consequence, fill-in in the stiffness matrix
[216, 52]. Systems stabilized by the ghost penalty have similar condition
number bounds as boundary-fitted formulations, and therefore larger systems
(_i.e._ , finer grids) yield larger condition numbers and require more
iterations in iterative solvers. This can be overcome by combining the ghost
penalty with a multigrid solver; see the discussion on this in Section 3.2.

#### Extension to other problems

The ghost penalty can also be applied to stabilize time-dependent problems,
for which generally the mass form needs to be stabilized as well. In that case
the scaling is chosen in such a way that the stabilization terms scale with
the mesh parameter in the same way as the mass matrix; see (4.13). Similar to
Schwarz preconditioning and aggregation, the ghost penalty can also be applied
to provide inf-sup stable discretization of mixed formulations such as that of
the Stokes problem; see, _e.g._ , [21, 28, 217, 52]. Furthermore, the ghost
penalty can also be applied to stabilize discretizations of codimensional
manifolds, such as partial differential equations posed on non-planar surfaces
discretized with three-dimensional meshes; see, _e.g._ , [218, 219, 220, 75].

###  4.5 A unified analysis of aggregation and stabilization

In this section, we present a framework for the analysis of stable immersed
finite element methods. We consider aggregated and ghost-penalty-stabilized
techniques at the same time, to identify the common parts and the main
differences. It is recalled from Section 4.1 that when Nitsche’s method is
applied without any stability-enhancing technique, the bilinear form is only
coercive with an _a priori_ unbounded Nitsche parameter $\beta$. In that case,
most of the properties and estimates presented in this section do not apply,
and error bounds cannot be provided. Remark 4.2 at the end of this section
discusses which specific analyses still hold and which no longer hold without
any form of stabilization.

#### Methods and formulations

We recall from (2.6) and (2.7) that the standard Nitsche finite element method
takes the form

$\displaystyle\begin{cases}\text{find }u_{h}\in V_{h}\text{ such that:}\\\
a_{h}(u_{h},v_{h})=l_{h}(v_{h})\qquad\forall v_{h}\in V_{h}\end{cases}$ (4.18)

with bilinear form $a_{h}(\cdot,\cdot)$ according to

$\displaystyle a_{h}(w_{h},v_{h})=(\nabla w_{h},\nabla
v_{h})_{\Omega}-(w_{h},\partial_{n}v_{h})_{\partial\Omega_{D}}-(\partial_{n}w_{h},v_{h})_{\partial\Omega_{D}}+(\beta
w_{h},v_{h})_{\partial\Omega_{D}}$ (4.19)

We note that this is consistent with the strong formulation in (2.1), _i.e._ ,
the exact solution $u$ to the corresponding continuous problem satisfies

$\displaystyle a_{h}(u,v_{h})=l_{h}(v_{h})\qquad v_{h}\in V_{h}$ (4.20)

provided that $u$ is regular enough, for instance if $u\in H^{2}(\Omega)$. To
accommodate both approaches in the same analysis we consider the formulation

$\displaystyle\begin{cases}\text{find }u_{h}\in V_{h,\bigstar}\text{ such
that:}\\\ a_{h,\bigstar}(u_{h},v_{h})=l_{h}(v_{h})\qquad\forall v_{h}\in
V_{h,\bigstar}\end{cases}$ (4.21)

where the finite dimensional function space is defined by

$\displaystyle V_{h,\bigstar}=\begin{cases}V_{h}^{\textrm{ag}}\subset
V_{h}&\text{(aggregation)}\\\ V_{h}&\text{(ghost-penalty
stabilization)}\end{cases}$ (4.22)

and the bilinear form is defined by

$a_{h,\bigstar}(w_{h},v_{h})=a_{\Omega}(w_{h},v_{h})-(w_{h},\partial_{n}v_{h})_{\partial\Omega_{D}}-(\partial_{n}w_{h},v_{h})_{\partial\Omega_{D}}+(\beta
w_{h},v_{h})_{\partial\Omega_{D}}$ (4.23)

with

$\displaystyle a_{\Omega}(w_{h},v_{h})=\begin{cases}(\nabla w_{h},\nabla
v_{h})_{\Omega}&\text{(aggregation)}\\\ (\nabla w_{h},\nabla
v_{h})_{\Omega}+s_{h}(w_{h},v_{h})&\text{(ghost-penalty
stabilization)}\end{cases}$ (4.24)

For the ghost-penalty method, the stabilization form is included in the
definition of $a_{\Omega}(\cdot,\cdot)$ because the stabilization form
$s_{h}(\cdot,\cdot)$ is naturally paired with the bulk or volumetric term
$(\nabla\cdot,\nabla\cdot)_{\Omega}$ to enable bounds of the form (4.7) and
(4.9).

#### Properties

For the analysis, we start from the following intrinsic properties of
aggregation methods and ghost-penalty formulations:

  * •

There is a constant (hidden in the binary operator $\lesssim$) such that

$\displaystyle\|\nabla v_{h}\|_{\Omega_{h}}\lesssim
a_{\Omega}(v_{h},v_{h})=\|v_{h}\|_{a_{\Omega}}\qquad v_{h}\in V_{h,\bigstar}$
(4.25)

where as before $\Omega_{h}=\cup_{T\in\mathcal{T}_{h}}T$ is the union of all
the active elements, such that $\Omega_{h}\supset\Omega$. See, respectively,
(4.3) and (4.10) for these estimates in the context of aggregated and weakly
stabilized methods. This property forms the basis of the effectiveness of
weakly stabilized and aggregated methods. Essentially, it indicates that the
solution on the extended domain $\Omega_{h}$, including the fictitious parts
outside $\Omega$ of active elements, is controlled by
$a_{\Omega}(\cdot,\cdot)$. Consequently, this precludes the very small
eigenvalues on small cut elements discussed in Section 3.1, and it provides
control of the normal gradients which enables bounding the Nitsche parameter
$\beta$ from above as discussed in Section 4.1 and 4.2. With aggregation this
is achieved by directly precluding basis functions with a very small support
in the physical domain as discussed in Section 4.3, while with the ghost
penalty such functions are given a contribution in $a_{\Omega}(\cdot,\cdot)$
by $s_{h}(\cdot,\cdot)$ as discussed in Section 4.4.

  * •

There is an interpolation operator $\pi_{h}:H^{s}(\Omega)\rightarrow
V_{h,\bigstar}$, $s\geq 2$, such that

$\displaystyle\|\nabla^{m}(v-\pi_{h}v)\|_{\Omega_{h}}\lesssim
h^{s-m}\|v\|_{H^{s}(\Omega)}\qquad 0\leq m\leq s\leq p+1$ (4.26)

we use a continuous extension operator $\mathcal{E}:v\in
H^{s}(\Omega)\rightarrow H^{s}(\mathbb{R}^{d})$ to extend $v$ outside of
$\Omega$. For ghost-penalty-stabilized methods, the interpolation operator is
generally constructed by first extending the function $v$ outside of $\Omega$
and then applying a standard interpolation operator, for instance the Clément
operator [221]. For the aggregated finite element space, a specific
construction of the extension is used to establish the interpolation estimate
[222, 192].

  * •

The stabilization form $s_{h}(\cdot,\cdot)$ satisfies

$\displaystyle\|\pi_{h}v\|_{s_{h}}^{2}=s_{h}(\pi_{h}v,\pi_{h}v)\lesssim
h^{s-1}\|v\|_{H^{s}(\Omega)}\qquad 1\leq s\leq p+1$ (4.27)

This is a consistency assumption on the stabilization form
$s_{h}(\cdot,\cdot)$, formulated in such a way that the stabilization form
only acts directly on functions in the finite element space, and acts on other
functions $v$ indirectly via their projection $\pi_{h}{}v$. Note that the
stabilization form $s_{h}(\cdot,\cdot)$ may for instance involve higher-order
derivatives, which are not generally well defined for the exact solution or
the solution to a dual problem used in an error analysis. This complication is
avoided by formulating the consistency condition (4.27) in terms of
projections onto the finite element space.

To provide an example we verify this consistency assumption (4.27) for the
face-based stabilization form (4.5) in the case of linear elements. Consider a
face $F\in\mathcal{F}_{h}$ shared by the two elements $T_{1}$ and $T_{2}$.
Then for any linear polynomial $w\in\mathcal{P}_{1}(T_{1}\cup T_{2})$ we have

$h\|[\partial_{n}\pi_{h}v]\|^{2}_{F}=h\|[\partial_{n}(\pi_{h}v-w)]\|^{2}_{F}\lesssim\sum_{i=1}^{2}\|\nabla(\pi_{h}v-w)\|^{2}_{T_{i}}$
(4.28)

where we used a trace inverse estimate to pass from the face to the elements
[97]. Using the Bramble-Hilbert lemma (see, _e.g._ , [223]) we can choose $w$
in such a way that

$\displaystyle\sum_{i=1}^{2}\|\nabla(\pi_{h}v-w)\|^{2}_{T_{i}}\leq\sum_{i=1}^{2}\|\nabla(\pi_{h}v-v)\|^{2}_{T_{i}}+\|\nabla(v-w)\|^{2}_{T_{i}}\lesssim\sum_{i=1}^{2}h^{2}\|v\|^{2}_{H^{2}(S(T_{i}))}$
(4.29)

where $S(T_{i})$ is the union of the set of elements in $\mathcal{T}_{h}$ that
share a node with $T$.

Element aggregation and ghost-penalty stabilization are conceptually different
in the manner in which they approach the conditions (4.25)–(4.27). Aggregation
methods aim to satisfy condition (4.25) by restricting the approximation to a
suitable subspace $V_{h}^{\textrm{ag}}\subset V_{h}$. On the other hand,
condition (4.26) demands optimal approximation properties of the (sequence of)
aggregated approximation spaces $V_{h}^{\textrm{ag}}$ and, hence, the
aggregated approximation spaces cannot be too small. Because the stabilization
operator $s_{h}(\cdot,\cdot)$ vanishes in aggregation methods, condition
(4.27) is trivially satisfied. Ghost-penalty stabilization requires that the
stabilization term $s_{h}(\cdot,\cdot)$ is strong enough to comply with
condition (4.25). On the other hand, condition (4.27) insists that the
stabilization cannot be too strong, as this will affect the solution. Because
ghost-penalty stabilization imposes no auxiliary conditions on the
approximation spaces, condition (4.26) is trivially satisfied.

#### Coercivity and Continuity

We define the following energy norms

$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|v|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}^{2}$
$\displaystyle=\|\nabla
v\|^{2}_{\Omega}+\|h^{-\frac{1}{2}}v\|^{2}_{\partial\Omega_{D}}+\|h^{\frac{1}{2}}\partial_{n}v\|^{2}_{\partial\Omega_{D}}$
(4.30)
$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}^{2}$
$\displaystyle=a_{\Omega}(v_{h},v_{h})+\|h^{-\frac{1}{2}}v_{h}\|^{2}_{\partial\Omega_{D}}+\|h^{\frac{1}{2}}\partial_{n}v_{h}\|^{2}_{\partial\Omega_{D}}$
(4.31)

and note that

$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h,\bigstar}=\begin{cases}|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h}&\text{aggregation}\\\
|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h}+\|v_{h}\|_{s_{h}}^{2}&\text{weak
stabilization}\end{cases}$ (4.32)

Note that these norms are different from the $\beta$-norm that was applied in
[1] and presented in (4.1) in Section 4.1. First,
$|\mspace{-1.0mu}|\mspace{-1.0mu}|\cdot|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}$
employs $h^{-1}$, which is uniformly bounded, instead of $\beta$, which is not
bounded in case no stabilization technique is applied. Second, in case of
ghost-penalty stabilization,
$|\mspace{-1.0mu}|\mspace{-1.0mu}|\cdot|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}$
also includes a ghost-penalty contribution. Consequently, whereas the
$\beta$-norm admits a natural analysis of unstabilized forms, the
$|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}$-norm
naturally accompanies weakly stabilized and aggregated formulations.

######  Lemma 4.1.

The form $a_{h,\bigstar}(\cdot,\cdot)$ is continuous on $V_{h,\bigstar}\oplus
H^{2}(\Omega)$

$\displaystyle
a_{h,\bigstar}(w,v)\leq|\mspace{-1.0mu}|\mspace{-1.0mu}|w|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}|\mspace{-1.0mu}|\mspace{-1.0mu}|v|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}\qquad
w,v\in V_{h,\bigstar}\oplus H^{2}(\Omega)$ (4.33)

and, if $\beta$ is chosen large enough, $a_{h,\bigstar}(\cdot,\cdot)$ is
coercive on $V_{h,\bigstar}$

$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}^{2}\lesssim
a_{h,\bigstar}(v_{h},v_{h})\qquad v_{h}\in V_{h,\bigstar}$ (4.34)

  * Proof.Under the notion that with aggregation or weak stabilization $\beta\lesssim h^{-1}$, the continuity follows directly from the Cauchy-Schwarz inequality together with the observation that $\partial_{n}v\in L^{2}(\partial\Omega_{D})$ for $v\in V_{h,\bigstar}\oplus H^{2}(\Omega)$ and therefore $\|\partial_{n}v\|_{\partial\Omega_{D}}$ is well defined. For the coercivity we use the element-wise trace inverse inequality [97]

$h\|\partial_{n}v\|^{2}_{T\cap\partial\Omega_{D}}\lesssim\|\nabla
v\|^{2}_{T}\qquad v\in\mathcal{P}_{p}(T)$ (4.35)

and (4.25) to conclude that

$\displaystyle\|h^{\frac{1}{2}}\partial_{n}v_{h}\|^{2}_{\partial\Omega_{D}}\lesssim\sum_{T\in\mathcal{T}_{h}^{\textrm{cut}}}\|\nabla
v_{h}\|^{2}_{T}\lesssim\sum_{T\in\mathcal{T}_{h}}\|\nabla
v_{h}\|^{2}_{T}\lesssim\|\nabla
v_{h}\|^{2}_{\Omega_{h}}\lesssim\|v_{h}\|^{2}_{a_{\Omega}}\qquad v_{h}\in
V_{h,\bigstar}$ (4.36)

where $\mathcal{T}_{h}^{\textrm{cut}}$ is the set of elements that intersect
the boundary. Next, for any $\delta>0$, we have

$\displaystyle 2(\partial_{n}v,v)_{\partial\Omega_{D}}\leq
2\|\partial_{n}v\|_{\partial\Omega_{D}}\|v\|_{\partial\Omega_{D}}\leq\|(h\delta)^{\frac{1}{2}}\partial_{n}v\|^{2}_{\partial\Omega_{D}}+\|(h\delta)^{-\frac{1}{2}}v\|^{2}_{\partial\Omega_{D}}$
(4.37)

which leads to

$\displaystyle a_{h,\bigstar}(v_{h},v_{h})$
$\displaystyle=\|v_{h}\|^{2}_{a_{\Omega}}-2(\partial_{n}v_{h},v_{h})_{\partial\Omega_{D}}+\|\beta^{\frac{1}{2}}v_{h}\|^{2}_{\partial\Omega_{D}}$
(4.38)
$\displaystyle\geq\|v_{h}\|^{2}_{a_{\Omega}}-\|(h\delta)^{\frac{1}{2}}\partial_{n}v_{h}\|^{2}_{\partial\Omega_{D}}-\|(h\delta)^{-\frac{1}{2}}v_{h}\|^{2}_{\partial\Omega_{D}}+\|\beta^{\frac{1}{2}}v\|^{2}_{\partial\Omega_{D}}$
$\displaystyle\geq(1-\delta
C)\|v_{h}\|^{2}_{a_{\Omega}}+\|(\beta-(h\delta)^{-1})^{\frac{1}{2}}v_{h}\|^{2}_{\partial\Omega_{D}}\qquad
v_{h}\in V_{h,\bigstar}$

where $C$ is the hidden constant between the first and final expressions in
(4.36). Taking $\delta$ small enough to guarantee that $\delta C\lesssim 1/2$
and $\beta$ large enough such that $\beta-(h\delta)^{-1}>h^{-1}/2$ on each
element in $\mathcal{T}_{h}^{\textrm{cut}}$, the coercivity follows. ∎

Assuming that $l_{h}(v_{h})$ is continuous, _i.e._ ,
$l_{h}(v_{h})\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}$
for $v_{h}\in V_{h,\bigstar}$, we may now apply the Lax–Milgram lemma to show
that there is a unique and stable solution $u_{h}\in V_{h,\bigstar}$ to
(4.21).

#### Error Estimates

Conditions (4.25)–(4.27) ensure that the corresponding immersed finite-element
approximation according to (4.18) has optimal approximation properties in the
energy norm.

######  Theorem 4.1.

Consider the immersed finite-element approximation (4.18). Assume that
conditions (4.25)–(4.27) hold. The approximate solution $u_{h}$ satisfies the
energy-norm error estimate

$|\mspace{-1.0mu}|\mspace{-1.0mu}|u-u_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}\lesssim
h^{p}\|u\|_{H^{p+1}(\Omega)}$ (4.39)

  * Proof.We start by splitting the error in two parts using the interpolation operator

$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|u-u_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|u-\pi_{h}u|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}+|\mspace{-1.0mu}|\mspace{-1.0mu}|\pi_{h}u-u_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$
(4.40)

We first consider the first term and show that

$|\mspace{-1.0mu}|\mspace{-1.0mu}|u-\pi_{h}u|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}\lesssim
h^{p}\|u\|_{H^{p+1}(\Omega)}$ (4.41)

To establish (4.41), we note that by virtue of a standard trace inequality
(see, _e.g._ , [223])

$\|v\|^{2}_{T\cap\partial\Omega_{D}}\lesssim h^{-1}\|v\|^{2}_{T}+h\|\nabla
v\|^{2}_{T}\qquad v\in H^{1}(T)$ (4.42)

it holds that

$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|v|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h}$
$\displaystyle=\|\nabla
v\|^{2}_{\Omega}+\|h^{-\frac{1}{2}}v\|^{2}_{\partial\Omega_{D}}+\|h^{\frac{1}{2}}\partial_{n}v\|^{2}_{\partial\Omega_{D}}$
(4.43) $\displaystyle\lesssim\|\nabla
v\|^{2}_{\Omega_{h}}+\|h^{-1}v\|^{2}_{\Omega_{h}}+\|h\nabla^{2}v\|^{2}_{\Omega_{h}}$

for $v\in V_{h,\bigstar}+H^{2}(\Omega)$. Setting $v=u-\pi_{h}u$ and using
(4.26), we obtain (4.41). Equation (4.43) applies the definition of
$|\mspace{-1.0mu}|\mspace{-1.0mu}|\cdot|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h}$
given in (4.30).

Considering the second term in (4.40), we note that the definition of
$|\mspace{-1.0mu}|\mspace{-1.0mu}|\cdot|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}$
and the coercivity property (4.34) imply

$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|\pi_{h}u-u_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|\pi_{h}u-u_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}\lesssim\sup_{w_{h}\in
V_{h,\bigstar}\setminus\\{0\\}}\frac{a_{h,\bigstar}(\pi_{h}u-u_{h},w_{h})}{|\mspace{-1.0mu}|\mspace{-1.0mu}|w_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}}$
(4.44)

In the case of the aggregation method, we have
$a_{h,\bigstar}(\cdot,\cdot)=a_{h}(\cdot,\cdot)$ and
$|\mspace{-1.0mu}|\mspace{-1.0mu}|\cdot|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}=|\mspace{-1.0mu}|\mspace{-1.0mu}|\cdot|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$
and it immediately follows that

$\displaystyle a_{h,\bigstar}(\pi_{h}u-u_{h},w_{h})$
$\displaystyle=a_{h}(\pi_{h}u-u_{h},w_{h})$ (4.45)
$\displaystyle=a_{h}(\pi_{h}u,w_{h})-l_{h}(w_{h})$
$\displaystyle=a_{h}(\pi_{h}u-u,w_{h})$
$\displaystyle\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|\pi_{h}u-u|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|w_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$
$\displaystyle\lesssim
h^{p}\|u\|_{H^{p+1}(\Omega)}|\mspace{-1.0mu}|\mspace{-1.0mu}|w_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$

where we used, consecutively, the consistency relation (4.20) to replace
$l_{h}(w_{h})$ by $a_{h}(u,w_{h})$, continuity of $a_{h}$ according to (4.33),
and the interpolation error estimate (4.41). In the case of ghost-penalty
stabilization, it holds that

$\displaystyle a_{h,\bigstar}(\pi_{h}u-u_{h},w_{h})$
$\displaystyle=a_{h,\bigstar}(\pi_{h}u,w_{h})-l_{h}(w_{h})$ (4.46)
$\displaystyle=a_{h}(\pi_{h}u,w_{h})-l_{h}(w_{h})+s_{h}(\pi_{h}u,w_{h})$
$\displaystyle=a_{h}(\pi_{h}u-u,w_{h})+s_{h}(\pi_{h}u,w_{h})$
$\displaystyle\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|\pi_{h}u-u|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|w_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}+\|\pi_{h}u\|_{s_{h}}\|w_{h}\|_{s_{h}}$
$\displaystyle\lesssim\Big{(}|\mspace{-1.0mu}|\mspace{-1.0mu}|\pi_{h}u-u|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h}+\|\pi_{h}u\|^{2}_{s_{h}}\Big{)}^{1/2}|\mspace{-1.0mu}|\mspace{-1.0mu}|w_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}$
$\displaystyle\lesssim
h^{p}\|u\|_{H^{p+1}(\Omega)}|\mspace{-1.0mu}|\mspace{-1.0mu}|w_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}$

where we used the interpolation error estimate (4.41) and the weak consistency
condition (4.27) on $s_{h}(\cdot,\cdot)$. From (4.44)–(LABEL:eq:energy-f) it
follows that the second term in (4.40) is also bounded in accordance with
(4.39). ∎

#### Condition number estimate

Under certain (non-restrictive) conditions, both aggregation and ghost-penalty
stabilization engender stiffness matrices that exhibit similar condition
numbers as conventional boundary-fitted finite element methods. To analyze the
condition numbers corresponding to both methods, we consider the basis
$\\{\phi_{i}\\}_{i=1}^{N}$ in $V_{h,\bigstar}$, so that any
$v_{h}\in{}V_{h,\bigstar}$ can be expanded as

$\displaystyle v_{h}=\sum_{i=1}^{N}y_{i}\phi_{i}$ (4.47)

where $\mathbf{y}\in\mathbb{R}^{N}$ is the coefficient vector. It is to be
noted that with ghost-penalty stabilization the basis $V_{h,\bigstar}=V_{h}$
is a standard nodal basis, while with aggregation
$V_{h,\bigstar}=V_{h}^{\textrm{ag}}$ is the basis that is modified at the
boundary in such a way that the support of each basis function has a
sufficiently large intersection with $\Omega$, more precisely
$|\text{supp}(\phi_{i})\cap\Omega|\sim h^{d}$. The stiffness matrix
$\mathbf{A}\in\mathbb{R}^{N\times N}$ is defined as

$\mathbf{z}^{T}\mathbf{A}\mathbf{y}=a_{h,\bigstar}(w_{h},v_{h})\qquad
w_{h},v_{h}\in V_{h,\bigstar}$ (4.48)

and for the considered symmetric bilinear forms, the condition number of the
stiffness matrix $\mathbf{A}$ corresponds to

$\kappa\left(\mathbf{A}\right)=\frac{\lambda_{\max}}{\lambda_{\min}}$ (4.49)

where $\lambda_{\max}$ and $\lambda_{\min}$ are the maximum and minimum
eigenvalues of $\mathbf{A}$, respectively. To estimate the condition number we
start from the following basic bounds, which hold by construction for both
aggregated and weakly stabilized methods:

  * •

There is a constant (hidden in the binary relation $\lesssim$) such that

$\displaystyle|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}\lesssim
h^{-1}\|v_{h}\|_{\Omega_{h}}\qquad v_{h}\in V_{h,\bigstar}$ (4.50)

Both for weakly stabilized and aggregated methods, this bound follows from
standard inverse inequalities; see (4.11) for the inverse estimate of the
ghost-penalty form.

  * •

There is a (hidden) constant such that

$\displaystyle\|v_{h}\|_{\Omega_{h}}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h,\bigstar}\qquad
v_{h}\in V_{h,\bigstar}$ (4.51)

This is a Poincaré inequality, requiring aggregation of the approximation
space or ghost-penalty stabilization of the bilinear form to gain control on
the approximate solution outside $\Omega$.

  * •

There are constants such that the following equivalence holds:

$\displaystyle h^{d}\|\mathbf{y}\|^{2}_{2}\sim\|v_{h}\|^{2}_{\Omega_{h}}\qquad
v_{h}\in V_{h,\bigstar}$ (4.52)

For weakly stabilized methods we have $V_{h,\bigstar}=V_{h}$ and this is a
standard estimate that follows from the quasi-uniformity of the mesh and the
use of a nodal basis. For the aggregated finite element space
$V_{h,\bigstar}=V_{h}^{\textrm{ag}}$, this estimate holds by virtue of the
construction of the extension operator.

######  Theorem 4.2.

Consider the stiffness matrix in (4.48). Assume that the bilinear form
$a_{h,\bigstar}(\cdot,\cdot)$ is coercive according to (4.34) and that
conditions (4.50)–(4.52) hold. The condition number
$\kappa\left(\mathbf{A}\right)$ of the stiffness matrix satisfies the
following bound

$\kappa\left(\mathbf{A}\right)\lesssim h^{-2}$ (4.53)

  * Proof.To estimate the condition number according to (4.49), we use the following Rayleigh quotient characterization of the eigenvalues to bound the maximum (resp. minimum) eigenvalue from above (resp. below):

$\displaystyle\lambda_{\max}=\max_{\mathbf{y}\in\mathbb{R}^{N}\setminus\\{0\\}}\frac{\mathbf{y}^{T}\mathbf{A}\mathbf{y}}{\|\mathbf{y}\|^{2}_{2}},\qquad\lambda_{\min}=\min_{\mathbf{y}\in\mathbb{R}^{N}\setminus\\{0\\}}\frac{\mathbf{y}^{T}\mathbf{A}\mathbf{y}}{\|\mathbf{y}\|^{2}_{2}}$
(4.54)

We bound the Rayleigh quotient from above as follows

$\displaystyle\frac{\mathbf{y}^{T}\mathbf{A}\mathbf{y}}{\|\mathbf{y}\|^{2}_{2}}\sim\frac{a_{h,\bigstar}(v_{h},v_{h})}{h^{-d}\|v_{h}\|^{2}_{\Omega_{h}}}=h^{d}\frac{a_{h,\bigstar}(v_{h},v_{h})}{|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h,\bigstar}}\frac{|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h,\bigstar}}{\|v_{h}\|^{2}_{\Omega_{h}}}\lesssim
h^{d-2}$ (4.55)

where we used, consecutively, the definition (4.48) of the stiffness matrix,
the equivalence (4.52), and continuity of $a_{h,\bigstar}(\cdot,\cdot)$
according to (4.33) and the inverse inequality (4.50). Conversely, to bound
the quotient from below, we infer from the equivalence relation (4.52),
coercivity condition (4.34), and the Poincaré inequality (4.51) that

$\displaystyle\frac{\mathbf{y}^{T}\mathbf{A}\mathbf{y}}{\|\mathbf{y}\|^{2}_{2}}\gtrsim\frac{|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|^{2}_{h,\bigstar}}{h^{-d}\|v_{h}\|^{2}_{\Omega_{h}}}\gtrsim
h^{d}$ (4.56)

Combining the upper and lower bounds, we obtain the desired estimate

$\kappa=\frac{\lambda_{\max}}{\lambda_{\min}}\lesssim\frac{h^{d-2}}{h^{d}}\lesssim
h^{-2}$ (4.57)

∎

######  Remark 4.2.

Comparison to Nitsche’s method without stabilization. We here briefly discuss
which aspects of Section 4.5 do, and do not, apply when Nitsche’s method is
applied in an immersed setting without an additional stability-enhancing
technique. First and foremost, the intrinsic property (4.25) will not be
satisfied. Consequently, the normal gradients on the boundary are not
sufficiently controlled by the gradients in the bulk, such that coercivity
only holds with locally very large (_i.e._ , _a priori_ unbounded) values for
the Nitsche parameter $\beta$. Nevertheless, with the Nitsche parameter large
enough, in the $\beta$-norm both (4.33) and (4.34) hold (_i.e._ , respectively
continuity and coercivity), and only the continuous dependence on the data in
the Lax–Milgram lemma can formally not be verified as the right-hand-side
operator $l_{h}(\cdot)$ is not bounded in the $\beta$-norm such that
$|\mspace{-1.0mu}|\mspace{-1.0mu}|u_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{\beta}$
cannot be _a priori_ bounded. The error estimate provided in (4.39) does not
hold for unstabilized formulations in combination with the $\beta$-norm, as
$\beta\sim h_{T_{\Omega}}^{-1}\gtrsim h^{-1}$ such that (4.43), and
consequently (4.41), do not apply. Note that (4.41) is also employed in the
proof of the bound on the second part of the right-hand-side of (4.40). As
discussed in Section 4.1, this conclusion is also obtained in [1], which shows
that a best approximation property in the $\beta$-norm can be demonstrated
following Céa’s lemma, but that the problem lies in the error estimates
because a bound on the error of the best approximation in the $\beta$-norm
cannot be provided. As already demonstrated in Section 3.1, condition number
bounds cannot be provided and tailored preconditioning is required. Related to
the analysis presented in this section, both (4.50) and (4.51) are contingent
on ghost-penalty stabilization or aggregation. For (4.50) the bound cannot be
provided as there is no _a priori_ upper bound on $\beta$. Nevertheless, this
lack of an upper bound to the eigenvalues is generally not the cause of the
conditioning problems of these methods, and this type of large eigenvalues
only occurs on highly irregularly shaped cut elements. The main problem with
the conditioning is that without additional stabilization, (4.51) does not
hold, and as demonstrated in Section 3.1 it is this lack of a lower bound to
the eigenvalues that generally causes the ill-conditioning of unstabilized
immersed methods.

##  5 Discussion and conclusion

In this review, we considered the problems resulting from small cut elements
in immersed finite element methods, as well as the solutions that have been
developed to counter these problems. This section assesses the current state
of the field and presents a discussion about future research directions.
Similarities and differences between the different techniques that are
discussed in this manuscript are treated in the respective sections, and we do
not provide a qualitative comparison here. This is because, first, the
presented techniques pursue different goals related to the different problems
on small cut elements, _i.e._ , stability and conditioning. Second, the
approaches are not mutually exclusive and combinations can be made. In Section
4 the combination of aggregation and ghost-penalty stabilization is already
discussed, and similarly it is possible to combine it with, _e.g._ , the
minimal stabilization procedure, which fixes the stability but not the
conditioning, with a preconditioning approach.

The field of immersed finite element methods has substantially grown and
developed over the past decade. The employed formulations have matured, and
are currently based on rigorous theoretical and mathematical foundations.
Furthermore, efficient parallel implementations have been developed that
enable the application to large problems into the billions of degrees of
freedom. Certain aspects are still challenging, however, which particularly
manifest in the application to real-world problems. These challenging aspects,
and the corresponding unresolved questions, form the basis of future research
directions that we believe are important for the further advancement of the
field of immersed finite element methods. In this section we discuss these
open research questions in a general sense, and indicate it if these are
specifically related to one of the techniques that are discussed in this
manuscript.

One aspect that requires further consideration to advance the field of
immersed finite element methods, is the detection of small features in the
geometry or the solution. The decoupling of the geometry description from the
finite element mesh makes the consideration of adaptive discretizations
natural to immersed formulations. In this regard, it is notable that grid
refinements in the vicinity of domain boundaries are commonly considered in
the $hp$-adaptive finite cell method [224, 225, 103, 49]. Such refinements are
not required along all boundaries, but mainly in regions where small solution
features such as large gradients are expected. For efficiency purposes, we
therefore consider the further advancement of error estimators and error-
estimation-based adaptivity schemes tailored to immersed formulations as
presented in [226] a desired development. One complication in regard of _a
posteriori_ error estimates is that the bilinear form is only coercive on
discrete spaces, and that the dependence of the bilinear form on the
discretization (both on the mesh size and on the order) complicates saturation
assumptions. Besides the accuracy considerations related to enriched
approximation spaces along small geometrical features, stabilized methods also
require a certain level of geometry resolution, which makes local adaptive
schemes well suited for the application of these techniques to problems on
(scan-based) porous domains, as the refinements reduce the ratio of elements
that are cut. Since element aggregation strongly ties degrees of freedom on
badly cut elements to well posed degrees of freedom, local refinement
procedures can be beneficial in cases where a large portion (or in certain
regions even the majority) of the active mesh consists of badly cut elements.
This is similar for ghost-penalty stabilization terms on domains with
(relative to the mesh) thin structures, which can potentially make the domain
artificially stiff in such cases. Note that in these cases the choice of the
parameters $\eta^{*}$ (with aggregation; see Section 4.3) and $\tau$ (with
ghost-penalty stabilization; see Section 4.4) is particularly critical, and
that in many cases problems can be avoided by a adequate parameter choice.
Nevertheless, the accurate approximation of problems with very thin features
in an unfitted setting remains a difficult problem for any method, and
deserves further attention. Schwarz preconditioners are not affected by the
matter of mesh resolution, but it should be noted that this only pertains to
the solving of the linear system. With a natural boundary condition on the
unfitted boundary and a preconditioner to resolve the conditioning problems, a
relatively low level of mesh resolvement can therefore lead to reasonable
solutions. It should be mentioned that this can still result in artificial
couplings as described in _e.g._ , [47, 163], however. With essential boundary
conditions on the unfitted boundary, such as with flow problems in porous
media, the application of a preconditioner without an additional stabilization
technique will lead to inaccurate approximations, as described in Section 4.

While this manuscript focuses on single-field problems, mixed formulations, in
particular flow problems, are an important application area of immersed finite
element methods. The techniques presented in this manuscript are all
applicable to such problems, and can retain or precondition inf-sup stable
formulations (see, _e.g._ , [198] for aggregation, [21, 52] for ghost penalty,
and [61] for preconditioning). An unresolved topic is the application of
(pointwise) divergence-conforming or geometry-preserving discretizations for
incompressible problems. Compatible pairs of approximation spaces that result
in pointwise divergence-free approximations for boundary-fitted
discretizations, generally lose this property in an unfitted setting.
Furthermore, stabilization techniques complicate the formulation of pointwise
divergence-free methods. Ghost-penalty terms interfere with the enforcement of
incompressibility in the bilinear form, and it has not been thoroughly
investigated how the existing discrete aggregation and extension operators
affect the compatibility of structure-preserving pairs of approximation
spaces. Besides these comments about $H$(div)-conforming discretizations for
the conservation of geometry, similar considerations can be made about
$H$(curl)-conforming discretizations when solving the Maxwell equations in
electromagnetism.

An aspect that relates to both of the previous topics of mesh adaptivity and
flow problems is the consideration of anisotropic grid refinements. To
accurately capture boundary layers at a reasonable computational cost,
boundary layer meshes are commonly anisotropic with a much finer resolution
normal to the boundary. While the aspect of boundary layers in immersed
methods has already been considered in multiple studies (see, _e.g._ , [227,
24]), present understanding of this aspect is incomplete, and this is a topic
that requires further investigation.

While only briefly considered in Section 2.2.3, a topic that is strongly
related to small cut elements is explicit time integration for problems in
dynamics. This is particularly interesting for crash simulations, as such
problems are commonly solved by explicit time steppers [150]. Although tools
exist that result in stable time steps that are not affected by the cut
elements (_i.e._ , element aggregation, ghost penalty terms in the mass
matrix, or lumping of the mass matrix [149, 150, 145]), the research on this
topic is not exhaustive. In particular, a comprehensive analysis of stable
time-step sizes with different choices for the mass and stiffness forms
(_i.e._ , boundary condition enforcement technique, ghost-penalty
stabilization of both forms, lumping of the mass matrix, etc.) does, to the
authors’ knowledge, not yet exist.

Another open question deals with the choice, or the optimization, of
parameters in the presented approaches. This does not only concern, _e.g._ ,
the parameters $\eta^{*}$ and $\tau$ that are explicitly contained in element
aggregation and the ghost penalty, but also hidden parameters (or choices) in
the setup such as the choice of the index blocks with Schwarz preconditioning.
While these parameters generally do not affect the asymptotic behavior, well-
chosen (or, conversely, poorly-chosen) parameters can significantly affect the
performance pre-asymptotically. We therefore consider an optimal choice of
(hidden) parameters an important topic for future research. This specifically
applies to the dependence of this optimality on the basis-function type and
the discretization order and, in particular for discretizations with local
refinements or local order elevations, to the option of choosing parameters
locally or element-wise (similar to the Nitsche parameter without
stabilization).

Over the past years, many well-developed, well-documented, and user-friendly
open source codes have become available that aid the application of immersed
finite element methods, a comprehensive list of which is beyond the scope of
this work. These toolboxes support all facets of immersed finite element
computations, such as pre-processing steps (in particular integration and the
application of boundary conditions), the solution process itself, and also
post-processing, where often a subtriangulation is employed like this is also
done for higher-order boundary-fitted finite elements. Something that is not
yet available is a graphical user interface that makes immersed finite element
methods accessible also to non-expert users. Besides the developments in open
source codes, also commercial software has adopted (aspects of) immersed
methods. To provide some examples: in Ansys Fluent the Immersed Boundary
Method (IBM) is available [63, 228]; Ansys LS-DYNA facilitates the application
of trimmed isogeometric analysis in both shell and solid models for static and
dynamic analyses [229, 230]; in Abaqus it is possible to (partially) include
and exclude elements from a computation with _element progressive activation_
[231, 232], which is mainly intended for the simulation of additive
manufacturing; and the platform Hyperganic [233] employs immersed finite
element techniques to perform simulations. The implementation of the
stabilization methods discussed in this work can require tools that are not
generally available in existing software and codes. For instance, face-based
ghost-penalty stabilization requires the computation of jumps in normal
gradients up to the order of approximation. For quadratic and higher-order
discretizations, this functionality is generally not available, even for
software frameworks that include the possibility to perform face loops and
compute jumps for standard Discontinuous Galerkin (DG) formulations. The
implementation of patch-based ghost-penalty stabilization or aggregation
require the implementation of an aggregation routine and the computation of
aggregate-wise terms. These terms can be implemented using standard element-
wise terms if the finite element software provides functionality to enforce
linear constraints. Many codes already provide this functionality as this is
essential for local grid refinements and periodic boundary conditions (_e.g._
, deal.ii [234], fenics [235], gridap [236], fempar [200], netgen [237], or
dune [238]). While (geometric) multigrid routines are intrusive to the finite
element software, the implementation of tailored Schwarz preconditioners is
exceptionally non-intrusive and only requires the matrix itself and the basis
function connectivity. Furthermore, as additive Schwarz preconditioning is a
well-established concept also in boundary-fitted finite element methods, these
preconditioners are available in many codes. The construction of the index
blocks is not generally the same in existing software, however, and it should
be noted that the effectivity for immersed methods relies on the specific
construction of these blocks. Additionally, advanced procedures in which
_e.g._ , index blocks are only devised for basis functions with support on cut
elements, require the availability of this information or a routine to extract
it.

By virtue of the stabilization, aggregation, and preconditioning techniques as
discussed in this review, immersed finite element methods have developed into
reliable and versatile computational analysis instruments. In our assessment,
the potential of immersed methods and, in particular, their ability to provide
new simulation workflows with minimal user intervention, has not been fully
exploited yet. Another recommendation therefore pertains to the exploration of
immersed methods in new application fields, such as hyperbolic systems,
diffuse interface models, poro-elasticity, computational aeroacoustics, and
phase field fracture.

#### Conflict of interest statement

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## References

  * [1] F. de Prenter, C. Lehrenfeld, and A. Massing.  A note on the stability parameter in Nitsche’s method for unfitted boundary value problems.  Computers & Mathematics with Applications, 75(12):4322–4336, 2018.
  * [2] A. Düster, J. Parvizian, Z. Yang, and E. Rank.  The finite cell method for three-dimensional problems of solid mechanics.  Computer Methods in Applied Mechanics and Engineering, 197(45-48):3768–3782, 2008.
  * [3] D. Schillinger, L. Dedè, M.A. Scott, J.A. Evans, M.J. Borden, E. Rank, and T.J.R. Hughes.  An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces.  Computer Methods in Applied Mechanics and Engineering, 249-252:116–150, 2012.
  * [4] D. Schillinger, M. Ruess, N. Zander, Y. Bazilevs, A. Düster, and E. Rank.  Small and large deformation analysis with the _p_ \- and B-spline versions of the Finite Cell Method.  Computational Mechanics, 50:445–478, 2012.
  * [5] D. Schillinger, A. Düster, and E. Rank.  The _hp_ -_d_ -adaptive finite cell method for geometrically nonlinear problems of solid mechanics.  International Journal for Numerical Methods in Engineering, 89(9):1171–1202, 2012.
  * [6] N. Zander, S. Kollmannsberger, M. Ruess, Z. Yosibash, and E. Rank.  The Finite Cell Method for linear thermoelasticity.  Computers & Mathematics with Applications, 64(11):3527–3541, 2012.
  * [7] M. Joulaian, S. Duczek, U. Gabbert, and A. Düster.  Finite and spectral cell method for wave propagation in heterogeneous materials.  Computational Mechanics, 54:661–675, 2014.
  * [8] M. Ruess, D. Schillinger, A.I. Özcan, and E. Rank.  Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries.  Computer Methods in Applied Mechanics and Engineering, 269:46–71, 2014.
  * [9] E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, and A. Düster.  Geometric modeling, isogeometric analysis and the finite cell method.  Computer Methods in Applied Mechanics and Engineering, 249-252:104–115, 2012.
  * [10] M. Ruess, D. Schillinger, Y. Bazilevs, V. Varduhn, and E. Rank.  Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method.  International Journal for Numerical Methods in Engineering, 95(10):811–846, 2013.
  * [11] E. Rank, S. Kollmannsberger, C. Sorger, and A. Düster.  Shell Finite Cell Method: A high order fictitious domain approach for thin-walled structures.  Computer Methods in Applied Mechanics and Engineering, 200(45-46):3200–3209, 2011.
  * [12] R. Schmidt, R. Wüchner, and K.-U. Bletzinger.  Isogeometric analysis of trimmed NURBS geometries.  Computer Methods in Applied Mechanics and Engineering, 241-244:93–111, 2012.
  * [13] Y. Guo and M. Ruess.  Weak Dirichlet boundary conditions for trimmed thin isogeometric shells.  Computers & Mathematics with Applications, 70(7):1425–1440, 2015\.
  * [14] A.M. Bauer, M. Breitenberger, B. Philipp, R. Wüchner, and K.-U. Bletzinger.  Embedded structural entities in NURBS-based isogeometric analysis.  Computer Methods in Applied Mechanics and Engineering, 325:198–218, 2017.
  * [15] Y. Guo, M. Ruess, and D. Schillinger.  A parameter-free variational coupling approach for trimmed isogeometric thin shells.  Computational Mechanics, 59:693–715, 2017.
  * [16] J.E. Dolbow and I. Harari.  An efficient finite element method for embedded interface problems.  International Journal for Numerical Methods in Engineering, 78(2):229–252, 2009.
  * [17] D. Schillinger and E. Rank.  An unfitted _hp_ -adaptive finite element method based on hierarchical B-splines for interface problems of complex geometry.  Computer Methods in Applied Mechanics and Engineering, 200(47-48):3358–3380, 2011.
  * [18] C. Annavarapu, M. Hautefeuille, and J.E. Dolbow.  A robust Nitsche’s formulation for interface problems.  Computer Methods in Applied Mechanics and Engineering, 225-228:44–54, 2012.
  * [19] M. Hautefeuille, C. Annavarapu, and J.E. Dolbow.  Robust imposition of Dirichlet boundary conditions on embedded surfaces.  International Journal for Numerical Methods in Engineering, 90(1):40–64, 2012.
  * [20] E. Wadbro, S. Zahedi, G. Kreiss, and M. Berggren.  A uniformly well-conditioned, unfitted Nitsche method for interface problems.  BIT Numerical Mathematics, 53:791–820, 2013.
  * [21] A. Massing, M.G. Larson, A. Logg, and M.E. Rognes.  A stabilized Nitsche fictitious domain method for the Stokes problem.  Journal of Scientific Computing, 61:604–628, 2014.
  * [22] B. Schott and W.A. Wall.  A new face-oriented stabilized XFEM approach for 2D and 3D incompressible Navier-Stokes equations.  Computer Methods in Applied Mechanics and Engineering, 276:233–265, 2014.
  * [23] P. Hansbo, M.G. Larson, and S. Zahedi.  A cut finite element method for a Stokes interface problem.  Applied Numerical Mathematics, 85:90–114, 2014.
  * [24] F. Xu, D. Schillinger, D. Kamensky, V. Varduhn, C. Wang, and M.-C. Hsu.  The tetrahedral finite cell method for fluids: Immersogeometric analysis of turbulent flow around complex geometries.  Computers and Fluids, 141:135–154, 2016.
  * [25] B. Schott, U. Rasthofer, V. Gravemeier, and W.A. Wall.  A face-oriented stabilized Nitsche-type extended variational multiscale method for incompressible two-phase flow.  International Journal for Numerical Methods in Engineering, 104(7):721–748, 2015.
  * [26] B. Schott, S. Shahmiri, R. Kruse, and W.A. Wall.  A stabilized Nitsche-type extended embedding mesh approach for 3D low- and high-Reynolds-number flows.  International Journal for Numerical Methods in Fluids, 82(6):289–315, 2016.
  * [27] M.-C. Hsu, C. Wang, F. Xu, A.J. Herrema, and A. Krishnamurthy.  Direct immersogeometric fluid flow analysis using B-rep CAD models.  Computer Aided Geometric Design, 43:143–158, 2016.
  * [28] T. Hoang, C.V. Verhoosel, F. Auricchio, E.H. van Brummelen, and A. Reali.  Mixed Isogeometric Finite Cell Methods for the Stokes problem.  Computer Methods in Applied Mechanics and Engineering, 316:400–423, 2017.
  * [29] A. Massing, B. Schott, and W.A. Wall.  A stabilized Nitsche cut finite element method for the Oseen problem.  Computer Methods in Applied Mechanics and Engineering, 328:262–300, 2018.
  * [30] M. Winter, B. Schott, A. Massing, and W.A. Wall.  A Nitsche cut finite element method for the Oseen problem with general Navier boundary conditions.  Computer Methods in Applied Mechanics and Engineering, 330:220–252, 2018.
  * [31] Y. Bazilevs, M.-C. Hsu, and M.A. Scott.  Isogeometric fluid–structure interaction analysis with emphasis on non-matching discretizations, and with application to wind turbines.  Computer Methods in Applied Mechanics and Engineering, 249:28–41, 2012.
  * [32] E. Burman and M.A. Fernández.  An unfitted Nitsche method for incompressible fluid–structure interaction using overlapping meshes.  Computer Methods in Applied Mechanics and Engineering, 279:497–514, 2014.
  * [33] T. Rüberg and F. Cirak.  A fixed-grid b-spline finite element technique for fluid–structure interaction.  International Journal for Numerical Methods in Fluids, 74(9):623–660, 2014.
  * [34] A. Massing, M.G. Larson, A. Logg, and M.E. Rognes.  A Nitsche-based cut finite element method for a fluid–structure interaction problem.  Communications in Applied Mathematics and Computational Science, 10(2):97–120, 2015.
  * [35] C. Kadapa, W.G. Dettmer, and D. Perić.  A fictitious domain/distributed Lagrange multiplier based fluid–structure interaction scheme with hierarchical B-spline grids.  Computer Methods in Applied Mechanics and Engineering, 301:1–27, 2016.
  * [36] C. Kadapa, W.G. Dettmer, and D. Perić.  A stabilised immersed boundary method on hierarchical b-spline grids for fluid–rigid body interaction with solid–solid contact.  Computer Methods in Applied Mechanics and Engineering, 318:242–269, 2017.
  * [37] Y. Wang, P.K. Jimack, and M.A. Walkley.  A one-field monolithic fictitious domain method for fluid–structure interactions.  Computer Methods in Applied Mechanics and Engineering, 317:1146–1168, 2017.
  * [38] M.C.H. Wu, D. Kamensky, C. Wang, A.J. Herrema, F. Xu, M.S. Pigazzini, A. Verma, A.L. Marsden, Y. Bazilevs, and M.-C. Hsu.  Optimizing fluid–structure interaction systems with immersogeometric analysis and surrogate modeling: Application to a hydraulic arresting gear.  Computer Methods in Applied Mechanics and Engineering, 316:668–693, 2017.
  * [39] C. Kadapa, W.G. Dettmer, and D. Perić.  A stabilised immersed framework on hierarchical b-spline grids for fluid–flexible structure interaction with solid–solid contact.  Computer Methods in Applied Mechanics and Engineering, 335:472–489, 2018.
  * [40] M.-C. Hsu, D. Kamensky, Y. Bazilevs, M.S. Sacks, and T.J.R. Hughes.  Fluid–structure interaction analysis of bioprosthetic heart valves: significance of arterial wall deformation.  Computational Mechanics, 54:1055–1071, 2014.
  * [41] D. Kamensky, M.-C. Hsu, D. Schillinger, J.A. Evans, A. Aggarwal, Y. Bazilevs, M.S. Sacks, and T.J.R. Hughes.  An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves.  Computer Methods in Applied Mechanics and Engineering, 284:1005–1053, 2015.
  * [42] M.-C. Hsu, D. Kamensky, F. Xu, J. Kiendl, C. Wang, M.C.H. Wu, J. Mineroff, A. Reali, Y. Bazilevs, and M.S. Sacks.  Dynamic and fluid–structure interaction simulations of bioprosthetic heart valves using parametric design with T-splines and Fung-type material models.  Computational Mechanics, 55:1211–1225, 2015.
  * [43] D. Kamensky, M.-C. Hsu, Y. Yu, J.A. Evans, M.S. Sacks, and T.J.R. Hughes.  Immersogeometric cardiovascular fluid–structure interaction analysis with divergence-conforming B-splines.  Computer Methods in Applied Mechanics and Engineering, 314:408–472, 2017.
  * [44] Z. Yang, S. Kollmannsberger, A. Düster, M. Ruess, E. Grande Garcia, R. Burgkart, and E. Rank.  Non-standard bone simulation: interactive numerical analysis by computational steering.  Computing and Visualization in Science, 14:207–216, 2011.
  * [45] M. Ruess, D. Tal, N. Trabelsi, Z. Yosibash, and E. Rank.  The finite cell method for bone simulations: verification and validation.  Biomechanics and Modeling in Mechanobiology, 11:425–437, 2012.
  * [46] A. Düster, H.-G. Sehlhorst, and E. Rank.  Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method.  Computational Mechanics, 50:413–431, 2012.
  * [47] C.V. Verhoosel, G.J. Van Zwieten, B. Van Rietbergen, and R. De Borst.  Image-based goal-oriented adaptive isogeometric analysis with application to the micro-mechanical modeling of trabecular bone.  Computer Methods in Applied Mechanics and Engineering, 284:138–164, 2015.
  * [48] V. Varduhn, M.-C. Hsu, M. Ruess, and D. Schillinger.  The tetrahedral finite cell method: Higher-order immersogeometric analysis on adaptive non-boundary-fitted meshes.  International Journal for Numerical Methods in Engineering, 107(12):1054–1079, 2016.
  * [49] M. Elhaddad, N. Zander, T. Bog, L. Kudela, S. Kollmannsberger, J. Kirschke, T. Baum, M. Ruess, and E. Rank.  Multi-level _hp_ -finite cell method for embedded interface problems with application in biomechanics.  International Journal for Numerical Methods in Biomedical Engineering, 34(4):e2951, 2017.
  * [50] S. Duczek, H. Berger, and U. Gabbert.  The Finite Pore Method: a new approach to evaluate gas pores in cast parts by combining computed tomography and the finite cell method.  International Journal of Cast Metals Research, 28(4):221–228, 2015\.
  * [51] M. Würkner, S. Duczek, H. Berger, H. Köppe, and U. Gabbert.  A software platform for the analysis of porous die-cast parts using the finite cell method, pages 327–341.  Springer, 2018.
  * [52] T. Hoang, C.V. C.V. Verhoosel, C.-Z. Qin, F. Auricchio, A. Reali, and E.H. van Brummelen.  Skeleton-stabilized immersogeometric analysis for incompressible viscous flow problems.  Computer Methods in Applied Mechanics and Engineering, 344:421–450, 2019.
  * [53] J. Parvizian, A. Düster, and E. Rank.  Topology optimization using the finite cell method.  Optimization and Engineering, 13:57–78, 2012.
  * [54] E. Nadal, J.J. Ródenas, J. Albelda, M. Tur, J.E. Tarancón, and F.J. Fuenmayor.  Efficient Finite Element Methodology Based on Cartesian Grids: Application to Structural Shape Optimization.  Abstract and Applied Analysis, page 953786, 2013.
  * [55] N.P. van Dijk, K. Maute, M. Langelaar, and F. van Keulen.  Level-set methods for structural topology optimization: a review.  Structural and Multidisciplinary Optimization, 48:437–472, 2013\.
  * [56] K. Bandara, T. Rüberg, and F. Cirak.  Shape optimisation with multiresolution subdivision surfaces and immersed finite elements.  Computer Methods in Applied Mechanics and Engineering, 300:510–539, 2016.
  * [57] J.P. Groen, M. Langelaar, O. Sigmund, and M. Ruess.  Higher-order multi-resolution topology optimization using the finite cell method.  International Journal for Numerical Methods in Engineering, 110(10):903–920, 2017.
  * [58] C.H. Villanueva and K. Maute.  CutFEM topology optimization of 3D laminar incompressible flow problems.  Computer Methods in Applied Mechanics and Engineering, 320:444–473, 2017.
  * [59] E. Burman, P. Hansbo, and M.G. Larson.  A cut finite element method with boundary value correction.  Mathematics of Computation, 87:633–657, 2018.
  * [60] J.N. Jomo, F. de Prenter, M. Elhaddad, D. D’Angella, C.V. Verhoosel, S. Kollmannsberger, J.S. Kirschke, V. Nübel, E.H. van Brummelen, and E. Rank.  Robust and parallel scalable iterative solutions for large-scale finite cell analyses.  Finite Elements in Analysis and Design, 163:14–30, 2019. Figures reproduced with permission from Elsevier.
  * [61] F. de Prenter, C.V. Verhoosel, and E.H. van Brummelen.  Preconditioning immersed isogeometric finite element methods with application to flow problems.  Computer Methods in Applied Mechanics and Engineering, 348:604–631, 2019. Figures reproduced with permission from Elsevier.
  * [62] C.S. Peskin.  Flow patterns around heart valves: A numerical method.  Journal of Computational Physics, 10(2):252–271, 1972.
  * [63] R. Mittal and G. Iaccarino.  Immersed Boundary Methods.  Annual Review of Fluid Mechanics, 37:239–261, 2005.
  * [64] J.M. Melenk and I Babuška.  The partition of unity finite element method: Basic theory and applications.  Computer Methods in Applied Mechanics and Engineering, 139(1-4):289–314, 1996.
  * [65] I. Babuška and J.M. Melenk.  The partition of unity method.  International Journal for Numerical Methods in Engineering, 40(4):727–758, 1997.
  * [66] C.A. Duarte, I. Babuška, and J.T. Oden.  Generalized finite element methods for three-dimensional structural mechanics problems.  Computers and Structures, 77(2):215–232, 2000.
  * [67] T. Strouboulis, I. Babuška, and K. Copps.  The design and analysis of the Generalized Finite Element Method.  Computer Methods in Applied Mechanics and Engineering, 181(1-3):43–69, 2000.
  * [68] T. Belytschko and T. Black.  Elastic crack growth in finite elements with minimal remeshing.  International Journal for Numerical Methods in Engineering, 45(5):601–620, 1999.
  * [69] N. Moës, J.E. Dolbow, and T. Belytschko.  A finite element method for crack growth without remeshing.  International Journal for Numerical Methods in Engineering, 46(1):131–150, 1999.
  * [70] A. Hansbo and P. Hansbo.  An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems.  Computer Methods in Applied Mechanics and Engineering, 191(47-48):5537–5552, 2002.
  * [71] J. Parvizian, A. Düster, and E. Rank.  Finite cell method.  Computational Mechanics, 41:121–133, 2007.
  * [72] D. Schillinger and M. Ruess.  The Finite Cell Method: A Review in the Context of Higher-Order Structural Analysis of CAD and Image-Based Geometric Models.  Archives of Computational Methods in Engineering, 22:391–455, 2015\.
  * [73] E. Burman and P. Hansbo.  Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method.  Applied Numerical Mathematics, 62(4):328–341, 2012.
  * [74] E. Burman, S. Claus, P. Hansbo, M.G. Larson, and A. Massing.  CutFEM: Discretizing geometry and partial differential equations.  International Journal for Numerical Methods in Engineering, 104(7):472–501, 2015.
  * [75] E. Burman, P. Hansbo, M.G. Larson, and A. Massing.  Cut finite element methods for partial differential equations on embedded manifolds of arbitrary codimensions.  ESAIM: Mathematical Modelling and Numerical Analysis (ESAIM: M2AN), 52(6):2247–2282, 2018.
  * [76] E. Burman.  Ghost penalty.  Comptes Rendus Mathematique, 348(21-22):1217–1220, 2010.
  * [77] S. Badia, F. Verdugo, and A.F. Martín.  The aggregated unfitted finite element method for elliptic problems.  Computer Methods in Applied Mechanics and Engineering, 336:533–553, 2018.
  * [78] S. Badia, E. Neiva, and F. Verdugo.  Linking ghost penalty and aggregated unfitted methods.  Computer Methods in Applied Mechanics and Engineering, 388:114232, 2022.
  * [79] O. Marco, R. Sevilla, Y. Zhang, J.J. Ródenas, and M. Tur.  Exact 3D boundary representation in finite element analysis based on cartesian grids independent of the geometry.  International Journal for Numerical Methods in Engineering, 103(6):445–468, 2015.
  * [80] K. Höllig, U. Reif, and J. Wipper.  Weighted extended B-spline approximation of Dirichlet problems.  SIAM Journal on Numerical Analysis, 39(2):442–462, 2001.
  * [81] K. Höllig, C. Apprich, and A. Streit.  Introduction to the Web-method and its applications.  Advances in Computational Mathematics, 23:215–237, 2005.
  * [82] R.A.K. Sanches, P.B. Bornemann, and F. Cirak.  Immersed b-spline (i-spline) finite element method for geometrically complex domains.  Computer Methods in Applied Mechanics and Engineering, 200(13-16):1432–1445, 2011.
  * [83] A. Johansson and M.G. Larson.  A high order discontinuous Galerkin Nitsche method for elliptic problems with fictitious boundary.  Numerische Mathematik, 123:607–628, 2013.
  * [84] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  Computer Methods in Applied Mechanics and Engineering, 194(39-41):4135–4195, 2005.
  * [85] B. Marussig and T.J.R. Hughes.  A Review of Trimming in Isogeometric Analysis: Challenges, Data Exchange and Simulation Aspects.  Archives of Computational Methods in Engineering, 25:1059–1127, 2018.
  * [86] A. Main and G. Scovazzi.  The shifted boundary method for embedded domain computations. Part II: Linear advection–diffusion and incompressible Navier–Stokes equations.  Journal of Computational Physics, 372:996–1026, 2018.
  * [87] A. Main and G. Scovazzi.  The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems.  Journal of Computational Physics, 372:972–995, 2018.
  * [88] N.M. Atallah, C. Canuto, and G. Scovazzi.  The second-generation Shifted Boundary Method and its numerical analysis.  Computer Methods in Applied Mechanics and Engineering, 372:113341, 2020.
  * [89] N.M. Atallah, C. Canuto, and G. Scovazzi.  Analysis of the shifted boundary method for the Poisson problem in domains with corners.  Mathematics of Computation, 90:2041–2069, 2021.
  * [90] J.H. Bramble, T. Dupont, and V. Thomée.  Projection methods for Dirichlet’s problem in approximating polygonal domains with boundary-value corrections.  Mathematics of Computation, 26:869–879, 1972.
  * [91] E. Burman and A. Ern.  An unfitted hybrid high-order method for elliptic interface problems.  SIAM Journal on Numerical Analysis, 56(3):1525–1546, 2018.
  * [92] S. Badia, J. Droniou, and L. Yemm.  Conditioning of a Hybrid High-Order scheme on Meshes with Small Faces.  Journal of Scientific Computing, 92(71), 2022.
  * [93] A. Ern and J.-L. Guermond.  Theory and Practice of Finite Elements.  Springer, 2004.
  * [94] L.C. Evans.  Partial differential equations.  American Mathematical Society, 2010.
  * [95] J. Nitsche.  Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind.  In Abhandlungen aus dem mathematischen Seminar der Universität Hamburg, pages 9–15, 1971.
  * [96] A. Embar, J.E. Dolbow, and I. Harari.  Imposing Dirichlet boundary conditions with Nitsche’s method and spline based finite elements.  International Journal for Numerical Methods in Engineering, 83(7):877–898, 2010.
  * [97] T. Warburton and J.S. Hesthaven.  On the constants in _hp_ -finite element trace inverse inequalities.  Computer Methods in Applied Mechanics and Engineering, 192(25):2765–2773, 2003.
  * [98] J.A. Evans and T.J.R. Hughes.  Explicit trace inequalities for isogeometric analysis and parametric hexahedral finite elements.  Numerische Mathematik, 123:259–290, 2013.
  * [99] J.A. Cottrell, T.J.R. Hughes, and Y. Bazilevs.  Isogeometric Analysis: Toward Integration of CAD and FEA.  Wiley, 2009.
  * [100] N. Korshunova, J. Jomo, G. Léko, D. Reznik, P. Balàzs, and S. Kollmannsberger.  Image-based material characterization of complex microarchitectured additively manufactured structures.  Computers & Mathematics with Applications, 80(11):2462–2480, 2020.
  * [101] A. Abedian, J. Parvizian, A. Düster, H. Khademyzadeh, and E. Rank.  Performance of different integration schemes in facing discontinuities in the Finite Cell Method.  International Journal of Computational Methods, 10(3):1350002, 2013\.
  * [102] T.-P. Fries and T. Belytschko.  The extended/generalized finite element method: An overview of the method and its applications.  International Journal for Numerical Methods in Engineering, 84(3):253–304, 2010.
  * [103] A. Düster, E. Rank, and B. Szabó.  The p-Version of the Finite Element and Finite Cell Methods, pages 1–35.  Wiley, 2017.
  * [104] A. Stavrev, L.H. Nguyen, R. Shen, V. Varduhn, M. Behr, S. Elgeti, and D. Schillinger.  Geometrically accurate, efficient, and flexible quadrature techniques for the tetrahedral finite cell method.  Computer Methods in Applied Mechanics and Engineering, 310:646–673, 2016.
  * [105] S.C. Divi, C.V. Verhoosel, F. Auricchio, A. Reali, and E.H. van Brummelen.  Error-estimate-based adaptive integration for immersed isogeometric analysis.  Computers & Mathematics with Applications, 80(11):2481–2516, 2020.
  * [106] T.-P. Fries and S. Omerović.  Higher-order accurate integration of implicit geometries.  International Journal for Numerical Methods in Engineering, 106(5):323–371, 2016.
  * [107] T.-P. Fries, S. Omerovi c, D. Schöllhammer, and J. Steidl.  Higher-order meshing of implicit geometries – Part I: Integration and interpolation in cut elements.  Computer Methods in Applied Mechanics and Engineering, 313:759–784, 2017.
  * [108] S. Omerović and T.-P. Fries.  Conformal higher-order remeshing schemes for implicitly defined interface problems.  International Journal for Numerical Methods in Engineering, 109(6):763–789, 2017.
  * [109] T.-P. Fries, S. Omerović, D. Schöllhammer, and J.W. Stanford.  A fictitious domain method with higher-order accurate integration in cut elements.  In Proceedings of the 6th European Conference on Computational Mechanics (ECCM 6) and the 7th European Conference on Computational Fluid Dynamics (ECFD 7), pages 3994–3405, 2018.
  * [110] R. Sevilla, S. Fernández-méndez, and A. Huerta.  NURBS-enhanced finite element method (NEFEM).  International Journal for Numerical Methods in Engineering, 76(1):56–83, 2008.
  * [111] R. Sevilla, S. Fernández-méndez, and A. Huerta.  NURBS-Enhanced Finite Element Method (NEFEM).  Archives of Computational Methods in Engineering, 18(441):441–484, 2011.
  * [112] C. Lehrenfeld.  High order unfitted finite element methods on level set domains using isoparametric mappings.  Computer Methods in Applied Mechanics and Engineering, 300:716–733, 2016.
  * [113] L. Kudela, N. Zander, T. Bog, S. Kollmannsberger, and E. Rank.  Efficient and accurate numerical quadrature for immersed boundary methods.  Advanced Modeling and Simulation in Engineering Sciences, 2(10), 2015.
  * [114] L. Kudela, N. Zander, S. Kollmannsberger, and E. Rank.  Smart octrees: Accurately integrating discontinuous functions in 3D.  Computer Methods in Applied Mechanics and Engineering, 306:406–426, 2016.
  * [115] S. Hubrich, P. Di Stolfo, L. Kudela, S. Kollmannsberger, E. Rank, A. Schröder, and A. Düster.  Numerical integration of discontinuous functions: moment fitting and smart octree.  Computational Mechanics, 60:863–881, 2017.
  * [116] S.E. Mousavi and N. Sukumar.  Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons.  Computational Mechanics, 47:535–554, 2011.
  * [117] M. Joulaian, S. Hubrich, and A. Düster.  Numerical integration of discontinuities on arbitrary domains based on moment fitting.  Computational Mechanics, 57:979–999, 2016.
  * [118] S. Hubrich and A. Düster.  Numerical integration for nonlinear problems of the finite cell method using an adaptive scheme based on moment fitting.  Computers & Mathematics with Applications, 77(7):1983–1997, 2019\.
  * [119] A. Düster and S Hubrich.  Adaptive Integration of Cut Finite Elements and Cells for Nonlinear Structural Analysis, pages 31–73.  Springer, 2020.
  * [120] G. Ventura.  On the elimination of quadrature subcells for discontinuous functions in the eXtended Finite-Element Method.  International Journal for Numerical Methods in Engineering, 66(5):761–795, 2006.
  * [121] A. Abedian and A Düster.  Equivalent Legendre polynomials: Numerical integration of discontinuous functions in the finite element methods.  Computer Methods in Applied Mechanics and Engineering, 343:690–720, 2019.
  * [122] G. Ventura, R. Gracie, and T. Belytschko.  Fast integration and weight function blending in the extended finite element method.  International Journal for Numerical Methods in Engineering, 77(1):1–29, 2009.
  * [123] T. Jonsson, M.G. Larson, and K. Larsson.  Cut finite element methods for elliptic problems on multipatch parametric surfaces.  Computer Methods in Applied Mechanics and Engineering, 324:366–394, 2017.
  * [124] R.I. Saye.  High-order quadrature methods for implicitly defined surfaces and volumes in hyperrectangles.  SIAM Journal on Scientific Computing, 37(2):A993–A1019, 2015\.
  * [125] E.B. Chin and N. Sukumar.  An efficient method to integrate polynomials over polytopes and curved solids.  Computer Aided Geometric Design, 82:101914, 2020.
  * [126] A. Longva, F. Löschner, T. Kugelstadt, J.A. Fernández-Fernández, and J. Bender.  Higher-order finite elements for embedded simulation.  ACM Transactions on Graphics, 39(6):Article 181, 2020.
  * [127] G. Strang and G.J. Fix.  An Analysis of the Finite Element Method.  Prentice-Hall, 1973.
  * [128] G. Strang.  Variational crimes in the finite element method, pages 689–710.  Academic Press, 1972.
  * [129] E. Burman, P. Hansbo, M.G. Larson, and S. Zahedi.  Cut finite element methods for coupled bulk–surface problems.  Numerische Mathematik, 133:203–231, 2016.
  * [130] E. Burman, P. Hansbo, M.G. Larson, and A. Massing.  A cut discontinuous Galerkin method for the Laplace–Beltrami operator.  IMA Journal of Numerical Analysis, 37(1):138–169, 2017.
  * [131] I. Babuška.  The finite element method with penalty.  Mathematics of Computation, 27:221–228, 1973.
  * [132] L. van Miegroet and P. Duysinx.  Stress concentration minimization of 2D filets using X-FEM and level set description.  Structural and Multidisciplinary Optimization, 33:425–438, 2007\.
  * [133] W. Jiang, C. Annavarapu, J.E. Dolbow, and I. Harari.  A robust Nitsche’s formulation for interface problems with spline-based finite elements.  International Journal for Numerical Methods in Engineering, 104(7):676–696, 2015.
  * [134] S. Fernández-Méndez and A. Huerta.  Imposing essential boundary conditions in mesh-free methods.  Computer Methods in Applied Mechanics and Engineering, 193(12-14):1257–1275, 2004.
  * [135] E. Burman.  A penalty-free nonsymmetric Nitsche-type method for the weak imposition of boundary conditions.  SIAM Journal on Numerical Analysis, 50(4):1959–1981, 2012.
  * [136] T. Boiveau and E. Burman.  A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity.  IMA Journal of Numerical Analysis, 36(2):770–795, 2016.
  * [137] D. Schillinger, I. Harari, M.-C. Hsu, D. Kamensky, S.K.F. Stoter, Y. Yu, and Y. Zhao.  The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements.  Computer Methods in Applied Mechanics and Engineering, 309:625–652, 2016.
  * [138] I. Babuška.  The finite element method with Lagrangian multipliers.  Numerische Mathematik, 20:179–192, 1973.
  * [139] D.P. Bertsekas.  Constrained optimization and Lagrange multiplier methods.  Athena Scientific, 1982.
  * [140] E. Burman and P. Hansbo.  Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method.  Computer Methods in Applied Mechanics and Engineering, 199(41-44):2680–2686, 2010.
  * [141] J. Baiges, R. Codina, F. Henke, S. Shahmiri, and W.A. Wall.  A symmetric method for weakly imposing Dirichlet boundary conditions in embedded finite element meshes.  International Journal for Numerical Methods in Engineering, 90(5):636–658, 2012.
  * [142] M. Tur, J. Albelda, E. Nadal, and J.J. Ródenas.  Imposing Dirichlet boundary conditions in hierarchical Cartesian meshes by means of stabilized Lagrange multipliers.  International Journal for Numerical Methods in Engineering, 98(6):399–417, 2014.
  * [143] S. Kollmannsberger, A. Özcan, J. Baiges, M. Ruess, E. Rank, and A. Reali.  Parameter-free, weak imposition of Dirichlet boundary conditions and coupling of trimmed and non-conforming patches.  International Journal for Numerical Methods in Engineering, 101(9):670–699, 2015.
  * [144] F. de Prenter, C.V. Verhoosel, G.J. van Zwieten, and E.H. van Brummelen.  Condition number analysis and preconditioning for the finite cell method.  Computer Methods in Applied Mechanics and Engineering, 316:297–327, 2017.
  * [145] E. Burman, P. Hansbo, and M.G. Larson.  Explicit time stepping for the wave equation using CutFEM with discrete extension.  SIAM Journal on Scientific Computing, 44(3):A1254–A1289, 2022\.
  * [146] I. Harari and U. Albocher.  Spectral investigations of Nitsche’s method.  Journal of Scientific Computing, 81:1472–1492, 2018.
  * [147] I. Harari and U. Albocher.  Complementary Solutions of Nitsche’s Method.  Finite Elements in Analysis and Design, 145:20–31, 2019.
  * [148] U. Albocher and I. Harari.  Spectral aspects of Nitsche’s method on nonconforming meshes.  Mechanics Research Communications, 112:103611, 2021.
  * [149] L.F. Leidinger, M. Breitenberger, A.M. Bauer, S. Hartmann, R. Wüchner, K.-U. Bletzinger, F. Duddeck, and L. Song.  Explicit dynamic isogeometric B-Rep analysis of penalty-coupled trimmed NURBS shells.  Computer Methods in Applied Mechanics and Engineering, 351:891–927, 2019.
  * [150] L.F. Leidinger.  Explicit Isogeometric B-Rep Analysis for Nonlinear Dynamic Crash Simulations.  PhD thesis, Technical University of Munich, 2020.
  * [151] G.H. Golub and C.F. van Loan.  Matrix computations.  John Hopkins University Press, 2013.
  * [152] Y. Saad.  Iterative Methods for Sparse Linear Systems.  SIAM, 2003.
  * [153] A. Greenbaum.  Iterative Methods for Solving Linear Systems.  SIAM, 1997.
  * [154] R. Barrett, M.W. Berry, T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. van der Vorst.  Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods.  SIAM, 1994.
  * [155] E. Béchet, H. Minnebo, N. Moës, and B. Burgardt.  Improved implementation and robustness study of the X-FEM for stress analysis around cracks.  International Journal for Numerical Methods in Engineering, 64(8):1033–1056, 2005.
  * [156] A. Menk and S.P.A. Bordas.  A robust preconditioning technique for the extended finite element method.  International Journal for Numerical Methods in Engineering, 85(13):1609–1632, 2011.
  * [157] B. Hiriyur, R.S. Tuminaro, H. Waisman, E.G. Boman, and D.E. Keyes.  A quasi-algebraic multigrid approach to fracture problems based on extended finite elements.  SIAM Journal on Scientific Computing, 34(2):A603–A626, 2012.
  * [158] L. Berger-Vergiat, H. Waisman, B. Hiriyur, R.S. Tuminaro, and D. Keyes.  Inexact Schwarz-algebraic multigrid preconditioners for crack problems modeled by extended finite element methods.  International Journal for Numerical Methods in Engineering, 90(3):311–328, 2012.
  * [159] H. Waisman and L. Berger-Vergiat.  An adaptive domain decomposition preconditioner for crack propagation problems modeled by XFEM.  International Journal for Multiscale Computational Engineering, 11(6):633–654, 2013.
  * [160] C. Lang, D. Makhija, and K. Doostan, A. Maute.  A simple and efficient preconditioning scheme for Heaviside enriched XFEM.  Computational Mechanics, 54:1357–1374, 2014.
  * [161] C. Lehrenfeld and A. Reusken.  Optimal preconditioners for Nitsche-XFEM discretizations of interface problems.  Numerische Mathematik, 135:313–332, 2017.
  * [162] T. Ludescher, S. Gross, and A. Reusken.  A multigrid method for unfitted finite element discretizations of elliptic interface problems.  SIAM Journal on Scientific Computing, 42(1):A318–A342, 2020.
  * [163] F. de Prenter, C.V. Verhoosel, E.H. van Brummelen, J.A. Evans, C. Messe, J.D. Benzaken, and K. Maute.  Multigrid solvers for immersed finite element methods and immersed isogeometric analysis.  Computational Mechanics, 65:807–838, 2020.
  * [164] S. Saberi, A. Vogel, and G. Meschke.  Parallel Finite Cell Method with Adaptive Geometric Multigrid.  In Euro-Par 2020: Parallel Processing, pages 578–593, 2020.
  * [165] J. Jomo, O. Oztoprak, F. de Prenter, N. Zander, S. Kollmannsberger, and E. Rank.  Hierarchical multigrid approaches for the finite cell method on uniform and multi-level _hp_ -refined grids.  Computer Methods in Applied Mechanics and Engineering, 386:114075, 2021.
  * [166] S. Badia and F. Verdugo.  Robust and scalable domain decomposition solvers for unfitted finite element methods.  Journal of Computational and Applied Mathematics, 344:740–759, 2018\.
  * [167] S. Badia, A.F. Martín, and J. Principe.  Multilevel balancing domain decomposition at extreme scales.  SIAM Journal on Scientific Computing, 38(1):C22–C52, 2016.
  * [168] H Kothari and R. Krause.  A Multigrid Method for a Nitsche-based Extended Finite Element Method.  International Journal of Computing and Visualization in Science and Engineering, 2021.
  * [169] A. Nüssing.  Fitted and unfitted finite element methods for solving the EEG forward problem.  PhD thesis, University of Münster, 2018.
  * [170] H. Kothari.  Multilevel Solution Strategies for Unfitted Finite Element Methods.  PhD thesis, Università della Svizzera italiana, 2020.
  * [171] H. Kothari and R. Krause.  Multigrid and saddle-point preconditioners for unfitted finite element modelling of inclusions.  In Proceedings of the 14th World Congress on Computational Mechanics (WCCM) European Community on Computational Methods in Applied Sciences (ECCOMAS) Congress 2020, 2021.
  * [172] H. Kothari and R. Krause.  A generalized multigrid method for solving contact problems in Lagrange multiplier based unfitted Finite Element Method.  Computer Methods in Applied Mechanics and Engineering, 392:114630, 2022.
  * [173] S. Gross and A. Reusken.  Optimal preconditioners for a Nitsche stabilized fictitious domain finite element method.  arXiv preprints, 2107.01182 [math.NA], 2021.
  * [174] S. Gross and A. Reusken.  Analysis of optimal preconditioners for CutFEM.  arXiv preprints, 2202.09069 [math.NA], 2022.
  * [175] B. Ayuso de Dios, K. Dunn, M. Sarkis, and S. Scacchi.  Dirichlet–Neumann Preconditioning for Stabilised Unfitted Discretization of High Contrast Problems.  In Domain Decomposition Methods in Science and Engineering XXV, pages 67–74, 2020.
  * [176] A.M. Matsokin and S.V. Nepomnyaschikh.  The Schwarz alternation method in a subspace.  Russian Mathematics (Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika), 29(10):78–84, 1985.
  * [177] P.L. Lions.  On the Schwarz Alternating Method. I.  In Proceedings of the First International Symposium on Domain Decomposition Methods for partial differential equations, pages 1–42, 1988.
  * [178] B.F. Smith, P.E. Bjørstad, and W.D. Gropp.  Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations.  Cambridge University Press, 1996.
  * [179] A. Toselli and O.B. Widlund.  Domain Decomposition Methods: Algorithms and Theory.  Springer, 2005.
  * [180] M. Dauge, A. Düster, and E. Rank.  Theoretical and Numerical Investigation of the Finite Cell Method.  Journal of Scientific Computing, 65:1039–1064, 2015.
  * [181] A. Buffa, R. Puppi, and R. Vazquez.  A minimal stabilization procedure for isogeometric methods on trimmed geometries.  SIAM Journal on Numerical Analysis, 58(5):2711–2735, 2020.
  * [182] P. Antolin, A. Buffa, R. Puppi, and X. Wei.  Overlapping multipatch isogeometric method with minimal stabilization.  SIAM Journal on Scientific Computing, 43(1):A330–A354, 2021.
  * [183] R. Puppi.  Isogeometric discretizations of the Stokes problem on trimmed geometries.  arXiv preprints, 2012.15582 [math.NA], 2021.
  * [184] J.M. Navarro-Jiménez, E. Nadal, M. Tur, J. Martìnez-Casas, and J.J. Ródenas.  On the use of stabilization techniques in the Cartesian grid finite element method framework for iterative solvers.  International Journal for Numerical Methods in Engineering, 121(13):3004–3020, 2020.
  * [185] D. Elfverson, M.G. Larson, and K. Larsson.  A new least squares stabilized Nitsche method for cut isogeometric analysis.  Computer Methods in Applied Mechanics and Engineering, 349:1–16, 2019.
  * [186] L. Larsson, S. Kollmannsberger, E. Rank, and M.G. Larson.  The finite cell method with least squares stabilized Nitsche boundary conditions.  Computer Methods in Applied Mechanics and Engineering, 393:114792, 2022.
  * [187] J.D. Hunt.  An adaptive three-dimensional Cartesian approach for the parallel computation of inviscid flow about static and dynamic configurations.  PhD thesis, University of Michigan, 2004.
  * [188] B. Müller, F. Kummer, and M. Oberlack.  Highly accurate surface and volume integration on implicit domains by means of moment-fitting.  International Journal for Numerical Methods in Engineering, 96(8):512–528, 2013.
  * [189] S. Badia, E. Neiva, and F. Verdugo.  Robust high-order unfitted finite elements by interpolation-based discrete extension.  arXiv preprints, 2201.06632 [math.NA], 2022.
  * [190] B. Marussig, J. Zechner, G. Beer, and T.-P. Fries.  Stable isogeometric analysis of trimmed geometries.  Computer Methods in Applied Mechanics and Engineering, 316:497–521, 2017.
  * [191] B. Marussig, R. Hiemstra, and T.J.R. Hughes.  Improved conditioning of isogeometric analysis matrices for trimmed geometries.  Computer Methods in Applied Mechanics and Engineering, 334:79–110, 2018.
  * [192] E. Burman, P. Hansbo, M.G. Larson, and K. Larsson.  Extension Operators for Trimmed Spline Spaces.  arXiv preprints, 2205.06543 [math.NA], 2022.
  * [193] S. Badia, H. Dilip, and F. Verdugo.  Space–time aggregated finite element methods for time-dependent problems on moving domains.  arXiv preprints, 2206.03626 [math.NA], 2022.
  * [194] F. Verdugo, A.F. Martín, and S. Badia.  Distributed-memory parallelization of the aggregated unfitted finite element method.  Computer Methods in Applied Mechanics and Engineering, 357:112583, 2019. Figures reproduced with permission from Elsevier.
  * [195] S. Badia, A.F. Martín, E. Neiva, and F. Verdugo.  The aggregated unfitted finite element method on parallel tree-based adaptive meshes.  SIAM Journal on Scientific Computing, 43(3):C203–C234, 2021.
  * [196] S. Badia, M.A. Caicedo, A.F. Martín, and J. Principe.  A robust and scalable unfitted adaptive finite element framework for nonlinear solid mechanics.  Computer Methods in Applied Mechanics and Engineering, 386:114093, 2021.
  * [197] E. Neiva and S. Badia.  Robust and scalable h-adaptive aggregated unfitted finite elements for interface elliptic problems.  Computer Methods in Applied Mechanics and Engineering, 380:113769, 2021.
  * [198] S. Badia, A.F. Martín, and F. Verdugo.  Mixed aggregated finite element methods for the unfitted discretization of the Stokes problem.  SIAM Journal on Scientific Computing, 40(6):B1541–B1576, 2018\.
  * [199] R. Codina, S. Badia, J. Baiges, and J. Principe.  Variational Multiscale Methods in Computational Fluid Dynamics.  Encyclopedia of Computational Mechanics Second Edition, pages 1–28, 2017.
  * [200] S. Badia, A.F. Martín, and J. Principe.  FEMPAR: An Object-Oriented Parallel Finite Element Framework.  Archives of Computational Methods in Engineering, 25:195–271, 2018\.
  * [201] S. Balay, S. Abhyankar, M.F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley, D.A. May, L.C. McInnes, R.T. Mills, T. Munson, K. Rupp, P. Sanan, B.F. Smith, S. Zampini, H. Zhang, and H. Zhang.  PETSc Users Manual.  Technical Report ANL-95/11 - Revision 3.12, Argonne National Laboratory, 2019.
  * [202] D.A. May, P. Sanan, K. Rupp, M.G. Knepley, and B.F. Smith.  Extreme-scale Multigrid Components within PETSc.  In PASC ’16: Proceedings of the Platform for Advanced Scientific Computing Conference, number 5, 2016.
  * [203] C. Burstedde, L.C. Wilcox, and O. Ghattas.  p4est: scalable algorithms for parallel adaptive mesh refinement on forests of octrees.  SIAM Journal on Scientific Computing, 33(3):1103–1133, 2011.
  * [204] W.G. Dettmer, C. Kadapa, and D. Perić.  A stabilised immersed boundary method on hierarchical b-spline grids.  Computer Methods in Applied Mechanics and Engineering, 311:415–437, 2016.
  * [205] P. Hansbo, M.G. Larson, and K. Larsson.  Cut Finite Element Methods for Linear Elasticity Problems.  In Geometrically Unfitted Finite Element Methods and Applications: Proceedings of the UCL Workshop 2016, pages 578–593, 2017.
  * [206] E. Burman, P. Hansbo, and M.G. Larson.  On the Design of Locking Free Ghost Penalty Stabilization and the Relation to CutFEM with Discrete Extension.  arXiv preprints, 2205.01340 [math.NA], 2022.
  * [207] J. Preuss.  Higher order unfitted isoparametric space–time FEM on moving domains.  Master’s thesis, University of Göttingen, 2018.
  * [208] C. Lehrenfeld and M. Olshanskii.  An Eulerian finite element method for PDEs in time-dependent domains.  ESAIM: Mathematical Modelling and Numerical Analysis (ESAIM: M2AN), 53(2):585–614, 2019.
  * [209] M.G. Larson and S. Zahedi.  Conservative Discontinuous Cut Finite Element Methods.  arXiv preprints, 2105.02202 [math.NA], 2021.
  * [210] R. Becker and M. Braack.  A finite element pressure gradient stabilization for the Stokes equations based on local projections.  CALCOLO, 38:173–199, 2001.
  * [211] S. Badia.  On stabilized finite element methods based on the Scott–Zhang projector. Circumventing the inf–sup condition for the Stokes problem.  Computer Methods in Applied Mechanics and Engineering, 247-248:65–72, 2012.
  * [212] A. Johansson, M.G. Larson, and A. Logg.  A Multimesh Finite Element Method for the Stokes Problem.  In Numerical Methods for Flows: FEF 2017 Selected Contributions, pages 189–198, 2020.
  * [213] A. Johansson, B. Kehlet, M.G. Larson, and A. Logg.  Multimesh finite element methods: Solving PDEs on multiple intersecting meshes.  Computer Methods in Applied Mechanics and Engineering, 343:672–689, 2019.
  * [214] A. Johansson, M.G. Larson, and A. Logg.  Multimesh finite elements with flexible mesh sizes.  Computer Methods in Applied Mechanics and Engineering, 372:113420, 2020.
  * [215] L. Nguyen, S.K.F. Stoter, T. Baum, J. Kirschke, M. Ruess, Z. Yosibash, and D. Schillinger.  Phase-field boundary conditions for the voxel finite cell method: Surface-free stress analysis of CT-based bone structures.  International Journal for Numerical Methods in Biomedical Engineering, 33(12):e2880, 2017.
  * [216] T. Hoang, C.V. Verhoosel, F. Auricchio, E.H. van Brummelen, and A. Reali.  Skeleton-stabilized IsoGeometric Analysis: High-regularity interior-penalty methods for incompressible viscous flow problems.  Computer Methods in Applied Mechanics and Engineering, 337:324–351, 2018.
  * [217] J. Guzmán and M. Olshanskii.  Inf–sup stability of geometrically unfitted Stokes finite elements.  Mathematics of Computation, 87:2091–2112, 2018.
  * [218] P. Hansbo, M.G. Larson, and S. Zahedi.  Characteristic cut finite element methods for convection–diffusion problems on time dependent surfaces.  Computer Methods in Applied Mechanics and Engineering, 293:431–461, 2015.
  * [219] E. Burman, P. Hansbo, and M.G. Larson.  A stabilized cut finite element method for partial differential equations on surfaces: The Laplace–Beltrami operator.  Computer Methods in Applied Mechanics and Engineering, 285:188–207, 2015.
  * [220] E. Burman, P. Hansbo, M.G. Larson, and A. Massing.  A cut discontinuous Galerkin method for the Laplace–Beltrami operator.  IMA Journal of Numerical Analysis, 37(1):138–169, 2017.
  * [221] P. Clément.  Approximation by finite element functions using local regularization.  Revue française d’automatique, informatique, recherche opérationnelle. Analyse numérique, 9(R2):77–84, 1975.
  * [222] E. Burman, P. Hansbo, and M.G. Larson.  CutFEM Based on Extended Finite Element Spaces.  arXiv preprints, 2101.10052 [math.NA], 2022.
  * [223] S.C. Brenner and L.R. Scott.  The Mathematical Theory of Finite Element Methods.  Springer, 2008.
  * [224] N. Zander, T. Bog, S. Kollmannsberger, D. Schillinger, and E. Rank.  Multi-level _hp_ -adaptivity: high-order mesh adaptivity without the difficulties of constraining hanging nodes.  Computational Mechanics, 55:499–517, 2015.
  * [225] D. D’Angella, N. Zander, S. Kollmannsberger, F. Frischmann, E. Rank, A. Schröder, and A. Reali.  Multi-level _hp_ -adaptivity and explicit error estimation.  Advanced Modeling and Simulation in Engineering Sciences, 3(33), 2016.
  * [226] S.C. Divi, P.H. van Zuijlen, T. Hoang, F. de Prenter, F. Auricchio, A. Reali, E.H. van Brummelen, and C.V. Verhoosel.  Residual-based error estimation and adaptivity for stabilized immersed isogeometric analysis using truncated hierarchical B-splines.  Jounal of Mechanics, 38:204–237, 2022.
  * [227] A. Gerstenberger and W. A. Wall.  An embedded Dirichlet formulation for 3D continua.  International Journal for Numerical Methods in Engineering, 82(5):537–563, 2010.
  * [228] C. Wolfe.  The immersed boundary approach to fluid flow simulation, 2009.  _https://www.ansys.com/resource-center/article/immersed-boundary-approach-to-fluid-flow-simulation-ansys-advantage-v3-i2_ , assessed on 15-5-2022.
  * [229] S. Hartmann, L.F. Leidinger, D.J. Benson, L. Li, A.P. Nagy, and M. Pigazzini.  Enabling the Analysis of Topologically Connected Multi-Patch Trimmed NURBS Shells in LS-DYNA.  In Proceedings of the 12th European LS-DYNA Conference, 2019.
  * [230] M. Meßmer, L.F. Leidinger, S. Hartmann, F. Bauer, F. Duddeck, R. Wüchner, and K.-U. Bletzinger.  Isogeometric Analysis on Trimmed Solids: A B-Spline-Based Approach Focusing on Explicit Dynamics.  In Proceedings of the 13th European LS-DYNA Conference, 2021.
  * [231] A.J. Favaloro, B. Brenken, E. Barocio, and R.B. Pipes.  Simulation of Polymeric Composites Additive Manufacturing using Abaqus.  In Proceedings of Dassault Systèmes’ Science in the Age of Experience, pages 103–114, 2017.
  * [232] B. Courter, V. Savane, J. Bi, S. Dev, and C.J. Hansen.  Finite Element Simulation of the Fused Deposition Modelling Process.  In Proceedings of the NAFEMS World Congress, pages 11–14, 2017\.
  * [233] Hyperganic.  _https://www.hyperganic.com_ , assessed on 15-5-2022.
  * [234] W. Bangerth, R. Hartmann, and G. Kanschat.  deal.II–A general-purpose object-oriented finite element library.  ACM Transactions on Mathematical Software, 33(4):Article 24, 2007\.
  * [235] M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M.E. Rognes, and G.N. Wells.  The FEniCS Project Version 1.5.  Archive of Numerical Software, 3(100):9–23, 2015.
  * [236] S. Badia and F. Verdugo.  Gridap: An extensible Finite Element toolbox in Julia.  Journal of Open Source Software, 5(52):2520, 2020.
  * [237] J. Schöberl.  NETGEN An advancing front 2D/3D-mesh generator based on abstract rules.  Computing and Visualization in Science, 1:41–52, 1997.
  * [238] P. Bastian, M. Blatt, A. Dedner, N.-A. Dreier, C. Engwer, R. Fritze, C. Gräser, C. Grüninger, D. Kempf, R. Klöfkorn, M. Ohlberger, and O. Sander.  The DUNE framework: Basic concepts and recent developments.  Computers & Mathematics with Applications, 81:75–112, 2021.