#  Robust and parallel scalable iterative solutions
for large-scale finite cell analyses

J. N. Jomo Corresponding author<EMAIL_ADDRESS>Chair for Computation in
Engineering, Technische Universität München  F. de Prenter  Department of
Mechanical Engineering, Eindhoven University of Technology  M. Elhaddad  Chair
for Computation in Engineering, Technische Universität München  D. D’Angella
Chair for Computation in Engineering, Technische Universität München
Institute for Advanced Study, Technische Universität München
C. V. Verhoosel  Department of Mechanical Engineering, Eindhoven University of
Technology  S. Kollmannsberger  Chair for Computation in Engineering,
Technische Universität München  J. S. Kirschke  Abteilung für Neuroradiologie,
Klinikum rechts der Isar, Technische Universität München
V. Nübel  Hilti Entwicklungsgesellschaft GmbH  E. H. van Brummelen  Department
of Mechanical Engineering, Eindhoven University of Technology  E. Rank  Chair
for Computation in Engineering, Technische Universität München  Institute for
Advanced Study, Technische Universität München

## Abstract

The finite cell method is a highly flexible discretization technique for
numerical analysis on domains with complex geometries. By using a non-boundary
conforming computational domain that can be easily meshed, automatized
computations on a wide range of geometrical models can be performed.
Application of the finite cell method, and other immersed methods, to large
real-life and industrial problems is often limited due to the conditioning
problems associated with these methods. These conditioning problems have
caused researchers to resort to direct solution methods, which significantly
limit the maximum size of solvable systems. Iterative solvers are better
suited for large-scale computations than their direct counterparts due to
their lower memory requirements and suitability for parallel computing. These
benefits can, however, only be exploited when systems are properly
conditioned. In this contribution we present an Additive-Schwarz type
preconditioner that enables efficient and parallel scalable iterative
solutions of large-scale multi-level $hp$-refined finite cell analyses.

Keywords: immersed methods, finite cell method, iterative solvers,
preconditioning, _hp_ -refinement, parallel computing, high performance
computing

##  1 Introduction

Immersed finite element methods, most prominently the finite cell method [1,
2, 3, 4, 5, 6] and cutFEM [7, 8], are useful techniques for performing
numerical computations on domains with complex geometries. Generating boundary
fitted meshes for these problems may be an error-prone, laborious and
computationally expensive operation. The general concept of immersed finite
element methods is not to directly mesh the complex geometry on which the
problem is posed, but to instead embed it in a larger, embedding domain of
simple shape that can be trivially discretized. An approximation space is then
defined on this structured mesh and restricted to the complex geometry of the
problem domain.

The finite cell method (FCM) combines the embedded mesh concept with high-
order finite elements to compute highly accurate results with a flexible
discretization scheme. The efficiency of immersed methods can be further
improved by using them in conjunction with $hp$-refined discretizations. The
multi-level $hp$-finite cell method [9, 10] has shown to be effective in
capturing local solution characteristics, such as steep gradients or stress
concentrations induced by (fine) geometrical features, e.g. [11, 12, 13]. An
unavoidable aspect of applying this method to large-scale real-life
applications is the necessity of solving large linear systems. The
computational cost of using a direct solver does not scale well with the size
of the system, however, and becomes increasingly expensive in terms of memory
consumption and execution time. Iterative solvers are the preferred method for
solving large systems, because their computational cost scales better with the
system size and due to their low memory requirement and ability to be easily
parallelized in comparison to direct solvers. Convergence of iterative
solvers, however, strongly depends on the conditioning of the system, e.g.
[14]. Without tailored preconditioning, FCM systems generally show severe ill-
conditioning, which practically prohibits application of iterative solvers
[15]. As a result, finite cell computations have been mainly restricted to the
use of direct solvers [2, 16, 17, 18, 19, 3].

Several techniques have been developed in order to improve the conditioning of
immersed finite elements and enable iterative solution techniques. Fictitious
domain stiffness assumes a soft material in the fictitious domain, and thereby
performs a volumetric stabilization through the evaluation of integrals in
both the problem and fictitious domain. An extensive analysis of this method
can be found in [20]. For FCM problems in the context of isogeometric
analysis, it is shown in [21, 22] that problematic basis functions with small
support inside the domain of computation can simply be removed from the
system. Ghost penalty stabilization [23] is another option which weakly
couples the degrees of freedom (DOFs) in cut elements to their neighboring
elements. This gives better control of the solution on cut elements and
improves the conditioning. Application of ghost penalty is customary in
methods referred to as cutFEM [7, 8]. Instead of weakly coupling DOFs in cut
elements to neighboring elements, these can also be constrained strongly as
originally introduced in [24, 25] and applied in e.g. [26, 27, 28, 29, 30,
31].

In addition to these stabilization methods multiple preconditioning techniques
have been proposed for immersed finite element methods and for similar
problems in the extended finite element method (X-FEM) [32, 33]. In X-FEM,
these preconditioners generally split the matrix in _(i)_ non-problematic
DOFs, which can be treated by standard preconditioning techniques and _(ii)_
problematic DOFs, that are treated separately by local Cholesky decompositions
[34], a tailored FETI-type method [35], a Schur complement based algebraic
multigrid preconditioner [36] or a Schwarz-type domain decomposition
preconditioner [37, 38]. The last two references are conceptually similar to
the approach applied in this work, but differ in the choice of subdomains as
the problematic DOFs in the finite cell method, i.e. all DOFs on the boundary,
are generally large connected sets. For immersed finite element methods it is
proposed in [39] to diagonally scale the problematic DOFs and treat the non-
problematic DOFs by an algebraic multigrid procedure, which is demonstrated to
be effective for linear basis functions and restrictions on the manner in
which elements can be cut. In [40] the scaling in a Balancing Domain
Decomposition by Constraints preconditioner (BDDC) is customized for cut basis
functions, which is shown to result in an effective method for linear bases. A
preconditioner that combines a diagonal scaling with a local
orthonormalization of the problematic DOFs is developed in [15]. In [41] it is
shown that this orthonormalization is very similar to Additive-Schwarz
preconditioning, e.g. [42], of the cut elements. The preconditioner developed
in [41] is demonstrated to be effective for immersed finite element methods
with higher order basis functions, but has so far only been tested on uniform
tensor product meshes without local refinements.

In this contribution, we further develop the preconditioning technique
presented in [41] through modifications that are essential to allow robust
application of iterative solvers for finite cell discretizations involving
local $hp$-refinements. The robustness of the preconditioner is improved by
stabilizing the local inverses in case the smallest eigenvalues approach
machine precision. Furthermore, the efficiency and scalability of the
preconditioner setup is investigated and optimized by (i) focusing on
problematic areas and (ii) the use of shared memory and distributed memory
parallelization. The developed technique is demonstrated to be able to
efficiently solve large FCM systems involving multi-level $hp$-refinement with
good parallel scalability using a preconditioned Conjugate Gradient (PCG)
solver [14]. The presented examples show the potential of using FCM in large-
scale real-life and industrial applications.

Section 2 introduces the finite cell method and multi-level $hp$-refinement.
Section 3 recapitulates the conditioning of the finite cell method and the
principal concepts of Additive-Schwarz preconditioning, which is followed by
the modifications required for preconditioning multi-level $hp$-refined
discretizations. Section 4 presents the details of the parallel implementation
and Section 5 investigates and demonstrates the effectiveness of the
preconditioner with a PCG solver in several numerical examples.

##  2 Discretization techniques

In this section we introduce the discretization methods that form the basis of
the proposed techniques, namely the finite cell method in Section 2.1 and
multi-level $hp$-refinement in Section 2.2 .

###  2.1 The finite cell method

The finite cell method (FCM) is an immersed high-order method suitable for
numerical analysis on domains of complex shape [1, 2]. As an underlying
numerical approach finite elements or spline-based shape functions can be
used. In the latter case FCM turns out to be a suitable technique for trimmed
isogeometric analysis [43, 44]. FCM extends the original problem domain
referred to as the physical domain, $\Omega_{\textrm{phys}}$, with a
surrounding fictitious domain, $\Omega_{\textrm{fict}}$, forming a
computational domain $\Omega$ with a simple structure that can be trivially
meshed (Figure 1). Distinction between points in the physical and fictitious
domains is made by an indicator function $\alpha$, which assumes the value
$\alpha=1$ for points lying in $\Omega_{\textrm{phys}}$ and $0\leq\alpha\ll 1$
for points outside $\Omega_{\textrm{phys}}$. The parameter $\alpha$ is
generally chosen as a small finite constant $0<\alpha\ll 1$ instead of
$\alpha=0$ to ensure numerical stability, e.g. [20].

(a) Physical domain $\Omega_{\text{phys}}$ (b) Fictitious extension:
$\Omega_{\textrm{phys}}\cup\Omega_{\textrm{fict}}$ (c) Structured
discretization

Figure 1: A schematic representation of the finite cell method [12].

Since the basis is generally not interpolatory on the physical domain
boundary, Dirichlet conditions are weakly imposed using e.g. Nitsche’s method
[45, 46] or the penalty method [47]. Because FCM results in conditioning
problems of a similar nature with both Nitsche’s method and the penalty
method, we apply the penalty method in all our examples because of its
simplicity and low computational cost. All presented results can be extended
to problems with Nitsche boundary conditions without significant changes.

The weak form in FCM is illustrated for a linear elasticity problem in Eq. (1)
where $\boldsymbol{u}$ is the unknown displacement field, $\boldsymbol{v}$ an
arbitrary test function, $\boldsymbol{\varepsilon}$ the infinitesimal strain
tensor, $\mathbb{C}$ the fourth-order elasticity tensor and $\beta$ the
penalty parameter. The vectors $\boldsymbol{b}$, $\boldsymbol{f}$ and
$\boldsymbol{g}$ represent the volumetric forces, surface traction and
Dirichlet data respectively.

$\int\limits_{\Omega}\alpha\,\boldsymbol{\varepsilon}(\boldsymbol{v}):\mathbb{C}:\boldsymbol{\varepsilon}(\boldsymbol{u})\,d\Omega\
+\
\int\limits_{\Gamma_{D}}\beta\,\boldsymbol{v}\cdot\boldsymbol{u}\,d\Gamma=\int\limits_{\Omega}\alpha\,\boldsymbol{v}\cdot\boldsymbol{b}\,d\Omega\
+\ \int\limits_{\Gamma_{N}}\boldsymbol{v}\cdot\boldsymbol{f}\,d\Gamma\
+\int\limits_{\Gamma_{D}}\beta\,\boldsymbol{v}\cdot\boldsymbol{g}\,d\Gamma\,.$
(1)

FCM has a large application field due to its ability to produce highly
flexible, non-conforming discretizations while maintaining the accuracy and
computational efficiency of high-order finite elements. For more details
regarding FCM the reader is referred [1, 2] or the review articles [3, 6] that
show various applications and recent developments.

###  2.2 Multi-level $hp$-discretizations

For problems with non-smooth solutions, a high quality finite element
approximation can be obtained with a relatively small number of unknowns by
the use of $hp$-refinement schemes. These schemes combine a local increase of
the spatial resolution ($h$-refinement) and an elevation of the element
polynomial order ($p$-refinement) in the crucial regions of the computational
domain. Often $hp$-refinement schemes perform spatial refinement by the
replacement of large elements having a high approximation error with a set of
smaller elements, e.g. [48, 49, 50]. Multi-level $hp$-refinement, introduced
in [9], performs $hp$-refinement by means of high-order hierarchical overlay
meshes. Instead of replacing large coarse elements, multi-level
$hp$-refinement superposes coarse elements with a hierarchy of overlay
elements in order to better capture solution characteristics in regions of
interest. The solution is the sum of the large-scale solution on the coarse
base mesh and the fine-scale solutions on the overlay meshes. This change of
paradigm from a replacement-driven to a superposition-driven spatial
refinement significantly reduces the complexity associated with the
implementation of $hp$-refinement schemes. The direct relation of degrees of
freedom to topological components (nodes, edges, faces and volumes) in high-
order finite elements allows a straightforward enforcement of both the
compatibility of the discretization and linear independence of the basis
functions in the multi-level $hp$-refinement setting c.f [9, 10] and [11, 13]
for applications.

(a) One-dimensional case (b) Two-dimensional case

problemformulation/pics/multiLevelhp3d-notext/_[rgb(0pt)=(0,0,0);
rgb(21.25bp)=(0,0,0); rgb(23.125bp)=(1,1,1); rgb(25bp)=(0,0,0);
rgb(50bp)=(0,0,0)]

(c) Three-dimensional case

Figure 2: Illustration of the multi-level $hp$-refinement scheme with two
refinement levels, $k=2$, in different spatial dimensions. The deactivation of
specific topological components following a simple rule-set ensures
compatibility and linear independence of the basis functions [10].

Figure 2 illustrates the superposition principle in multi-level
$hp$-refinement. The scheme is not restricted to 1-irregular meshes since
arbitrary levels of hanging nodes are treated automatically through the use of
homogeneous Dirichlet boundary conditions at the boundary of the overlay
meshes. This deactivates specific topological components and enforces
$C^{0}$-continuity over the complete refinement hierarchy.
$C^{\infty}$-continuity is maintained in unrefined base elements and overlay
elements on the highest refinement level. Such elements are collectively
referred to as _active or leaf elements_. Numerical integration is performed
on these leaf elements.

##  3 Preconditioning

As mentioned in Section 1 , iterative solution methods are preferred over
direct solvers for large systems due to their lower memory requirements and
suitability for parallel computing. It is well known that the convergence
speed of iterative methods generally depends on the conditioning of the
system, e.g. [14]. The conditioning of FCM systems is generally rather poor,
as remarked in e.g. [16, 17, 18, 19, 3] and shown in detail in [15], which
impedes solving large FCM systems. In [41] an Additive-Schwarz scheme is
proposed that is efficient for FCM systems on uniform grids, but that does not
account for the local refinements in multi-level $hp$-discretizations. Section
3.1 revisits the conditioning problems associated with the finite cell method
and discusses how Additive-Schwarz preconditioners can alleviate these
problems. Section 3.2 presents the required modifications to the Additive-
Schwarz scheme that allow efficient preconditioning of immersed multi-level
$hp$-discretizations.

###  3.1 FCM conditioning and Additive-Schwarz preconditioning

This manuscript focuses on problems in linear elasticity, for which FCM yields
a symmetric positive definite (SPD) system of the form

$\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}\,,$ (2)

with system matrix $\boldsymbol{A}$, solution vector $\boldsymbol{x}$ and
right hand side vector $\boldsymbol{b}$. The condition number
$\kappa\left(\boldsymbol{A}\right)$ of SPD systems is defined as

$\kappa\left(\boldsymbol{A}\right)=\large\|\boldsymbol{A}\large\|\,\large\|\boldsymbol{A}^{-1}\large\|=\frac{\lambda_{\textsc{max}}}{\lambda_{\textsc{min}}}=\frac{\max_{\|\boldsymbol{v}\|=1}\boldsymbol{v}^{T}\boldsymbol{A}\boldsymbol{v}}{\min_{\|\boldsymbol{u}\|=1}\boldsymbol{u}^{T}\boldsymbol{A}\boldsymbol{u}}\,,$
(3)

with $\lambda_{\textsc{max}}$ and $\lambda_{\textsc{min}}$ denoting the
largest and smallest eigenvalue and $\boldsymbol{v}$ and $\boldsymbol{u}$
(eigen)vectors maximizing and minimizing the Rayleigh quotient with
$\boldsymbol{A}$. From [15] we know that FCM systems are generally ill-
conditioned because of the occurrence of elements with a very small
intersection with the physical domain. Basis functions on such elements can be
arbitrarily small in magnitude and/or basis functions that are sufficiently
linearly independent on the full element can become almost linearly dependent
when the element is cut. The latter is caused firstly by higher order
contributions of basis functions, that can become arbitrarily small compared
to linear (or lower order) contributions, and secondly by dependencies on a
(e.g. horizontal or vertical) parametric dimension, that can become
arbitrarily small compared to other parametric dimensions. This reduces a
higher order or multivariate function to an almost linear or univariate
function. When this occurs for multiple functions, these functions can become
almost linearly dependent. The magnitude of a function in the approximation
space that is only supported on a small cut element — which strongly depends
on the volume of the intersection between the element and the physical domain
— can therefore be arbitrarily small compared to the magnitude of the vector
that represents the function in the isomorphic vector space — which is
independent of this volume. Because of that, FCM systems generally have
$\min_{\|\boldsymbol{u}\|=1}\boldsymbol{u}^{T}\boldsymbol{A}\boldsymbol{u}\ll
1$, such that $\kappa\left(\boldsymbol{A}\right)\gg 1$. When the smallest cut
element is assumed to be shape regular and for $\alpha=0$, it is shown in [15]
that the condition number of second order FCM systems can be estimated by

$\kappa\left(\boldsymbol{A}\right)\propto C\eta^{-(2p+1-2/d)}\,,$ (4)

with $p$ the polynomial order of the discretization, $d$ the number of
dimensions and $\eta$ the relative volume of the intersection between the
smallest cut element and the physical domain

$\eta=\min_{T\in\mathcal{T}}\frac{|T\cap\Omega_{\text{phys}}|}{|T|}\,,$ (5)

with $\mathcal{T}$ the set of elements, $|T\cap\Omega_{\text{phys}}|$ the
volume of the intersection between element $T$ and the physical domain and
$|T|$ the full volume of element $T$. Since the volume of the intersection can
become arbitrarily small, the condition number in FCM can be arbitrarily
large.

###### Remark 1

It is stipulated that the ill-conditioning of FCM is not only due to basis
functions that become very small. This could easily be resolved by Jacobi or
Gauss-Seidel preconditioning. Notwithstanding that Jacobi and Gauss-Seidel
preconditioners improve the conditioning, these do not repair conditioning
problems related to almost linearly dependent basis functions on small cut
elements. Therefore these simple preconditioners do not result in robust
preconditioning for the finite cell method, as is shown in [15] and
demonstrated further in Section 5.2 .

#### Additive-Schwarz preconditioning

Ill-conditioning in FCM can be effectively resolved by the use of Additive-
Schwarz preconditioners [41]. A rich literature exists on these
preconditioners in a general setting, e.g. [51, 52], and specifically for
finite element methods, e.g. [53, 42]. The general idea of Additive-Schwarz
preconditioners is the construction of a sparse approximation of
$\boldsymbol{A}^{-1}$ through the inversion and summation of sub-matrices of
matrix $\boldsymbol{A}$ as per the formula

$\boldsymbol{S}=\sum_{B\in\mathcal{B}}\boldsymbol{R}_{B}^{T}\underbrace{\left(\boldsymbol{R}_{B}\boldsymbol{A}\boldsymbol{R}_{B}^{T}\right)^{-1}}_{\boldsymbol{A}_{B}^{-1}}\boldsymbol{R}_{B}\,.$
(6)

In this formulation $B$ denotes a _block_ , which is a set of indices
referring to basis functions. The set of all blocks is denoted by
$\mathcal{B}$. A basis function $\phi_{k}$ can be in multiple (i.e.
overlapping) blocks, but each basis function must be in at least one block for
$\boldsymbol{S}$ to be nonsingular. Premultiplying matrix $\boldsymbol{A}$
with the (non-square) restriction matrix $\boldsymbol{R}_{B}$ restricts
$\boldsymbol{A}$ to the row indices in $B$, and postmultiplying matrix
$\boldsymbol{A}$ with the restriction matrix $\boldsymbol{R}_{B}$ restricts
$\boldsymbol{A}$ to the column indices in $B$. The inverse of sub-matrix
$\boldsymbol{A}_{B}$ is placed at the indices of block $B$ by pre- and
postmultiplying with $\boldsymbol{R}_{B}^{T}$ and $\boldsymbol{R}_{B}$
respectively.

Properly choosing the blocks $\mathcal{B}$ is essential for the effectiveness
of Additive-Schwarz preconditioning, because different blocks yield varying
preconditioners with different properties. For example one can assign
$\mathcal{B}$ such that every function is in a single, separate block. This
results in the standard Jacobi preconditioner which is not robust to cut
elements, as demonstrated in Section 5.2 . On the other extreme one can
theoretically assign only one block containing all functions such that
$\boldsymbol{S}=\boldsymbol{A}^{-1}$. Such a preconditioner is optimal in
terms of spectral properties, but is prohibitively expensive as it involves
inverting the full system. For FCM, it is useful to assign blocks based on the
Additive-Schwarz lemma [54, 55, 51, 52]:

$\boldsymbol{v}^{T}\boldsymbol{S}^{-1}\boldsymbol{v}=\min_{\boldsymbol{v}=\sum\limits_{B\in\mathcal{B}}\boldsymbol{R}_{B}^{T}\boldsymbol{v}_{B}}\sum_{B\in\mathcal{B}}\boldsymbol{v}_{B}^{T}\boldsymbol{A}_{B}\boldsymbol{v}_{B}\,.$
(7)

In this formulation $\boldsymbol{v}_{B}$ denotes a vector whose length
corresponds to the size of block $B$.
$\boldsymbol{R}_{B}^{T}\boldsymbol{v}_{B}$ is a prolongated vector whose
length corresponds to the size of the full system and which only has nonzero
entries at the indices of block $B$. When each index is in at least one block,
for every vector $\boldsymbol{v}$ there exists a set
$\\{\boldsymbol{v}_{B}\\}$ such that
$\boldsymbol{v}=\sum_{B\in\mathcal{B}}\boldsymbol{R}_{B}^{T}\boldsymbol{v}_{B}$.
With overlapping blocks, different sets $\\{\boldsymbol{v}_{B}\\}$ with this
property exist. The lemma states that the inner product of a vector
$\boldsymbol{v}$ with $\boldsymbol{S}^{-1}$ is equal to the minimum of the sum
of inner products over all sets $\\{\boldsymbol{v}_{B}\\}$ that sum up to
$\boldsymbol{v}$.

Next, recall the mechanism of small eigenvalues in FCM when basis functions on
small cut elements are small or almost linearly dependent. The Additive-
Schwarz lemma establishes that a vector $\boldsymbol{v}$ that has a small
inner product with $\boldsymbol{A}$ because of this mechanism, will also have
a small inner product with $\boldsymbol{S}^{-1}$ provided that functions that
are almost linearly dependent are aggregated in one block. This clustering can
be achieved by assigning one block for every element, consisting of all basis
functions that are supported on that element, since basis functions need
intersecting supports in order to be almost linearly dependent. When the
blocks are set up in this manner, it so follows from the Additive-Schwarz
lemma that $\boldsymbol{S}^{-1}$ inherits the small eigenvalues of
$\boldsymbol{A}$ caused by small cut elements. The preconditioner can
therefore effectively treat the fundamental cause of ill-conditioning of the
finite cell method, which has already been demonstrated for uniform
discretizations in [41].

###  3.2 Preconditioning the multi-level $hp$-refined finite cell method

It is not straightforward to apply this principle to multi-level $hp$-refined
grids. To begin with, the notion of an element is less trivial, because of the
different levels of elements. Considering the requirement that functions that
are almost linearly dependent need to be in the same block, as stated in
Section 3.1 , it is sufficient to only set blocks for the active or leaf
elements, i.e. elements that do not have overlay elements on a finer level.
This still does not yield optimal results however. Therefore the following
three modifications are made to the preconditioning technique.

#### Selecting suitable blocks

Uniform grids, as used in [41], and multi-level $hp$-grids, as applied here,
have two fundamental differences with regard to choosing Additive-Schwarz
blocks in an element-wise manner based on the functions supported on a leaf
element. First, on uniform grids the number of basis functions supported on an
element is constant, i.e. $(p+1)^{d}$ for scalar problems with $p$ denoting
the polynomial order and $d$ the number of dimensions 111 For the full tensor
product space, e.g [56]. . For multi-level $hp$-grids the total number of
basis functions supported on a leaf element is, in general, not bounded. This
is attributed to the superposed linear hat functions that are present on every
refinement level. This superposition results in multiply-defined linear basis
functions within a leaf element as shown in Figure 2 . Directly applying the
strategy from [41] to leaf elements in multi-level $hp$-grids therefore yields
large blocks that are computationally expensive in the setup, storage and
application of the preconditioner. Second, on uniform grids the number of
elements that a basis function is supported on is bounded, i.e. $\leq 2^{d}$
for $C^{0}$-finite elements. With multi-level $hp$-discretizations, it is
possible that a (coarse) linear basis function is supported on a large number
of (fine) leaf elements, Figure 2 . Consider the case that a basis function
$\phi_{k}$ is in $n$ different blocks with $n\gg 1$. The unit vector
$\boldsymbol{e}_{k}$ corresponding to $\phi_{k}$ has a Rayleigh quotient
$A_{kk}$ with $\boldsymbol{A}$. Next we construct a set
$\\{\boldsymbol{v}_{B}\\}$ with the value $1/n$ for all vector entries
corresponding to $\phi_{k}$ and the value $0$ for all other entries. Clearly
$\sum_{B\in\mathcal{B}}\boldsymbol{R}_{B}^{T}\boldsymbol{v}_{B}=\boldsymbol{e}_{k}$
and from Eq. (7) it follows that the Rayleigh quotient of $\boldsymbol{e}_{k}$
with $\boldsymbol{S}^{-1}$ is bounded from above by $A_{kk}/n$. This shows
that a basis function that is in many different blocks yields a small
eigenvalue in $\boldsymbol{S}^{-1}$, relative to $\boldsymbol{A}$, and
consequently a large eigenvalue in $\boldsymbol{S}$. This reduces the
efficiency of the preconditioner, which is demonstrated in Sections 5.1 and
5.3 .

In order to obtain a robust and efficient Additive-Schwarz preconditioner for
the multi-level $hp$-finite cell method, the blocks in $\mathcal{B}$ have to
be carefully selected to ensure that _(i)_ the number of blocks a given basis
function belongs to is bounded and that _(ii)_ the size of each block is
small. As already noted in [41] it is _sufficient_ but not _necessary_ to
devise a block for every element with _all_ basis functions supported on it.
In fact it is enough to overcome the conditioning problems of FCM, due to
small eigenvalues on small cut elements, if the blocks contain the basis
functions that can potentially become almost linearly dependent. This
obviously reduces the size of the blocks, and thereby repairs the large
computational cost of setting up, storing and applying the preconditioner.
Furthermore, the large linear basis functions that are supported on many leaf
elements, and cause the second aforementioned efficiency problem, are not the
basis functions that suffer from almost linear dependencies. Therefore,
truncating the blocks to only the basis functions that can become almost
linearly dependent when the element is cut resolves both mentioned issues
related to the preconditioning of the multi-level $hp$-finite cell method.

Figure 3: Block selection for multi-level $hp$-grids with the full and
truncated block of element $T_{3}$.

To describe the truncation, we first consider a one dimensional multi-level
$hp$-grid with elements of polynomial order $p=3$. Figure 3a illustrates that
more than $p+1$ basis functions can be supported on a leaf element due to
multiply-defined linear basis functions. As mentioned previously, these large
sets of supported basis functions need to be truncated in order to result in
efficient element-wise Additive-Schwarz blocks. To achieve this, it is first
noted that although more then $p+1$ basis functions can be supported on a
single element, each element only has $p+1$ unique polynomial degrees of
freedom. Since eigenvalues with problematically small eigenmodes are only
supported on a small cut element, efficient preconditioning only requires the
Additive-Schwarz block of an element to span these $p+1$ polynomial degrees of
freedom. This set of $p+1$ functions must be made up of the $p-1$ higher order
basis functions accompanied by two linear basis functions. These linear basis
functions depend on the position and refinement depth of a leaf element as
illustrated in Figure 3a . In the grid shown, two different cases are
distinguished: _(i)_ the element simply has two linear functions on the
highest level e.g. $T_{1}$ and $T_{5}$; _(ii)_ the element has one linear
function on the highest level and the second linear function is selected by
traversing down the element hierarchy and taking the linear function that is
linearly independent of the first selected linear function. This can be from a
refinement level immediately below the highest level, as is the case for
$T_{3}$ and $T_{4}$, or from an even lower level in the element hierarchy, as
is the case for $T_{2}$. Figures 3b and Figure 3c show the block $B_{3}$ on
element $T_{3}$ before and after truncation, respectively. Note that when the
domain in Figure 3a is cut either from the left or the right on any of the
elements, the truncated Additive-Schwarz block described here always spans all
possible functions that are only supported on the cut element and form the
problematic eigenmodes for the conditioning of FCM. The truncation in multiple
dimensions is achieved by the tensor product of the procedure in one
dimension, yielding blocks of $(p+1)^{d}$ basis functions. Moreover, this
choice of blocks not only reduces the size of each block, but also ensures
that the number of blocks a single basis function belongs to is at most
$2^{d}$. For clarity it is indicated in all numerical examples in Section 5
whether the full or truncated blocks are applied.

###### Remark 2

Multi-level $hp$-bases are well-conditioned for boundary fitted meshes, due to
the orthogonality properties of the superposed linear basis functions and the
integrated Legendre basis functions [10]. In an immersed setting these
properties diminish, and when no fictitious stiffness is applied (i.e.
$\alpha=0$) even full linear dependencies of the linear basis functions can
occur in specific cut scenarios. This generates a nullspace resulting in a
singular matrix $\mathbf{A}$, and therefore restricts direct solvers to
strictly $\alpha>0$. Iterative solution methods based on Krylov subspaces
neglect these nullspaces. The PCG solver that is applied here minimizes the
energy in the range $R(\boldsymbol{A})$, e.g. [57, 14]. Even when
$\boldsymbol{A}$ is singular, the right hand side $\boldsymbol{b}$ lies in the
range $R(\boldsymbol{A})$, see Eq. (1). Therefore (preconditioned) iterative
solvers are robust to setting $\alpha=0$.

#### Optimizing the number of blocks

The bulk of the computational effort in the construction of the preconditioner
consists of inverting the sub-matrices in Eq. (6) and storing the inverses,
and is proportional to the number of blocks in $\mathcal{B}$. When a block is
devised for _every_ element in the mesh, this leads to a large computational
cost. Repairing the specific ill-conditioning effects of FCM does not require
a block for _every_ element, as can be explained by the Additive-Schwarz lemma
in Section 3.1 and is demonstrated in [41] by only devising blocks for cut
elements. In this work we explore whether the computational cost can be
further reduced by only devising blocks for cut elements $T$ whose volume
fraction $\eta_{T}$

$\eta_{T}=\frac{|T\cap\Omega_{\text{phys}}|}{|T|}\,,$ (8)

is smaller than a threshold $\bar{\eta}$. Note that this generally results in
basis functions in the physical domain that are not contained in any of the
blocks, for which Jacobi preconditioning is applied. In this case, the Jacobi
preconditioning can technically be interpreted as devising a separate block
for every such function. Choosing a smaller value of $\bar{\eta}$ reduces the
number of blocks, which yields a more sparse preconditioner and thereby
reduces the computational cost. However, smaller values of $\bar{\eta}$ allow
smaller untreated elements and with that large condition numbers in accordance
with Eq. (4). This reduces the effectiveness of the preconditioner and
increases the required number of iterations in an iterative solver. The effect
of $\bar{\eta}$ on the memory usage, number of iterations and computation time
is investigated in Section 5.2 .

#### Stabilization using pseudo inverses

The construction of the preconditioner in Eq. (6) requires the inversion of
sub-matrices. All these sub-matrices should theoretically be symmetric
positive definite, because of the coercive and symmetric weak form in Eq. (1).
Since cut elements can be arbitrarily small, these sub-matrices can contain
arbitrarily small eigenvalues, however. When the difference between the
largest and smallest eigenmode in a sub-matrix is of the order of the machine
precision or smaller, a small eigenvalue that should theoretically be positive
can become negative because of rounding errors due to finite precision.
Inverting a sub-matrix that contains a negative eigenvalue with an arbitrarily
small magnitude yields an inverse with an arbitrary large negative eigenmode.
As will be studied numerically in more detail in Section 5.2 , Conjugate
Gradient solvers may diverge when such a mode becomes a dominant part of the
error. This problem can be remedied by replacing the inverse of
$\boldsymbol{A}_{B}$ with a _pseudo_ inverse. Because of the symmetric
positive definite nature of $\boldsymbol{A}_{B}$ we can write

$\boldsymbol{A}_{B}=\sum_{k=1}^{m_{B}}\lambda_{k}\boldsymbol{v}_{k}\otimes\boldsymbol{v}_{k}\hskip
8.53581pt,\hskip
8.53581pt\boldsymbol{A}_{B}^{-1}=\sum_{k=1}^{m_{B}}\frac{1}{\lambda_{k}}\boldsymbol{v}_{k}\otimes\boldsymbol{v}_{k}\,,$
(9)

with $m_{B}$ the dimension of $\boldsymbol{A}_{B}$ and $\lambda_{k}$ and
$\boldsymbol{v}_{k}$ the $k$-th eigenvalue and eigenvector of
$\boldsymbol{A}_{B}$, respectively. The pseudo inverse that replaces
$\boldsymbol{A}_{B}^{-1}$ and stabilizes the preconditioning technique is
defined as

$\boldsymbol{A}_{B+}^{-1}=\sum_{k=1}^{m_{B}}\lambda_{k+}^{-1}\boldsymbol{v}_{k}\otimes\boldsymbol{v}_{k}\,,$
(10)

with

$\lambda_{k+}^{-1}=\begin{cases}\lambda_{k}^{-1}&\text{for
}\lambda_{k}>\lambda_{\text{tres}}\,,\\\ 0&\text{for
}\lambda_{k}\leq\lambda_{\text{tres}}\,,\end{cases}$ (11)

and $\lambda_{\text{tres}}=\varepsilon\max_{k}(\lambda_{k})$ with
$\varepsilon$ here selected as $=10^{-13}$ such that the smallest eigenvalue
that is inverted is sufficiently larger than machine precision. Although the
proposed algorithm requires an eigenvalue decomposition, see Eq. (10), which
is more computationally intensive than the computation of a standard inverse,
this routine is not a bottleneck and the setup can be easily accelerated
through parallelization as elaborated in Section 4 . In Section 5.2 the
stabilized and unstabilized preconditioner are compared, and it is
demonstrated that the stabilized preconditioner does not suffer from the
convergence problems of the Conjugate Gradient solver that the unstabilized
preconditioner suffers from. All other computations in Section 5.2 only apply
the stabilized preconditioner.

###### Remark 3

It should be noted that the stabilized preconditioner with pseudo inverses
generally contains a nullspace. Therefore iterative solvers based on Krylov
subspaces such as the Conjugate Gradient method theoretically do not fully
converge to the exact solution, unless $\boldsymbol{b}$ is exactly in the
range $R(\mathbf{A})$, e.g. [57, 14]. However, even if these modes are a part
of the solution, this effect is very small in practice because the eigenvalues
of the modes disregarded are at least a factor $\varepsilon$ smaller than the
largest eigenvalues in the system. To quantify this effect: _(i)_ the relative
residual can still converge to at least $\varepsilon$ (which is set to be
smaller than the tolerance of the solver), _(ii)_ the relative energy error as
defined in Section 5 still converges to at least $\sqrt{\varepsilon}$ and
_(iii)_ the relative preconditioned residual can theoretically still fully
converge to $0$.

##  4 Implementation and parallelization

A major advantage of the proposed preconditioning technique is its suitability
for parallel computing, as operations for each block $B$ can be performed
independently with minimal synchronization. This can be done using both shared
and distributed memory parallelism, as described in the following sections.

###  4.1 Shared memory parallelism

Different shared memory paradigms can be used to accelerate the computation of
$\boldsymbol{S}$ on shared memory systems. We opt for the use of OpenMP
pragmas [58] to speed up the construction by parallelizing the loop over the
basis function blocks shown in lines 3 to 7 of Algorithm 1 . A short study
showing the OpenMP scalability of the computation of $\boldsymbol{S}$ will be
discussed in Section 5.2 .

###  4.2 Distributed memory parallelism

Distributed memory systems are essential for overcoming the time and memory
bottlenecks related to the computation of large engineering problems [59]. In
such computations, the computational mesh $\mathcal{T}$ is divided over a set
of processes that communicate with one another at well-defined synchronization
points via message passing libraries such as the Message Passage Interface
(MPI) [60]. The locality of the developed preconditioner allows it to be
easily used in a distributed memory setting. Each MPI-process is identified by
its _rank_ and assigned a portion $\mathcal{T}_{\textrm{rank}}$ of the
computational mesh, with $\mathcal{T}_{\textrm{rank}}\subset\mathcal{T}$. We
adopt a hybrid approach combining MPI and OpenMP, where each MPI-process can
use OpenMP threads to speed up local computations on every
$\mathcal{T}_{\textrm{rank}}$.

Although iterative solvers are easier to parallelize than direct solvers, e.g.
[61, 14], they face the challenge of robustness w.r.t. the number of processes
used [59]. A typical caveat of parallel preconditioned iterative solvers is an
increase in the number of iterations needed for convergence when the number of
processes utilized in a simulation increases. Because of the sensitivity of
Additive-Schwarz preconditioning for the finite cell method caused by the
inversion of small cut element contributions, we apply two layers of ghost
elements (elements present on a process but owned by another process) to
ensure that the exact theoretical preconditioner is obtained regardless of the
number of processes or the topology of the partitions
$\mathcal{T}_{\textrm{rank}}$. This has the additional advantage that the
number of iterations needed for convergence is independent of the number of
processes used and thus favors good parallel scalability. It should be noted
that the evaluation of the ghost element layers does not significantly
increase the setup time of the linear system, since the number of ghost
elements on a process is much smaller than the number of active elements.
Moreover, since the bulk of the execution time is spent in solving the linear
system, the iterative solver’s scalability largely governs the overall
scalability of the simulation.

Figure 4 depicts the partitioning of a computational mesh between two MPI-
processes. The red and blue elements represent active elements on
$\mathcal{T}_{0}$ and $\mathcal{T}_{1}$ owned by process 0 and process 1
respectively. The hatched elements constitute the first layer of ghost
elements. This layer allows a process to compute the full contributions in
$\mathbf{A}$ associated with DOFs it owns without communicating with other
processes. This is illustrated in Figure 4b where process 0 can independently
compute all entries in $\mathbf{A}$ associated with DOF 11 by integrating all
the elements around this node. Likewise, the second layer of ghost elements
ensures that the exact theoretical preconditioner $\boldsymbol{S}$ can be
constructed without communication. Figure 4b also illustrates the importance
of this second ghost-element layer. Since DOF 11 shares a support with DOFs 9,
10 and 12, process 0 requires the complete entries in $\mathbf{A}$ associated
with these DOFs to ensure that all entries in $\mathbf{S}$ related to DOF 11
are computed independent of the partitioning as per Eq. (6). This guarantees
that the convergence of an iterative solution method is independent of the
number of processes and how the elements are divided over these. The above
explanation considers uniform meshes, but can trivially be extended to locally
refined grids. The parallel performance of a preconditioned parallel Conjugate
Gradient solver using the described approach will be investigated in Section
5.3.3 .

Figure 4: Partitioning of the computational mesh in parallel computations for
the setup of $\boldsymbol{A}$ and $\boldsymbol{S}$.

###  4.3 Algorithm

The preconditioner construction is summarized in Algorithm 1 . The volume
fraction $\eta_{T}$ of an element is used as a criterion to decide whether a
block $B$ is devised for the basis functions supported on it. Note that the
algorithm does not explicitly devise a separate block of one basis function
for the basis functions that are not present in any of the truncated blocks
$B\in\mathcal{B}$ — as the diagonal scaling of these basis functions can be
interpreted theoretically, Section 3 — but that a separate diagonal scaling
step is done for these basis functions in lines $8$ to $14$.

Algorithm 1 Computation the preconditioner $\boldsymbol{S}$

1:function computePreconditioner ( $\boldsymbol{A}$, $\\{\eta_{T}\\}$,
$\bar{\eta}$ )

2: $\mathcal{B}$ = getBlocks( $\boldsymbol{A}$, $\\{\eta_{T}\\}$, $\bar{\eta}$
)

3: for $B\in\mathcal{B}$ do $\triangleright$ ( Loop over basis function blocks
)

4: $\boldsymbol{A}_{B}$ = extractSubMatrix( $\mathbf{A}$, $B$ )

5: $\boldsymbol{A}_{B+}^{-1}$ = computePseudoInverse( $\boldsymbol{A}_{B}$ )

6: $\boldsymbol{S}$ += prolongateMatrix(
$\boldsymbol{S},\boldsymbol{A}_{B+}^{-1}$)

7: end

8: $n$ = size$(\boldsymbol{A})$

9: $l=1$

10: for $l<n+1$ do $\triangleright$ ( Diagonal scaling )

11: if ( $\boldsymbol{S}(l,l)==0$ ) do

12: $\boldsymbol{S}(l,l)$ = $\dfrac{1.0}{\boldsymbol{A}(l,l)}$

13: $l$ += $1$

14: end

15: return $\boldsymbol{S}$

##  5 Numerical examples

This section demonstrates the suitability of the proposed preconditioner for
finite cell analysis of real-life problems. The numerical examples are
specifically chosen to highlight various aspects presented in the previous
sections. Focus is placed on image-based geometries (in Sections 5.2 and 5.3),
an application field where the finite cell method is particularly advantageous
in circumventing laborious mesh generation procedures. The first example
addresses the extension of the Additive-Schwarz preconditioning technique to
FCM problems involving multi-level $hp$-refinement by studying a simple
example of a cube with a spherical cavity subjected to uniaxial loading. The
second example, from the field of biomechanics, considers the loading of a
lumbar vertebral body and studies the effect of the pseudo inverse and the
volume fraction threshold $\bar{\eta}$ on the computational cost of a
preconditioned Conjugate Gradient solver (PCG). The final example brings
together FCM, multi-level $hp$-refinement, Additive-Schwarz preconditioning
and parallel computing in the image-based analysis of a die cast gearbox
housing.

The convergence of the PCG solver applied in this section is assessed by
monitoring the development of the relative preconditioned residual

$\frac{\|\mathbf{S}\boldsymbol{r}_{i}\|}{\|\mathbf{S}\boldsymbol{b}\|}=\frac{\|\mathbf{S}\boldsymbol{b}-\mathbf{S}\boldsymbol{A}\boldsymbol{x}_{i}\|}{\|\mathbf{S}\boldsymbol{b}\|}\,,$
(12)

and when possible the relative error in the energy norm

$e_{i}=\frac{\|\boldsymbol{x}_{\text{ref}}-\boldsymbol{x}_{i}\|_{\boldsymbol{A}}}{\|\boldsymbol{x}_{\text{ref}}\|_{\boldsymbol{A}}}=\frac{\sqrt{\left(\boldsymbol{x}_{\text{ref}}-\boldsymbol{x}_{i}\right)\boldsymbol{A}\left(\boldsymbol{x}_{\text{ref}}-\boldsymbol{x}_{i}\right)^{T}}}{\sqrt{\boldsymbol{x}_{\text{ref}}\boldsymbol{A}\boldsymbol{x}_{\text{ref}}^{T}}}\,,$
(13)

with a reference solution $\boldsymbol{x}_{\text{ref}}$ obtained from the
parallel direct solver Intel® Pardiso, contained in the Intel Math Kernel
Library [62]. The subscript $i$ in Equations (12) and (13) denotes quantities
in the $i$-th Conjugate Gradient iteration. In [41] the effectiveness of the
preconditioning technique is assessed by focusing on the condition number of
the system. We do not adopt this approach here but rather focus of the
residual and energy convergence since the cost of computing condition numbers
for the large systems under investigation is prohibitive.

###  5.1 Compression of a cube with a spherical exclusion

The following example illustrates the importance of selecting suitable blocks
for Additive-Schwarz preconditioning in FCM simulations with multi-level
$hp$-refinement as described in Section 3.2 . To this end, a simple example
comprised of a cube with a spherical cavity under compressive loading is
considered. Although this setup yields a relatively small system, with a small
number of DOFs ($<17\,000$), it serves as a good starting point to illustrate
the effectiveness of the preconditioning technique developed in this
manuscript as it already shows the conditioning problems associated with FCM
and multi-level $hp$-refinement.

Figure 5: Cube with a spherical cavity under compressive loading.

#### Problem setup

A cube of unit length with a spherical cavity of radius $r=0.01$ is subjected
to a homogeneous pressure load $P$ as shown in Figure 5 . The cube has a
Young’s modulus of $70$ GPa and Poisson’s ratio $\nu=0.34$. Homogeneous
Dirichlet boundary conditions in the normal direction are applied using a
penalty parameter $\beta=10^{10}$. The base mesh consists of $8\times 8\times
8$ elements and is refined toward the surface of the exclusion with a
refinement depth $k\in\\{0,1,2,3,4\\}$. A polynomial order $p=3$ is used,
resulting in the numbers of unknowns for the different refinement depths given
in Table 1 . Numerical integration is performed using octree partitioning with
a tree-depth of 3 relative to the local element size, see e.g. [63] for
details.

$k$ | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---
DOFs | 13 851 | 14 541 | 15 231 | 15 921 | 16 611
Table 1: Number of degrees of freedom for different refinement depths $k$.

#### Convergence behavior

The convergence of the PCG solver preconditioned with the full and the
truncated blocks is shown in Figure 6 . The results show that the convergence
of the residual closely follows the convergence of the energy. This is
generally not the case for ill-conditioned systems such as unpreconditioned
FCM, e.g. [15], and indicates that both preconditioners are robust w.r.t.
small cut cells. Figures 6a and 6c , however, show the effect of local
refinements on the convergence speed of the solver when no truncation is used
in the construction of $\mathbf{S}$. This indicates that the preconditioner
with full blocks is not robust w.r.t. the refinement depth — i.e. a small
increase in the number of unknowns, by approximately a factor of 1.2, leads to
a significant increase in the number of iterations needed to reach
convergence, which increases by a factor of around 3.6. The intensity of this
effect is significantly reduced when the preconditioner with truncated blocks
is used as demonstrated in Figures 6b and 6d . The truncated preconditioner is
therefore not only computationally less expensive than the full preconditioner
(due to the smaller block sizes), but also yields better conditioning. Similar
results are obtained using a more complex geometry in the third example
(Section 5.3).

(a) Full preconditioner: residual convergence.

(b)  Truncated preconditioner: residual convergence.

(c) Full preconditioner: energy convergence.

(d) Truncated preconditioner: energy convergence.

Figure 6:  Convergence behavior of the PCG solver with full and truncated
blocks for multi-level $hp$-refinement. Note that the convergence of the
residual is not monotone as the PCG solver minimizes the energy and the
residual is only loosely bound to this, e.g. [14] .

###  5.2 Image-based simulation of a lumbar vertebra

The applicability of the preconditioner in biomechanical simulations is
highlighted in the following image-based simulation of a lumbar vertebra
subjected to compressive loading, with the computational setup following [13].
Elastostatic computations are performed with the described setup, utilizing a
PCG solver for the solution of the linear system of equations with either no
preconditioner, diagonal scaling, or the presented Additive-Schwarz
preconditioner. Various aspects of the iterative solver’s convergence behavior
are monitored during the simulation for different solver and preconditioner
configurations. This example does not consider local mesh refinement, which
yields a relatively simple problem on a complex geometry, which is a good test
case to benchmark the stability of the preconditioning technique and the
efficiency of different threshold volume fractions $\bar{\eta}$, Section 3.2 .
Note that on the uniform mesh without local refinements truncation of the
blocks is not applicable.

#### Problem setup

The geometry of the lumbar vertebra is obtained via a high-resolution micro-CT
scan of the specimen with a voxel size of $80\times 80\times 80\ \mu m^{3}$. A
subsequent segmentation of the scan using ITK-SNAP [64] yields a model of the
vertebral body without the surrounding soft tissue, Figure 7a . It is noted
that the segmentation only considers the outer boundary of the vertebra and no
distinction is made between cortical and trabecular bone.

A non-boundary conforming discretization is generated using the finite cell
method with $10^{3}$ voxels contained in each element resulting in 170 982
elements. Elements completely outside the physical domain are removed in a
preprocessing step yielding 75 821 active elements. A polynomial order of
$p=3$ is chosen, resulting in approximately 1.7 million degrees of freedom.
Material properties within each element are defined on a voxel-level by
applying a threshold on a voxel’s CT-number. The CT-number is a measure of the
x-ray absorption coefficient that can be related to the density and
consequently the stiffness of the material in a voxel. Material identified as
inside the bone structure is assigned a Young’s modulus $E_{\text{bone}}=12$
GPa and Poisson’s ratio $\nu=0.3$, corresponding to values commonly used in
literature [65]. Fictitious material, on the other hand, is assigned the
material properties $E_{\text{fict}}=10^{-11}$ GPa and $\nu=0.3$ in order to
avoid peak stresses in cut elements [66]. The surface of the entire vertebral
body is recovered using a marching cubes algorithm [67] and thereafter trimmed
to provide a triangulation for the superior end-plate, where a uniformly
distributed vertical load of 800N is applied, and the inferior end-plate,
which is weakly clamped using the penalty method with $\beta=10^{6}$.

Material | vertebral bone
---|---
Young’s modulus, $E_{\text{bone}}$ | 12 GPa
Poisson’s ratio, $\nu$ | 0.3
$l\times w\times h$ [mm] | $47.3\times 55.3\times 34.1$

(a) Geometric and material properties of the lumbar vertebra.

(b) Computational mesh of active elements.

(c) Cut showing the porous bone structure and distribution of the CT-number.

Figure 7: Setup of the lumbar vertebra example.

#### Convergence behavior and stability of $\mathbf{S}$

Figure 8 shows the values of the relative preconditioned residual and the
energy error plotted against the number of PCG iterations performed. The
solver is terminated when either a cut-off tolerance in the relative
preconditioned residual of $10^{-10}$ or the maximum number of iterations (300
000) is reached. One can observe that the unpreconditioned system, denoted by
the curves labeled $\mathbf{A}$, suffers from ill-conditioning characterized
by the slow residual convergence rates and high energy norm errors even after
a significant number of iterations. The diagonally preconditioned system,
$\mathbf{D}\mathbf{A}$, shows an improved convergence behavior but also fails
to reach the desired cut-off tolerance within the specified maximum number of
iterations.

(a) Convergence of the relative residual

(b) Convergence of the error in the energy norm

Figure 8:  Convergence of the relative residual and the energy norm errors for
different preconditioners in the lumbar vertebra example.

The curves $\boldsymbol{S}^{*}\boldsymbol{A}$ and
$\boldsymbol{S}\boldsymbol{A}$ depict the systems preconditioned with the
preconditioning scheme developed here using the full inverse
($\boldsymbol{S}^{*}$) and stabilized pseudo inverse ($\boldsymbol{S}$)
respectively, with cut-off tolerance $\varepsilon=10^{-13}$. The importance of
filtering out eigenvalues smaller than the machine precision in the local
inverses $\mathbf{A}_{B+}^{-1}$, as described in Section 3.2 , is clearly
visible from these curves in Figures 8a and 8b . Both
$\boldsymbol{S}^{*}\boldsymbol{A}$ and $\boldsymbol{S}\boldsymbol{A}$ exhibit
the same convergence until the eigenvectors of
$\boldsymbol{S}^{*}\boldsymbol{A}$ with negative eigenvalues become the
dominant part of the error. It is clearly visible that the relative energy
norm error stagnates and the relative residual diverges. These problematic
eigenvalues are discarded by $\boldsymbol{S}\boldsymbol{A}$, resulting in
convergence up to the desired tolerance as indicated in the figures. It should
be noted that potential loss of accuracy due to the discarding of extremely
small eigenmodes is not observed in this example.

It is customary in FCM to preprocess the system and discard all functions that
completely lie outside of the physical domain. By stabilizing with the pseudo
inverse and setting the fictitious stiffness $\alpha=0$, this preprocessing
step is theoretically not required since the pseudo inverse will ignore these
functions. Discarding these functions in a preprocessing step still reduces
the computational cost however. It is noted that with $\alpha$ large enough
the system does not contain eigenvalues of the order of the machine precision,
and the pseudo inverse is equivalent to the full inverse.

#### Optimization: Preconditioning only severely cut cells

In [41] it is mentioned that the number of blocks used in the setup of the
preconditioner can be reduced by considering only elements below a certain
volume fraction. This idea is followed in Section 3.2 where it is proposed to
only invert blocks of basis functions for elements with a volume fraction
smaller than a threshold value $\bar{\eta}\in[0.0,1.0]$. Basis functions that
are not present in any of these elements are instead diagonally scaled, to
reduce the computational cost.

In this section we study the effect of the threshold value $\bar{\eta}$ on the
performance of the preconditioner for the lumbar vertebra model described
above. To this end, the linear system is solved by a PCG solver for different
values of $\bar{\eta}$. A value $\bar{\eta}=1.0$ corresponds to the
preconditioner with a block for every element while $\bar{\eta}=0.0$ results
in a Jacobi preconditioner. It is again stipulated that in this computation
the mesh does not contain local refinements, such that the blocks are not
truncated as described in Section 3.2.

Figure 9 depicts the influence of $\bar{\eta}$ on the convergence of the PCG
solver, the solution time and memory requirements of setting up the
preconditioner. Figures 9a \- 9b confirm that the effectiveness of the
preconditioner, as described in Section 3.2 , is maintained as long as
problematic basis functions are taken into account when constructing
$\mathbf{S}$. This holds for $\bar{\eta}\in[0.4,1.0]$, while values of
$\bar{\eta}\leq 0.4$ allow smaller untreated elements and thus reduce the
effectiveness of $\mathbf{S}$. Reducing the value of $\bar{\eta}$ leads to a
sparser preconditioner which requires less storage space as shown in Figure 9c
, where the size of the preconditioner scales linearly with the number of
blocks in $\mathcal{B}$. The increased sparsity of $\mathbf{S}$ for values of
$\bar{\eta}<1.0$ leads to a reduction in the computational cost of a single
PCG iteration, since the matrix-vector multiplication involving $\mathbf{S}$
and $\boldsymbol{r}$ takes less time. This leads to an overall faster
execution time as shown in Figure 9d , where the number of iterations needed
until convergence is plotted alongside the corresponding solution time for
different values of $\bar{\eta}$. The solver converges in less time for
$0.6\leq\bar{\eta}<1.0$, with the solver converging 1.6 times faster for
$\bar{\eta}=0.7$ than for $\bar{\eta}=1.0$. This study shows that values of
$\bar{\eta}\in[0.6,0.7]$ yield an effective preconditioner which is sparser
and consequently faster than the _complete_ preconditioner for
$\bar{\eta}=1.0$ with a block for _every_ element.

To enable comparison of all presented examples, in the remainder of this
manuscripts the non-optimized choice of $\bar{\eta}=1.0$ will be considered.

(a) Convergence of the relative residual for various $\bar{\eta}$

(b) Convergence of the energy error for various $\bar{\eta}$

(c) Preconditioner size for different values of $\bar{\eta}$

$\cdot 10^{3}$$\cdot 10^{2}$

(d) Convergence behavior for different values of $\bar{\eta}$

Figure 9:  Preconditioning of cut cells with $\eta_{T}\leq\bar{\eta}$.

#### Computational cost of setting up $\mathbf{S}$

Section 4.1 explains how the construction of $\mathbf{S}$ can be accelerated
by the use of shared memory parallelism. This is illustrated in a short study
that shows the OpenMP speed up achieved for a variable number of threads
involved in the setup of the preconditioner. Three different preconditioners
with different numbers of blocks in $\mathcal{B}$, and consequently different
computational costs, are analyzed. The results of this study are shown in
Figure 10 . Although good parallel scalability is achieved for the different
number of blocks, the parallel efficiency decreases gradually with an increase
of the number of threads. This is attributed to the synchronization needed
when assembling the entries of $\mathbf{A}_{B+}^{-1}$ into $\mathbf{S}$, Eq.
(6). This effect can be resolved by the use of graph coloring algorithms [68].

Figure 10: Setup of $\mathbf{S}$ — OpenMP scalability for different numbers of
blocks in $\mathcal{B}$.

###  5.3 Loading of a die cast gearbox housing

Finally we consider an industrial example to portray the suitability of the
preconditioner for the solution of large-scale embedded problems involving
multi-level $hp$-refinement. The example focuses on the loading of an aluminum
gearbox housing produced by die casting, a manufacturing process by which
molten metal is forced under high pressure into a mold cavity producing a net-
shape structure upon solidification of the molten metal. Die casting is widely
used for the production of small- to medium-sized non-ferrous metal parts due
to its high production rates, the quality and dimensional consistency of the
produced parts and its design freedom, that allows the production of parts
with high geometric complexity [69]. The extensive use of die cast parts is,
however, limited by defects such as shrinkage cavities and entrapped gas
bubbles that lead to increased porosity in the parts. Pores negatively affect
the mechanical properties of die cast parts and are the preferred sites for
fatigue-crack initiation.

Several studies have been conducted on the influence of pore size, position
and orientation on the stress state, crack growth and fatigue behavior of
aluminum die cast parts [70, 71, 72]. Finite element analysis of die cast
parts can help quantify the stress concentration in the vicinity of pores
through the use of image-based computations [73]. The limitations of boundary
conforming methods in capturing the complex pore morphology can be overcome by
the use of the finite cell method as investigated in [74, 75, 76].

In this section, we revisit the die cast example considered in [74]. The first
major modification in contrast to the original example is the use of a
preconditioned iterative solver for solving the linear system of equations
instead of a direct solver. This modification allows, due to the smaller
memory requirement, for computations with a higher resolution of the die cast
part. Moreover, distributed memory parallelism is exploited to accelerate the
performance of the preconditioned iterative solver. The second major
modification is the use of multi-level $hp$-refinement to achieve a higher
spatial resolution around the pores in the specimen leading to a better
resolution of the complex stress state in these regions.

#### Problem setup

The geometry of the gearbox housing is obtained from a CT scan with a voxel
size of $198.7\times 198.7\times 198.7\ \mu m^{3}$. A non-boundary conforming
mesh is generated with $10^{3}$ voxels grouped to form a single finite cell.
Elements completely inside the fictitious domain are eliminated yielding a
total of 34 230 active elements in the base mesh before refinement. The
classification of the material properties is done by applying a threshold on
the intensity values (CT-number) of the CT scan. The threshold is chosen such
that the metal material and cavities or internal pores (which are part of the
fictitious domain) can be distinguished, as illustrated in Figure 11b and
Figure 11c . Metal material is assigned a Young’s modulus
$E_{\text{metal}}=70$ GPa and Poisson’s ratio $\nu=0.3$ while fictitious
material is assigned the values $E_{\text{fict}}=10^{-11}$ GPa and $\nu=0.3$.
Figure 12a depicts the computational mesh of the gearbox housing. The
discretization is refined towards the surfaces of the interior pores with a
refinement depth $k\in\\{0,1,2,3\\}$ and a polynomial order $p=4$ is chosen,
resulting in the number of DOFs listed in Table 2 . The numerical integration
in this example is is performed by Gauss quadrature over integration sub-cells
instead of standard octrees, see e.g. [77, 13] for details. Homogeneous
Dirichlet boundary conditions are applied on the dark green surfaces while a
displacement $\bar{u}_{x}=0.1$cm is applied on the red cylindrical surfaces as
depicted in Figure 12b . All boundary conditions are enforced weakly using the
penalty method with the penalty parameter $\beta=10^{6}$.

Material | aluminum alloy
---|---
EN AC-46000
Young’s modulus, $E$ | 70 GPa
Poisson’s ratio, $\nu$ | 0.29
$l\times w\times h$ [mm] | $175\times 101\times 80$

(a) Geometric and material properties of the gearbox housing.

(b) Cut-plane C of the housing CT scan

depicting pores in the specimen.

(c) Distribution of the CT-number

along the cut-line $A$-$B$.

Figure 11: Setup of the die cast gearbox housing example.

(a) Discretization of the gearbox housing.

(b) Application of boundary conditions.

Figure 12: Mesh and boundary conditions of the gearbox housing example. All displacements on the dark green surfaces are fixed while a displacement field $\bar{u}_{x}=0.1$cm is prescribed on the red surfaces. $k$ | 0 | 1 | 2 | 3
---|---|---|---|---
DOFs | 1 662 666 | 1 670 478 | 1 683 792 | 1 716 102
Table 2: Number of degrees of freedom for different refinement depths $k$.

####  5.3.1 Influence of pores on the stress distribution

Before studying the convergence behavior of the iterative solver for different
refinement depths $k$, we highlight the necessity of local refinements in
resolving stress states in the vicinity of small geometric features with an
economical number of DOFs. Different values of $k$ yield comparable results of
the overall displacement and stress fields in the gearbox housing with an
exemplary result given in Figure 13 . Differences arise, however, in the
stress fields in the vicinity of the pores as shown in Figures 14 \- 15 .
Without local mesh refinement, the chosen mesh is unable to accurately capture
the expected stress concentration around the pore boundary, as shown in Figure
15a . The use of a globally $h$-refined mesh would help capture stress
concentration around the pores, but is not a viable option since the number of
DOFs would increase drastically. Moreover, performing $h$-refinement would not
significantly affect the stress distribution away from the pores, since this
state is already adequately represented by the coarse mesh as shown in Figure
14 . The stated arguments motivate the use of $hp$-refinement to capture local
solution stress distributions in a computationally efficient manner. The
results of the Von Mises stress available in Figure 15b and Figure 14 show
that the location of the peak stress shifts to the boundary of the pores upon
application of multi-level $hp$-refinement, which is physically realistic
behavior.

(a) Displacement magnitude.

(b) Von Mises stress distribution.

Figure 13: Results of the displacement and von Mises stress fields of the
gearbox housing. Figure 14: Comparison of the von Mises stress in the vicinity
of the pores along the cut-line $A$-$B$ shown in Figure 11b .

(a) Distribution of the von Mises stress in cut-plane $C$ for $k=0$.

(b) Distribution of the von Mises stress in cut-plane $C$ for $k=4$.

Figure 15: Influence of the refinement depth on the stress distribution around
pores in cut-plane $C$.

####  5.3.2 Convergence behavior

An analogous study to that in Section 5.1 is carried out to analyze the
influence of multi-level $hp$-refinement on the convergence of the PCG solver.
We refine the computational mesh towards the internal pores with a refinement
depth $k\in\\{0,1,2,3\\}$. The development of the relative preconditioned
residual and energy error is monitored in computations involving the full and
truncated preconditioners introduced in Section 3.2 . Figure 16 portrays the
convergence behavior of the full and truncated preconditioner for different
levels of multi-level $hp$-refinement of the gearbox housing mesh. The results
show similar convergence behavior as in the numerical example considered in
Section 5.1 . When the full preconditioner is used, the number of iterations
drastically increases with the refinement depth. The truncated preconditioner
performs better than the full preconditioner since it only takes into account
basis functions that can potentially become almost linearly dependent as
described in Section 3.2 . Table 3 demonstrates that $n_{\textrm{max}}$, the
maximum number of blocks a basis function belongs to, is bounded and equals
$2^{d}=8$ for the truncated preconditioner.

(a) Full preconditioner: residual convergence.  (b) Truncated preconditioner:
residual convergence.  (c) Full preconditioner: energy convergence.  (d)
Truncated preconditioner: energy convergence.

Figure 16:  Comparison of the convergence behavior of the full and truncated preconditioners for multi-level $hp$-refinement. $k$ | 0 | 1 | 2 | 3
---|---|---|---|---
Full blocks | 8 | 64 | 274 | 810
Truncated blocks | 8 | 8 | 8 | 8
Table 3: Maximum overlap $n_{\textrm{max}}$ of the full and truncated
Additive-Schwarz blocks for different refinement depths $k$ in the gearbox
housing example.

####  5.3.3 Analysis of the parallel scalability

Section 4.2 describes the use of the preconditioner in a hybrid (distributed
and shared memory) parallel setting. The study at hand aims to show the
performance of a parallel PCG solver in large-scale finite cell computations.
We use the parallel PCG solvers available in Trilinos [78] in conjunction with
the presented preconditioner. Two different discretizations of the gearbox
housing are considered. The first discretization is the mesh described in
Section 5.3 consisting of 34 230 elements and approximately 1.6 million
degrees of freedom. The second mesh is a finer discretization of the gearbox
that consists of 83 757 elements formed by grouping $7^{3}$ voxels to form a
single finite cell. A polynomial order of $p=4$ is chosen here as well,
resulting in approximately 3.98 million DOFs. Boundary conditions are applied
for both discretizations as illustrated in Figure 12b with a penalty value
$\beta=10^{6}$ . The numerical computations are performed on the CoolMAC
cluster at the Technical University of Munich equipped with four AMD Bulldozer
Opteron 6274 CPUs and 256 GB memory per node. These hybrid computations are
carried out by varying the number of nodes from 1 to 16 for the coarser
discretization and from 2 to 16 for the fine discretization. Each node is made
up of 64 processing units (cores) i.e. one MPI-process $\times$ 64 OpenMP
threads.

Figure 17 shows the execution time of the PCG solver plotted against the
number of processing units used in the hybrid computations for both
discretizations of the gearbox housing. The preconditioned solver shows
excellent parallel scalability for both discretizations up to 256 cores where
superlinear speed up is achieved due to cache effects when moving to more
cores. This trend can be maintained as long as the computational work in every
PCG iteration — the matrix vector multiplications — is significantly larger
than the communication overhead. This is shown by the loss of efficiency in
the computations involving 1.7 million degrees of freedom and more than 256
cores. The fine discretization provides adequate computational work and
maintains a parallel efficiency of 115% up to 1024 processing units.

Figure 17: Parallel scalability study of the PCG solver for two different
problem sizes.

##  6 Conclusion

The presented work describes a preconditioning technique that enables
iterative solutions of large finite cell method (FCM) computations with multi-
level $hp$-refined discretizations. An Additive Schwarz preconditioning
technique based on overlapping element-wise blocks is developed and tailored
to the hierarchical structure of the $hp$-discretization scheme. Furthermore,
a robust preconditioner is obtained through the use of a pseudo inverse and
the efficiency of the preconditioner improved by utilizing parallelism. The
results show that the preconditioner is robust to local refinements and to how
elements are cut. It is verified that the computational cost of setting up and
applying the preconditioner scales well with the number of processors in a
parallel setting. Therefore the preconditioning technique increases tolerable
sizes of FCM systems, which extends the scope of problems that can be solved
with FCM to real-life applications.

The current work is restricted to problems which yield symmetric positive
definite matrices. Notwithstanding that this covers a large range of
engineering problems, this does not include multivariate problems with
indefinite matrices or problems involving convection with nonsymmetric
matrices. Applying or modifying the preconditioner for such applications, as
done for simple meshes in [41], is therefore an important topic of future
research. Furthermore, the current setup is tailored for the multi-level
$hp$-refinement scheme introduced in [9]. A relevant topic for further
research is the development and testing of similar techniques for different
refinement schemes.

Finally the preconditioner does not resolve the conditioning effects that can
be expected from increasing the number of elements, e.g. [79]. Be it that the
preconditioner is robust to cut elements and thereby extends the numbers of
DOFs that FCM systems can currently be solved with, we do not expect it to be
efficient for exascale systems. In order to further extend the range of
solvable systems, it is an interesting topic for future research to apply the
preconditioner together with a multi-grid or domain-decomposition approach.

## Acknowledgements

The authors gratefully acknowledge the financial support of the NWO under the
Graduate program Fluid & Solid Mechanics and the European Research Council
under Grant ERC-2014- StG 637164.

## References

  * [1] J. Parvizian, A. Düster, and E. Rank, “Finite cell method,” Computational Mechanics, vol. 41, no. 1, pp. 121–133, 2007.
  * [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, vol. 197, no. 45–48, pp. 3768–3782, 2008.
  * [3] 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, vol. 22, no. 3, pp. 391–455, 2014.
  * [4] 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 & Fluids, vol. 141, pp. 135 – 154, 2016.
  * [5] V. Varduhn, M. 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, vol. 107, no. 12, pp. 1054–1079, 2016.
  * [6] A. Düster, E. Rank, and B. A. Szabó, “The p-version of the finite element method and finite cell methods,” in Encyclopedia of Computational Mechanics, vol. 2, pp. 1–35, Chichester, West Sussex: John Wiley & Sons, 2017.
  * [7] E. Burman and P. Hansbo, “Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method,” Applied Numerical Mathematics, vol. 62, no. 4, pp. 328 – 341, 2012.
  * [8] 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, vol. 104, no. 7, pp. 472–501, 2014.
  * [9] 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, vol. 55, no. 3, pp. 499–517, 2015.
  * [10] N. Zander, T. Bog, M. Elhaddad, F. Frischmann, S. Kollmannsberger, and E. Rank, “The Multi-Level Hp-Method for Three-Dimensional Problems: Dynamically Changing High-Order Mesh Refinement with Arbitrary Hanging Nodes,” Computer Methods in Applied Mechanics and Engineering, vol. 310, pp. 252–277, 2016.
  * [11] N. Zander, M. Ruess, T. Bog, S. Kollmannsberger, and E. Rank, “Multi-Level Hp-Adaptivity for Cohesive Fracture Modeling,” International Journal for Numerical Methods in Engineering, vol. 109, no. 13, pp. 1723–1755, 2016\.
  * [12] T. Bog, N. Zander, S. Kollmannsberger, and E. Rank, “Weak Imposition of Frictionless Contact Constraints on Automatically Recovered High-Order, Embedded Interfaces Using the Finite Cell Method,” Computational Mechanics, vol. 61, pp. 385–407, Apr 2018.
  * [13] 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, vol. 34, no. 4, p. e2951, 2018.
  * [14] Y. Saad, Iterative Methods for Sparse Linear Systems.  Philadelphia, PA, USA: Society for Industrial and Applied Mathematics, 2nd ed., 2003.
  * [15] F. de Prenter, C.V. Verhoosel, G.J. van Zwieten, and E.H. van Brummelen, “Condition number analysis and preconditioning of the finite cell method,” Computer Methods in Applied Mechanics and Engineering, vol. 316, pp. 297 – 327, 2017.
  * [16] D. Schillinger, L. Dedè, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, and T. J. 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, vol. 249-252, pp. 116–150, 2012.
  * [17] 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, vol. 249-252, pp. 104–115, 2012.
  * [18] 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, vol. 95, no. 10, pp. 811–846, 2013.
  * [19] 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, vol. 269, pp. 46–71, 2014.
  * [20] M. Dauge, A. Düster, and E. Rank, “Theoretical and Numerical Investigation of the Finite Cell Method,” Journal of Scientific Computing, vol. 65, no. 3, pp. 1039–1064, 2015.
  * [21] 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, vol. 284, pp. 138–164, 2015.
  * [22] D. Elfverson, M. G. Larson, and K. Larsson, “CutIGA with basis function removal,” Advanced Modeling and Simulation in Engineering Sciences, vol. 5, no. 1, p. 6, 2018.
  * [23] E. Burman, “Ghost penalty,” Comptes Rendus Mathematique, vol. 348, no. 21, pp. 1217 – 1220, 2010.
  * [24] K. Höllig, U. Reif, and J. Wipper, “Weighted Extended B-Spline Approximation of Dirichlet Problems,” SIAM Journal on Numerical Analysis, vol. 39, no. 2, pp. 442–462, 2001.
  * [25] K. Höllig, C. Apprich, and A. Streit, “Introduction to the Web-method and its applications,” Advances in Computational Mathematics, vol. 23, no. 1, pp. 215–237, 2005.
  * [26] R. Sanches, P. Bornemann, and F. Cirak, “Immersed b-spline (i-spline) finite element method for geometrically complex domains,” Computer Methods in Applied Mechanics and Engineering, vol. 200, no. 13, pp. 1432 – 1445, 2011.
  * [27] T. Rüberg and F. Cirak, “Subdivision-stabilised immersed b-spline finite elements for moving boundary flows,” Computer Methods in Applied Mechanics and Engineering, vol. 209-212, pp. 266 – 283, 2012.
  * [28] 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, vol. 74, no. 9, pp. 623–660, 2013.
  * [29] T. Rüberg, F. Cirak, and J. M. García Aznar, “An unstructured immersed finite element method for nonlinear solid mechanics,” Advanced Modeling and Simulation in Engineering Sciences, vol. 3, no. 1, p. 22, 2016.
  * [30] 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, vol. 336, pp. 533–553, 2018.
  * [31] B. Marussig, R. Hiemstra, and T. J. Hughes, “Improved conditioning of isogeometric analysis matrices for trimmed geometries,” Computer Methods in Applied Mechanics and Engineering, vol. 334, pp. 79 – 110, 2018.
  * [32] T. Belytschko and T. Black, “Elastic crack growth in finite elements with minimal remeshing,” International Journal for Numerical Methods in Engineering, vol. 45, no. 5, pp. 601–620, 1999.
  * [33] N. Moës, J. Dolbow, and T. Belytschko, “A finite element method for crack growth without remeshing,” International Journal for Numerical Methods in Engineering, vol. 46, no. 1, pp. 131–150, 1999.
  * [34] 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, vol. 64, no. 8, pp. 1033–1056, 2005.
  * [35] A. Menk and S. Bordas, “A robust preconditioning technique for the extended finite element method,” International Journal for Numerical Methods in Engineering, vol. 85, no. 13, pp. 1609–1632, 2011.
  * [36] B. Hiriyur, R. Tuminaro, H. Waisman, E. Boman, and D. Keyes, “A Quasi-algebraic Multigrid Approach to Fracture Problems Based on Extended Finite Elements,” Siam Journal on Scientific Computing, vol. 34, no. 2, pp. 603–626, 2012.
  * [37] L. Berger-Vergiat, H. Waisman, B. Hiriyur, R. 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, vol. 90, no. 3, pp. 311–328, 2012.
  * [38] H. Waisman and L. Berger-Vergiat, “An adaptive domain decomposition preconditioner for crack propagation problems modeled by XFEM,” International Journal for Multiscale Computational Engineering, vol. 11, no. 6, pp. 633–654, 2013.
  * [39] C. Lehrenfeld and A. Reusken, “Optimal preconditioners for Nitsche-XFEM discretizations of interface problems,” Numerische Mathematik, vol. 135, no. 2, pp. 313–332, 2017.
  * [40] S. Badia and F. Verdugo, “Robust and scalable domain decomposition solvers for unfitted finite element methods,” Journal of Computational and Applied Mathematics, 2017.  In Press.
  * [41] F. de Prenter, C.V. Verhoosel, and E.H. van Brummelen, “Preconditioning immersed isogeometric finite element methods with application to flow problems,” ArXiv e-prints, 2017.
  * [42] S. Brenner and L. Scott, The Mathematical Theory of Finite Element Methods.  Springer, 2008.
  * [43] 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, vol. 200, no. 47-48, pp. 3358–3380, 2011.
  * [44] D. Kamensky, M.-C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, and T. J. Hughes, “An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves,” Computer Methods in Applied Mechanics and Engineering, vol. 284, pp. 1005 – 1053, 2015.
  * [45] J. Nitsche, “Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind,” Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg, vol. 36, no. 1, pp. 9–15, 1971.
  * [46] A. Embar, J. Dolbow, and I. Harari, “Imposing Dirichlet boundary conditions with Nitsche’s method and spline based finite elements,” International Journal for Numerical Methods in Engineering, vol. 83, no. 7, pp. 877–898, 2010.
  * [47] I. Babuška, “The Finite Element Method with Penalty,” Mathematics of Computation, vol. 27, no. 122, p. 221, 1973.
  * [48] P. Šolín, Higher-Order Finite Element Methods.  Studies in advanced mathematics, Boca Raton: Chapman & Hall/CRC, 2004\.
  * [49] W. Rachowicz, D. Pardo, and L. Demkowicz, “Fully automatic hp-adaptivity in three dimensions,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 37–40, pp. 4816–4842, 2006.
  * [50] L. Demkowicz, Computing with Hp-Adaptive Finite Elements, Vol. 1: One and Two Dimensional Elliptic and Maxwell Problems.  Applied mathematics and nonlinear science series, Boca Raton: Chapman & Hall/CRC, 2007.
  * [51] B. Smith, P. Bjørstad, and W. Gropp, Domain Decomposition: parallel multilevel methods for elliptic partial differential equations.  Cambridge university press, 1996.
  * [52] A. Toselli and O. Widlund, Domain Decomposition Methods: Algorithms and Theory.  Springer, 2005.
  * [53] R. Ferencz and T. Hughes, Solutions in Nonlinear Solid Mechanics.  Elsevier, 1998.
  * [54] A. Matsokin and S. Nepomnyaschikh, “The Schwarz alternation method in a subspace,” Soviet Mathematics (Izv. Vyssh. Uchebn. Zaved. Mat.), 1985.
  * [55] P. Lions, “On the Schwarz alternating method. I,” in First international symposium on domain decomposition methods for partial differential equations, pp. 1–42, 1988.
  * [56] B. A. Szabó, A. Düster, and E. Rank, “The P-Version of the Finite Element Method,” in Encyclopedia of Computational Mechanics, Chichester, West Sussex: John Wiley & Sons, 2004.
  * [57] E. F. Kaasschieter, “Preconditioned Conjugate Gradients for Solving Singular Systems,” J. Comput. Appl. Math., vol. 24, no. 1-2, pp. 265–275, 1988\.
  * [58] OpenMP Architecture Review Board, “OpenMP Application Program Interface Version 3.0,” 2008.
  * [59] M. Ferronato, “Preconditioning for Sparse Linear Systems at the Dawn of the 21st Century: History, Current Developments, and Future Perspectives,” International Scholarly Research Notices Applied Mathematics, vol. 2012, p. 49, 2012.
  * [60] M. P. Forum, “MPI: A Message-Passing Interface Standard,” tech. rep., Knoxville, TN, USA, 1994.
  * [61] J. J. Dongarra, L. S. Duff, D. C. Sorensen, and H. A. V. Vorst, Numerical Linear Algebra for High Performance Computers.  Philadelphia, PA, USA: Society for Industrial and Applied Mathematics, 1998.
  * [62] Intel, “Intel Math Kernel Library.” https://software.intel.com/en-us/mkl.
  * [63] 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, vol. 50, no. 4, pp. 445–478, 2012.
  * [64] P. A. Yushkevich, J. Piven, H. Cody Hazlett, R. Gimpel Smith, S. Ho, J. C. Gee, and G. Gerig, “User-Guided 3D Active Contour Segmentation of Anatomical Structures: Significantly Improved Efficiency and Reliability,” Neuroimage, vol. 31, no. 3, pp. 1116–1128, 2006.
  * [65] T. M. Keaveny, E. F. Morgan, G. L. Niebur, and O. C. Yeh, “Biomechanics of Trabecular Bone,” Annual Review of Biomedical Engineering, vol. 3, no. 1, pp. 307–333, 2001.
  * [66] F. de Prenter, C. Leherenfeld, and A. Massing, “A note on the stability parameter in Nitsche’s method for unfitted boundary value problems,” Computers & Mathematics with Applications, 2018.
  * [67] C. Maple, “Geometric design and space planning using the marching squares and marching cube algorithms,” in 2003 International Conference on Geometric Modeling and Graphics, 2003. Proceedings, pp. 90–95, 2003.
  * [68] A. H. Gebremedhin, F. Manne, and A. Pothen, “What Color Is Your Jacobian? Graph Coloring for Computing Derivatives,” SIAM Rev., vol. 47, no. 4, pp. 629–705, 2005.
  * [69] E. Vinarcik, High Integrity Die Casting Processes.  Wiley, 2002.
  * [70] H. Ammar, A. Samuel, and F. Samuel, “Porosity and the fatigue behavior of hypoeutectic and hypereutectic aluminum–silicon casting alloys,” International Journal of Fatigue, vol. 30, no. 6, pp. 1024 – 1035, 2008.
  * [71] H. Mayer, M. Papakyriacou, B. Zettl, and S. Stanzl-Tschegg, “Influence of porosity on the fatigue limit of die cast magnesium and aluminium alloys,” International Journal of Fatigue, vol. 25, no. 3, pp. 245 – 256, 2003.
  * [72] J. Z. Yi, Y. X. Gao, P. D. Lee, H. M. Flower, and T. C. Lindley, “Scatter in fatigue life due to effects of porosity in cast A356-T6 aluminum-silicon alloys,” Metallurgical and Materials Transactions A, vol. 34, no. 9, p. 1879, 2003.
  * [73] G. Nicoletto, G. Anzelotti, and R. Konečná, “X-ray computed tomography vs. metallography for pore sizing and fatigue of cast Al-alloys,” Procedia Engineering, vol. 2, no. 1, pp. 547 – 554, 2010.  Fatigue 2010.
  * [74] M. Monavari, The Finite Cell Method for the analysis of die casts based on X-ray computed tomography.  Master’s thesis, Technische Universität München, Munich, 2011\.
  * [75] 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, vol. 28, no. 4, pp. 221–228, 2015.
  * [76] M. Würkner, S. Duczek, H. Berger, H. Köppe, U. Gabbert, E. Carrera, and G. Kulikov, “A software platform for the analysis of porous die-cast parts using the finite cell method,” in Analysis and Modelling of Advanced Structures and Smart Systems, pp. 327–341, Springer Singapore, 2018.
  * [77] Z. Yang, M. Ruess, S. Kollmannsberger, A. Düster, and E. Rank, “An Efficient Integration Technique for the Voxel-Based Finite Cell Method,” International Journal for Numerical Methods in Engineering, vol. 91, no. 5, pp. 457–471, 2012.
  * [78] M. A. Heroux, R. A. Bartlett, V. E. Howle, R. J. Hoekstra, J. J. Hu, T. G. Kolda, R. B. Lehoucq, K. R. Long, R. P. Pawlowski, E. T. Phipps, A. G. Salinger, H. K. Thornquist, R. S. Tuminaro, J. M. Willenbring, A. Williams, and K. S. Stanley, “An overview of the Trilinos project,” ACM Transactions on Mathematical Software, vol. 31, no. 3, pp. 397–423, 2005.
  * [79] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method.  Dover Books on Mathematics Series, Dover Publications, Incorporated, 2012\.