# Space-time unfitted finite elements on moving explicit geometry
representations

Santiago Badia1,2 ,  Pere A. Martorell3,∗ and  Francesc Verdugo4

###### Abstract.

This work proposes a novel variational approximation of partial differential
equations on moving geometries determined by explicit boundary
representations. The benefits of the proposed formulation are the ability to
handle large displacements of explicitly represented domain boundaries without
generating body-fitted meshes and remeshing techniques. For the space
discretization, we use a background mesh and an unfitted method that relies on
integration on cut cells only. We perform this intersection by using clipping
algorithms. To deal with the mesh movement, we pullback the equations to a
reference configuration (the spatial mesh at the initial time slab times the
time interval) that is constant in time. This way, the geometrical
intersection algorithm is only required in 3D, another key property of the
proposed scheme. At the end of the time slab, we compute the deformed mesh,
intersect the deformed boundary with the background mesh, and consider an
exact transfer operator between meshes to compute jump terms in the time
discontinuous Galerkin integration. The transfer is also computed using
geometrical intersection algorithms. We demonstrate the applicability of the
method to fluid problems around rotating (2D and 3D) geometries described by
oriented boundary meshes. We also provide a set of numerical experiments that
show the optimal convergence of the method.

1 School of Mathematics, Monash University, Clayton, Victoria, 3800,
Australia.
2 Centre Internacional de Mètodes Numèrics a l’Enginyeria, Campus Nord, 08034,
Barcelona, Spain.
3 Department of Civil and Environmental Engineering, Universitat Politècnica
de Catalunya, Edifici C1, Campus Nord UPC, Jordi Girona 1-3, 08034 Barcelona,
Spain.
4 Department of Computer Science, Vrije Universiteit Amsterdam, De Boelelaan
1111, 1081 HV Amsterdam, The Netherlands.
∗ Corresponding author.
E-mails<EMAIL_ADDRESS>(SB<EMAIL_ADDRESS>(PM),
<EMAIL_ADDRESS>(FV)

Keywords: Unfitted finite elements, embedded finite elements, computer-aided
design, computational geometry, immersed boundaries, boundary representations,
space-time discretizations, large displacements.

##  1\. Introduction

Space-time formulations for finite element methods are valuable techniques for
solving transient numerical problems. Unlike standard time stepping schemes,
which employ different discretizations for space and time, space-time
formulations discretize the problem simultaneously in both space and time.
However, generating 4D space-time body-fitted meshes becomes extremely
challenging on moving complex geometries. There are no general-purpose tools
to generate 4D meshes. Besides, the re-meshing process due to moving domains
introduces projection errors when transferring between meshes.

The bottleneck of mesh generation can be addressed by employing unfitted FEMs.
Unfitted FEMs, also known as embedded or immersed FEMs, offers a solution that
eliminates the need for body-fitted mesh generation, relying instead on a
simple background grid, such as a uniform or adaptive Cartesian grid. For
large-scale simulations on distributed-memory platforms, one can replace
expensive (both in terms of computational time and memory) unstructured mesh
partitioning [1] by efficient tree-mesh partitioners with load balancing [2,
3]. Unfitted FEMs have gained increasing popularity in various applications,
including fluid-structure interaction (FSI) [4, 5, 6], fracture mechanics [7,
8], additive manufacturing [9, 10], and stochastic geometry problems [11].
Traditionally, unfitted problems utilize level sets to describe geometries.
However, recent work [12, 13] has extended unfitted discretizations to
simulate geometries described by computer-aided design (CAD) models, i.e.,
explicit geometry representations.

The small cut cell problem, commonly discussed in the literature [14], is a
significant limitation of unfitted FEMs. The intersection between the physical
domain and the background cells can become arbitrarily small, thereby giving
rise to ill-conditioning issues. While several authors have attempted to
address this problem, only a few techniques have demonstrated robustness and
optimal convergence. The ghost penalty methods [15], employed within the so-
called CutFEM framework [16], represent one such approach. Alternatively, cell
agglomeration techniques have emerged as viable options to ensure robustness
in the presence of cut cells, naturally applied on discontinuous Galerkin (DG)
methods [17]. Extensions to the $\mathcal{C}^{0}$ Lagrangian finite element
(FE) have been introduced in [18], while mixed methods have been explored in
[19], where the aggregated finite element method (AgFEM) term was coined.
AgFEM exhibits good numerical qualities, including stability, bounds on
condition numbers, optimal convergence, and continuity concerning data.
Distributed implementations have been exploited in [20, 21], while AgFEM has
also been extended to $h$-adaptive meshes [22] and higher-order FE [23]. In
[24], a novel technique combining ghost penalty methods with AgFEM was
proposed, offering reduced sensitivity to stabilization parameters compared to
standard ghost methods.

In the body-fitted case, one can consider variational space-time formulations
[25, 26, 27] to avoid the generation of 4D meshes. In these formulations, a
body-fitted mesh is extruded using a geometric mapping technique to represent
the temporal evolution of the boundary displacements. However, when dealing
with large displacements, the geometric mapping process becomes ill-posed,
often necessitating re-meshing. Similar challenges arise in body-fitted
arbitrary Lagrangian Eulerian (ALE) [28, 29] schemes, where frequent re-
meshing becomes necessary due to large changes in topology. Consequently, none
of these formulations are robust when confronted with large displacements.
Furthermore, the re-meshing process introduces bottlenecks within the
simulation procedure and incurs in projection errors when transferring between
meshes.

Despite the notable advantages of unfitted FEMs, its application in the space-
time domain [30, 31] is currently limited to implicit level set geometrical
representations (in which 4D geometrical treatment is attainable), 2D explicit
representations (in which 3D geometrical tools, e.g., in [12], can readily be
used for space-time) or piecewise constant approximations of the geometry in
time [32]. This limitation stems from the complexity involved in developing 4D
geometrical tools for unfitted methods.

The present work introduces a novel variational space-time formulation for
unfitted FEMs that eliminates the need for complex 4D geometrical algorithms.
Instead, our formulation relies on time-extruded FE spaces. Furthermore, we
incorporate an _exact_ inter-slab transfer mechanism through an intersection
algorithm, which removes projection errors between time slabs. The main
contributions of our work are as follows:

  * •

We develop space-time formulations for unfitted FEMs on moving 3D geometries
described explicitly.

  * •

We propose to pull back the problem to a reference configuration that is
constant in each time slab to avoid 4D geometrical algorithms.

  * •

We propose and implement an exact time-slab mechanism transfer that relies on
3D intersection algorithms.

  * •

We compute the solution of transient problems on complex unfitted domains with
large displacements.

The outline of this paper is as follows. In Section 2, we introduce the
geometry description and the space-time FE spaces employed in this work. In
Section 3, we present the variational formulation through a model problem. We
also define the integration measures and the inter-slab transfer mechanism. In
Section 4, we describe the intersection algorithm for time slab transfer.
Then, in Section 5, we present numerical results for space and time
convergence, condition number tests, and a numerical example of a moving
domain. Finally, in Section 6, we present our conclusions and future work.

##  2\. Space time unfitted finite element method

In this section, we introduce the geometry description by utilizing a
geometrical map defined in a space-time FE space. We also define the unfitted
space-time FE spaces that are necessary to support our proposed formulation.

###  2.1. Geometry description for moving domains

Let us consider an initial Lipschitz domain $\Omega_{0}$ in a description
suitable for unfitted FEM. Following the methodology described in [12], we
assume that $\Omega_{0}$ is represented by a parameterized surface, e.g., a
stereolithography (STL) mesh or a CAD model. For simplicity, we consider an
oriented surface mesh $\mathcal{B}_{h}^{0}$ where the interior corresponds to
the domain $\Omega_{0}$. Let $[0,T]$ be the time interval of interest.

We consider the coordinates of the mesh nodes of $\mathcal{B}_{h}$ to be time-
depedendent. Consequently, $\mathcal{B}_{h}(t)$ and its corresponding interior
$\Omega(t)$. This variation of coordinates can be represented by a map
$\boldsymbol{D}:\mathcal{B}_{h}(0)\times[0,T]\rightarrow\mathcal{B}_{h}(t)$
such that $\boldsymbol{D}(\cdot,0)$ is the identity.

Next, we define the space-time domain $Q=\\{x\in\Omega(t):t\in[0,T]\\}$, in
accordance with the nomenclature in [30], see Fig. 1. Additionally, we
partition the boundary of the domain $\partial\Omega(t)$ into Dirichlet and
Neumann boundaries, $\partial\Omega_{D}(t)$ and $\partial\Omega_{N}(t)$, resp.
This partition is such that
$\partial\Omega(t)=\partial\Omega_{D}(t)\cup\partial\Omega_{N}(t)$ and
$\partial\Omega_{D}(t)\cap\partial\Omega_{N}(t)=\emptyset$. The space-time
boundaries are defined as $\partial
Q_{*}\doteq\bigcup_{t\in[0,T]}\partial\Omega_{*}(t)\times\\{t\\}$ for
$*\in\\{N,D\\}$. Consequently, the boundary of $Q$ is given by $\partial
Q\doteq\Omega(0)\cup\Omega(T)\cup\partial Q_{D}\cup\partial Q_{N}$.

$\Omega(t)$
---
$\Omega_{\mathrm{art}}$
---
$\Omega(0)$
---
$t^{n+1}$
---
$Q^{n}_{\mathrm{art}}$
---
$Q$
---
$Q_{\mathrm{art}}$
---
$Q^{n}$
---
$\partial Q^{n}$
---
$Q^{n}_{\mathrm{art}}$
---
$Q^{n}_{\mathrm{art}}$
---
$t^{n}$
---
$\bar{\mathcal{T}}^{n}$
---
$x$
---
$y$
---
$t$
---
$y$
---
$x$
---
Figure 1. Representation of the space-time domain $Q$ embedded within an
artificial space-time domain $Q_{\mathrm{art}}$. The spacial domains
$\Omega(t)$ and $\Omega_{\mathrm{art}}$ in 2D are extruded in the time
dimension. The solution is computed in each time-slab $J^{n}=(t^{n},t^{n+1})$
on the space-time domain $Q^{n}$, which is embedded in $Q_{\mathrm{art}}$.

In addition, we introduce an artificial spatial domain
$\Omega_{\mathrm{art}}$, such that $\Omega(t)\subset\Omega_{\mathrm{art}}$ for
all $t\in[0,T]$. This artificial domain can be a simple geometric shape, such
as a bounding box, that can conveniently be meshed with a regular grid, such
as a Cartesian mesh. Subsequently, we define the artificial space-time domain
$Q_{\mathrm{art}}\doteq\Omega_{\mathrm{art}}\times[0,T]$, such that $Q\subset
Q_{\mathrm{art}}$.

We define a partition in time, $\\{J^{n}\\}_{n=1}^{N}$ for $[0,T]$, in which
the time domain is divided into time slabs. The time slabs are defined as
$J^{n}=(t^{n},t^{n+1}),1<n<N$ where $t^{n}<t^{n+1},\forall n\in 1,...,N$.
Within each time slab, the time-step size is defined as
$\tau^{n}=t^{n+1}-t^{n}$ and the domain time-step size is defined as
$\tau=\max_{n=1,...,N}\tau^{n}$. In each time-slab, we define an artificial
space-time domain $Q^{n}_{\mathrm{art}}\doteq\Omega_{\mathrm{art}}\times
J^{n}$. In a time slab, the space-time domain is determined by $Q^{n}\doteq
Q^{n}_{\mathrm{art}}\cap Q$. Similarly, the space-time boundary is denoted as
$\partial Q^{n}_{\\{D,N\\}}$. We also introduce the notation
$\Omega^{n}\doteq\Omega(t^{n}),\ n=1,...,N+1$.

Finally, we define an undeformed space-time domain at each time slab as
$\hat{Q}^{n}\doteq\Omega^{n}\times J^{n}$, see Fig. 2. Here, we note that
$\Omega^{n}$ is the initial spatial domain of the time slab
$J^{n}=(t^{n},t^{n+1})$. In order to define the reference configuration, we
require a map $\boldsymbol{\varphi}^{n}_{h}:\hat{{Q}}^{n}\rightarrow{Q}^{n}$
that must satisfy

$\boldsymbol{\varphi}^{n}_{h}(\partial\Omega(t))=\boldsymbol{D}(t)(\mathcal{B}_{h}^{0}),$
(1)

i.e., respect the boundary position determined by $\boldsymbol{D}$.

$\hat{Q}^{n}$
---
$Q^{n}$
---
$Q^{n}_{\mathrm{art}}$
---
$\boldsymbol{\varphi}^{n}_{h}$
---
$Q^{n}_{\mathrm{art}}$
---
Figure 2.  Representation of the deformation map
$\boldsymbol{\varphi}^{n}_{h}$ resulting from extending the variation of the
surface map $\boldsymbol{D}$ in the time slab domain $Q^{n}$ (or in the time
slab artifical domain $Q^{n}_{\mathrm{art}}$). The undeformed configuration
$\hat{Q}^{n}$ will be used for the FE analysis.

Let $\bar{\mathcal{T}}_{\mathrm{art}}^{n}$ represent a simple partition of
$\Omega_{\mathrm{art}}$ at time $t^{n}$. Even though $\Omega_{\mathrm{art}}$
does not evolve in time, the partition $\bar{\mathcal{T}}_{\mathrm{art}}^{n}$
can differ across time slabs, e.g., when using adaptive mesh refinement (AMR)
techniques. We define the partition of $Q^{n}_{\mathrm{art}}$ as a Cartesian
product of $\bar{\mathcal{T}}_{\mathrm{art}}^{n}$ and $J^{n}$. Specifically,
$\mathcal{T}_{\mathrm{art}}^{n}\doteq\\{\bar{K}\times
J^{n}:\bar{K}\in\bar{\mathcal{T}}^{n}_{\mathrm{art}}\\}$. We can classify the
spatial artificial cells $K\in\bar{\mathcal{T}}^{n}_{\mathrm{art}}$ as
interior, cut, and exterior depending on the relative position concerning the
domain boundary $\partial\Omega^{n}$. Since there is a one-to-one map between
cells in $\bar{\mathcal{T}}^{n}_{\mathrm{art}}$ and
${\mathcal{T}}^{n}_{\mathrm{art}}$, it already provides a classification of
the cells in the space-time artificial mesh $\mathcal{T}^{n}_{\mathrm{art}}$,
which corresponds to the in-out classification for the extruded space-time
domain boundary $\partial\hat{Q}^{n}$.

In the context of unfitted FEM, we are interested in solving partial
differential equations using the active portion of
$\mathcal{T}_{\mathrm{art}}^{n}$, i.e., the cells touching $\Omega^{n}$ at
$t^{n}$. Thus, we remove the exterior cells
$\bar{\mathcal{T}}^{n}_{\mathrm{out}}\doteq\\{\bar{K}\in\bar{\mathcal{T}}^{n}_{\mathrm{art}}:\bar{K}\cap\Omega^{n}=\emptyset\\}$
from the artificial spatial partition $\bar{\mathcal{T}}^{n}_{\mathrm{art}}$,
leading to an active spatial partition
$\bar{\mathcal{T}}^{n}\doteq\bar{\mathcal{T}}^{n}_{\mathrm{art}}\setminus\bar{\mathcal{T}}^{n}_{\mathrm{out}}$.
From the spatial active partition, we compute the active space-time partition
$\mathcal{T}^{n}\doteq\\{\bar{K}\times
J^{n}:\bar{K}\in\bar{\mathcal{T}}^{n}\\}$ by extrusion.

###  2.2. Space-time finite element spaces

In this section, we define space-time unfitted FE spaces for the
discretization of the model problem on moving domains. However, a main
difference with respect to [30], the original problem on $Q^{n}$ will be re-
state in the undeformed domain $\hat{Q}^{n}$. As a result, in this section, we
only need a construction that works on time-extruded domains. Thus, the active
mesh $\mathcal{T}^{n}$ is the one that intersects $\Omega^{n}$ at $t^{n}$, as
defined above. We note that $\mathcal{T}^{n}$ may change across time slabs in
evolving domains.

Firstly, we define the FE space $\mathcal{X}^{n}_{h}$ in the spatial domain.
Then, for each degree of freedom (DOF) in $\mathcal{X}_{h}$, we introduce a 1D
FE basis $\mathcal{Y}^{n}_{\tau}$. Thus, we construct the space-time FE space
composed as the tensor product of these two spaces,
$\mathcal{V}^{n}_{h}\doteq\mathcal{X}^{n}_{h}\otimes\mathcal{Y}^{n}_{\tau}$.

In each space-time cell $T^{n}=\bar{T}^{n}\times J^{n}\in\mathcal{T}^{n}$, we
can define the FE local interpolation, which consists of tensor product DOFs
$\Sigma\doteq\Sigma_{X}\otimes\Sigma_{Y}$ and shape functions
$\Phi\doteq\Phi_{X}\otimes\Phi_{Y}$. Specifically, any shape function
$\phi\in\Phi$ can be expressed as
$\phi^{(\alpha_{X},\alpha_{Y})}(\boldsymbol{x},t)=\phi^{\alpha_{X}}_{X}(\boldsymbol{x})\otimes\phi^{\alpha_{Y}}_{Y}(t)$.

To address the ill-conditioning issues in unfitted FEM, we will utilize as a
model example the AgFEM, even though the proposed numerical framework can
readily be used for other techniques, e.g., ghost penalty stabilization (see
[30] for space-time unfitted formulations). The AgFEM eliminates problematic
DOFs by constraining them to well-posed DOFs using a discrete extension
operator $\mathcal{E}$. We define the spatial extension operator
$\bar{\mathcal{E}}:\bar{\mathcal{V}}_{h,\mathrm{in}}\mapsto\bar{\mathcal{V}}_{h}$
between the the spatial interior $\bar{\mathcal{V}}_{h,\mathrm{in}}$ and
active $\bar{\mathcal{V}}_{h}$ FE spaces . The spatial aggregated finite
element (AgFE) space is defined as
$\bar{\mathcal{V}}_{h,\mathrm{ag}}\doteq\bar{\mathcal{E}}(\bar{\mathcal{V}}_{h,\mathrm{in}})$

We define a slab-wise space-time discrete operator between
$\mathcal{V}^{n}_{h,\mathrm{in}}$ and $\mathcal{V}^{n}_{h}$. In each time
slab, we aggregate the cut cells in the spatial active mesh
$\bar{\mathcal{T}}^{n}$ using the aggregation techniques from the AgFEM. For a
given time slab $J^{n}$, the space-time extension operator
$\mathcal{E}^{n}:\mathcal{V}^{n}_{h,\mathrm{in}}\mapsto\mathcal{V}^{n}_{h}$ is
defined as

$\mathcal{E}^{n}(\mathcal{V}^{n}_{h,\mathrm{in}})=\mathcal{E}^{n}(\bar{\mathcal{V}}^{n}_{h,\mathrm{in}}\otimes\mathcal{Y}^{n}_{\tau})=\bar{\mathcal{E}}^{n}(\bar{\mathcal{V}}^{n}_{h,\mathrm{in}})\otimes\mathcal{Y}^{n}_{\tau},$
(2)

where
$\bar{\mathcal{E}}^{n}:\bar{\mathcal{V}}^{n}_{h,\mathrm{in}}\mapsto\bar{\mathcal{V}}^{n}_{h}$
represents the spatial extension operator, and $\mathcal{Y}_{\tau}^{n}$ the FE
space in time at the time slab $J^{n}$. We define the space-time AgFE space on
the time-slab. Now, we can define the global AgFE space as
$\mathcal{V}_{h,\mathrm{ag}}\doteq\mathcal{V}_{h,\mathrm{ag}}^{1}\times\dots\times\mathcal{V}_{h,\mathrm{ag}}^{N}$.
Note that no continuity is imposed across time slabs, i.e., the discrete time
space is discontinuous.

In ghost penalty formulations, one can readily use the non-aggregated spaces
on the active mesh and add stabilization terms (see, e.g., [33]) to make the
problem robust with respect to cut locations. Thus, the space-time FE spaces
to be used in this case is simply $\mathcal{V}_{h}^{n}$.

###  2.3. Extension of the deformation map

In the case in which the domain is represented in time by its explicit
boundary representation and its boundary displacement $\boldsymbol{D}(t)$, we
must extend $\boldsymbol{D}(t)$ defined on $\mathcal{B}(t)$ to the space-time
domain $\hat{Q}^{n}$ of each time slab $J^{n}$. For this purpose, we solve the
following linear elasticity problem in $\hat{Q}^{n}$, even though other
extension operators could be considered. Find
$\hat{\boldsymbol{u}}\in\mathbb{R}^{d}$ such that,

$\begin{cases}\hat{\boldsymbol{\nabla}}_{x}\cdot\boldsymbol{\sigma}(\hat{\boldsymbol{u}})=0&\mathrm{in}\
\hat{Q}^{n},\\\ \hat{\boldsymbol{u}}=\hat{\boldsymbol{u}}_{D}^{n}&\mathrm{on}\
\partial\hat{Q}^{n}_{D},\\\ \hat{\boldsymbol{u}}=0&\mathrm{on}\
\partial\hat{Q}^{n}_{D_{0}},\\\
\mathbf{n}_{x}\cdot\boldsymbol{\sigma}(\hat{\boldsymbol{u}})=0&\mathrm{on}\
\partial\hat{Q}^{n}_{N}.\end{cases}$

where $\boldsymbol{\sigma}$ is the stress tensor,
$\hat{\boldsymbol{\nabla}}_{x}$ is the spatial gradient, $\mathbf{n}_{x}$ is
the spatial component of outward normal vector to $\partial\hat{Q}^{n}$, and

$\hat{\boldsymbol{u}}_{D}^{n}(t)=\left(\boldsymbol{D}(t)-\boldsymbol{D}(t^{n})\right)\circ\left(\boldsymbol{D}(t^{n})\right)^{-1}$

is the Dirichlet boundary condition on
$\partial\hat{Q}^{n}_{D}=\mathcal{B}_{h}(t^{n})\times J^{n}$, which imposes
the deformation near the unfitted boundary (see Fig. 3). On the boundary of
the artificial domain and at the initial spatial domain,
$\partial\hat{Q}^{n}_{D_{0}}=\partial\Omega_{\mathrm{art}}\times
J^{n}\cup\Omega^{n}\times\\{t^{n}\\}$, the deformation is fixed to zero. The
Neumann boundary is defined on
$\partial\hat{Q}^{n}_{N}=\Omega^{n}\times\\{t^{n+1}\\}$.

$\boldsymbol{D}(t^{n})$
---
$\boldsymbol{D}(t^{n+1})$
---
$t^{n+1}$
---
$t^{1}=0$
---
$t^{n}$
---
$\boldsymbol{\varphi}^{n}_{h}(t^{n+1})$
---
$\mathcal{B}_{h}(0)$
---
$\partial\Omega_{I}^{n+1}$
---
$\partial\Omega_{I}^{n}$
---
$\Omega(0)$
---
$\Omega^{n}$
---
$\tilde{\Omega}^{n+1}$
---
Figure 3. Representation of the deformation map $\boldsymbol{D}$ and the
extended map $\boldsymbol{\varphi}_{h}^{n}$ in the time slab
$J^{n}=(t^{n-1},t^{n})$. This example represents the external domain. The
spatial deformation map
$\boldsymbol{D}(t):\mathcal{B}_{h}(0)\mapsto\mathcal{B}_{h}(t),\ t\in[0,T]$,
is defined on the surface mesh, and the extended map
$\boldsymbol{\varphi}_{h}^{n}(t):\Omega^{n}\mapsto\tilde{\Omega}(t),\ t\in
J^{n}$, is defined on the spatial domain. At $t^{n+1}$, the description of the
inner boundary $\partial\Omega_{I}^{n+1}$ may differ across maps:
$\boldsymbol{\varphi}_{h}^{n}(t^{n+1})(\partial\Omega_{I}^{n})\approx\boldsymbol{D}(t^{n+1})(\mathcal{B}_{h}(0))$,
where $\partial\Omega_{I}^{n}=\boldsymbol{D}(t^{n})(\mathcal{B}_{h}(0))$. The
approximation error decreases with the spatial discretization size $h$.

This problem is approximated with a continuous Galerkin (CG) method and weak
imposition of the Dirichlet boundary conditions on the unfitted boundary
$\partial\hat{Q}_{D}$ using Nitsche’s method [34]. In the weak formulation, we
find $u\in\mathcal{V}^{n}_{h,\mathrm{ag}}$ such that,

$a(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})=l(\hat{\boldsymbol{v}}),\quad\forall\hat{\boldsymbol{v}}\in\mathcal{V}^{n}_{h,\mathrm{ag}}$

where the bilinear form and the linear form are defined as,

$a(\hat{\boldsymbol{u}},\hat{\boldsymbol{v}})=\int_{\hat{Q}^{n}}\boldsymbol{\varepsilon}(\hat{\boldsymbol{v}}):\boldsymbol{\sigma}(\boldsymbol{\varepsilon}(\hat{\boldsymbol{u}}))d\hat{\boldsymbol{x}}dt+\int_{\partial\hat{Q}_{D}}\left(\tau_{D}\hat{\boldsymbol{v}}\cdot\hat{\boldsymbol{u}}-v\cdot\left(\mathbf{n}_{x}\cdot\boldsymbol{\sigma}(\boldsymbol{\varepsilon}(\hat{\boldsymbol{u}}))\right)\right)d\hat{\boldsymbol{x}}dt,$

and

$l(\hat{\boldsymbol{v}})=\int_{\partial\hat{Q}_{D}}\left(\hat{\boldsymbol{v}}\cdot\hat{\boldsymbol{u}}_{D}-\mathbf{n}_{x}\cdot\boldsymbol{\sigma}(\boldsymbol{\varepsilon}(\hat{\boldsymbol{v}}))\cdot\hat{\boldsymbol{u}}_{D}\right)d\hat{\boldsymbol{x}}dt,$

resp., where $\boldsymbol{\varepsilon}$ is the symmetric (spatial) gradient
operator.

From the approximated solution
$\hat{\boldsymbol{u}}^{n}_{h}\in\mathcal{V}^{n}_{h,\mathrm{ag}}$ we extract
the map
$\boldsymbol{\varphi}^{n}_{h}(\hat{\boldsymbol{x}},t)=\hat{\boldsymbol{x}}+\hat{\boldsymbol{u}}^{n}_{h}(\hat{\boldsymbol{x}},t)$.
One can alternatively use a ghost penalty stabilization to solve the
elasticity problem. Since the deformation is imposed weakly, the deformation
$\hat{\boldsymbol{u}}^{n}_{h}$ is not equal to $\hat{\boldsymbol{u}}_{D}^{n}$
on $\mathcal{B}_{h}(t)$. Thus, the deformed space-time domain $Q^{n}$ is
approximated by $\tilde{Q}^{n}=\boldsymbol{\varphi}^{n}_{h}(\hat{Q}^{n})$. In
any case, the geometrical error is expected to decrease with the mesh size
$h$. For given time $t\in J^{n}$, we can extract the spatial map
$\boldsymbol{\varphi}^{n}_{h}(t):\Omega^{n}\mapsto\tilde{\Omega}(t)$ (see Fig.
3). Here, the spatial domain is also approximated,
$\Omega(t)\approx\tilde{\Omega}(t)$.

We assume that the computed map is one-to-one, i.e.,
$\mathrm{det}(\boldsymbol{\nabla}\boldsymbol{\varphi}_{h}^{n})>0$ at all
times. Assuming that $\boldsymbol{D}$ is such that it admits an extension that
is diffeomorphic, this requirement can be attained for small enough time steps
in the time-discrete problem. One can also move back to the Cartesian mesh
after several time steps, as soon as the map remains bijective.

###  2.4. Extended active mesh

As discussed above, in the case in which the deformation map has to be
computed, $\boldsymbol{\varphi}^{n}_{h}$ is not equal to $\boldsymbol{D}$ on
$\mathcal{B}_{h}(t)$, since unfitted methods usually make use of weak
imposition of boundary conditions.

At the end of the time slab $J^{n}$, we can define the solution in the
approximated domain $\tilde{Q}^{n}$ as

${u}_{h}^{n+1,-}(\boldsymbol{x})=\hat{u}_{h}^{n}(\boldsymbol{\varphi}^{n}_{h}(t^{n+1})^{-1}(\boldsymbol{x}),t^{n+1}),\quad\forall\boldsymbol{x}\in\tilde{\Omega}^{n+1},$

which is a function defined on the undeformed domain $\Omega^{n}$; since we
integrate the forms on $\hat{Q}^{n}$, the inverse map is never computed in
practice. On the other side, with the method proposed below, we have to
compute an inter-slab integral on $\Omega^{n+1}$ that involves
$u^{n+1,-}_{h}$.

In order for $u_{h}^{n+1,-}$ to be defined on $\Omega^{n+1}$, we proceed as
follows. First, we extend the active mesh by
$\bar{\mathcal{T}}_{\mathrm{ext}}^{n}=\bar{\mathcal{T}}^{n}\cup\\{K\in\bar{\mathcal{T}}^{n}_{\mathrm{art}}:\boldsymbol{\varphi}^{n}(t^{n+1})(K)\cap\Omega^{n+1}\neq\emptyset\\}$.
Assuming that the geometrical map is one-to-one, the inverse
$\boldsymbol{\varphi}^{n}_{h}(t^{n+1})^{-1}$ is well-defined on
$\Omega^{n+1}$.

To accommodate a FE space $\bar{\mathcal{V}}_{h,\mathrm{ext}}^{n}$ in
$\bar{\mathcal{T}}_{\mathrm{ext}}^{n}$, one can consider the modification of
the extension operator. In the AgFEM framework, we can simply redefine the
extension operator
$\bar{\mathcal{E}}^{n}_{\mathrm{ext}}:\bar{\mathcal{V}}^{n}_{h,\mathrm{in}}\mapsto\bar{\mathcal{V}}_{h,\mathrm{ext}}^{n}$.
In non-aggregated methods like CutFEM, we can utilize the same extension on
the solution
$\bar{\mathcal{E}}^{n}_{\mathrm{ext}}(\bar{\mathcal{V}}^{n}_{h})$.

We note that the implementation of the extended triangulation and FE space can
be computed a posteriori, on demand, whenever it is needed. One can mark the
additional cells that are needed for the extension, and then extend the
aggregates to these cells.

##  3\. Variational formulation on a model problem

In this section, we establish the space-time variational formulation by
employing a model problem, specifically the heat equation. Although we use the
heat equation for demonstration purposes, it is essential to note that a
similar approach can be applied to other PDEs.

###  3.1. Weak formulation

In order to define the space-time variational formulation, we first define the
convection-diffusion equation in a space-time domain $Q$, as find $u$ such
that,

$\begin{cases}\partial_{t}u+(\boldsymbol{w}\cdot\boldsymbol{\nabla}_{x})u-\mu\Delta_{x}u=f\quad&\text{in
}Q,\\\ u=u_{D}\quad&\text{on }\partial Q_{D},\\\
\mu\mathbf{n}_{x}\cdot\boldsymbol{\nabla}_{x}u=g_{N}\quad&\text{on }\partial
Q_{N},\\\ u=u_{0}\quad&\text{on }Q(0).\end{cases}$ (3)

where $\mu$ is the diffusion coefficient, $\boldsymbol{w}$ the advection
velocity, $f$ the source term, $u_{D}$ the Dirichlet boundary condition, and
$g_{N}$ boundary flux on $\partial Q_{N}$. $\boldsymbol{\nabla}_{x}$ and
$\Delta_{x}$ denote the spatial gradient and spatial Laplacian, resp. Well-
posedness requires that $\boldsymbol{w}\cdot\mathbf{n}_{x}+n_{t}\geq 0$ on the
Neumann boundary $\partial Q_{N}$. Here, $\mathbf{n}=(\mathbf{n}_{x},n_{t})$
is the outward normal to $\partial Q$.

We discretize this problem with the AgFEM in space (or a ghost penalty
stabilization ) and a DG method in time. We weakly impose the Dirichlet
boundary conditions using Nistche’s method [34]. Since the coupling between
time slabs respects causality, we analyze the problem on a single time slab,
assuming we know the solution of the previous one (see also [30]). To analyze
each time slab $J^{n}=(t^{n},t^{n-1})$, we define the problem in the reference
domain $\hat{Q}^{n}=(\boldsymbol{\varphi}^{n}_{h})^{-1}(\tilde{Q}^{n})$ as:
find $\hat{u}\in\mathcal{V}^{n}_{h,\mathrm{ag}}$ such that

$B_{h}^{n}(\hat{u},\hat{v})=L_{h}^{n}(\hat{v}),\quad\forall\hat{v}\in\mathcal{V}^{n}_{h,\mathrm{ag}},$
(4)

with $u=\hat{u}\circ(\boldsymbol{\varphi}^{n}_{h})^{-1}$. The bilinear form
reads as:

$\displaystyle B_{h}^{n}(\hat{u},\hat{v})$
$\displaystyle=\int_{\hat{Q}^{n}}\hat{v}\partial_{t}^{n}\hat{u}|J_{\tilde{Q}^{n}}|d\hat{\boldsymbol{x}}dt+\int_{\Omega^{n}}\hat{v}(t^{n})\hat{u}(t^{n})|J_{\Omega^{n}}|d\hat{\boldsymbol{x}}+a_{h}(\hat{u},\hat{v}),$
(5) $\displaystyle a_{h}(\hat{u},\hat{v})$
$\displaystyle=\int_{\hat{Q}^{n}}\left(\mu\boldsymbol{\nabla}^{n}_{x}\hat{u}\cdot\boldsymbol{\nabla}^{n}_{x}\hat{v}\
+\
\hat{v}(\boldsymbol{w}\circ\boldsymbol{\varphi}_{h}^{n}\cdot\boldsymbol{\nabla}^{n}_{x})\hat{u}\right)|J_{\tilde{Q}^{n}}|d\hat{\boldsymbol{x}}dt$
$\displaystyle+\int_{\partial\hat{Q}^{n}_{D}}\left(\beta_{h}\hat{v}\hat{u}-\hat{v}\left(\mathbf{n}^{n}_{x}\cdot\mu\boldsymbol{\nabla}^{n}_{x}\hat{u}\right)-\left(\mathbf{n}^{n}_{x}\cdot\mu\boldsymbol{\nabla}^{n}_{x}\hat{v}\right)\hat{u}\right)|J_{\partial\tilde{Q}^{n}_{D}}|d\hat{\boldsymbol{x}}dt,$

and the linear form reads as,

$\displaystyle L^{n}_{h}(\hat{v})$
$\displaystyle=\int_{\Omega^{n}}\hat{v}(t^{n})\hat{u}^{n-1}(t^{n})|J_{\Omega^{n}}|d\hat{\boldsymbol{x}}+l_{h}(\hat{v}),$
(6) $\displaystyle l_{h}(\hat{v})$
$\displaystyle=\int_{\hat{Q}^{n}}\hat{v}\left(f\circ\boldsymbol{\varphi}_{h}^{n}\right)|J_{\tilde{Q}^{n}}|d\hat{\boldsymbol{x}}dt+\int_{\partial\hat{Q}^{n}_{N}}\hat{v}\left(g_{N}\circ\boldsymbol{\varphi}_{h}^{n}\right)|J_{\partial\tilde{Q}^{n}_{N}}|d\hat{\boldsymbol{x}}dt$
$\displaystyle+\int_{\partial\hat{Q}^{n}_{D}}\left(\beta_{h}\hat{v}\left(u_{D}\circ\boldsymbol{\varphi}_{h}^{n}\right)-\left(\mathbf{n}^{n}_{x}\cdot\mu\boldsymbol{\nabla}^{n}_{x}\hat{v}\right)\left(u_{D}\circ\boldsymbol{\varphi}_{h}^{n}\right)\right)|J_{\partial\tilde{Q}^{n}_{D}}|d\hat{\boldsymbol{x}}dt,$

Let us define the norms used in [30] to prove stability and convergence
results, which we will computed in the numerical experiments. In [30] the
space-time accumulated DG norm of $\mathcal{V}_{h,\mathrm{ag}}$ is defined as
follows,

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{n,*}\doteq\|v^{n}(t^{n+1})\|^{2}_{L^{2}(\Omega^{n})}+\sum^{n}_{i=1}\|v^{i}(t^{i})-v^{i-1}(t^{i})\|^{2}_{L^{2}(\Omega^{i})}+c_{\mu}\int^{t^{n+1}}_{0}\|v\|^{2}_{\bar{\mathcal{V}}^{n}(h)}dt,$
(7)

where $c_{\mu}$ is the coercivity constant and
$\bar{\mathcal{V}}^{n}(h)\doteq\bar{\mathcal{V}}^{n}_{h,\mathrm{ag}}+H^{2}(\Omega(t))$
is the norm of the FE space at a time step $t$ given by

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\bar{\mathcal{V}}^{n}(h)}\doteq\mu\|\boldsymbol{\nabla}v\|^{2}_{L^{2}(\Omega(t))}+\sum_{\bar{T}\in\bar{\mathcal{T}}_{h,\mathrm{act}}}\beta_{\bar{T}}\|v\|^{2}_{L^{2}(\bar{T}\cap\partial\Omega_{D}(t))}+\sum_{\bar{T}\in\bar{\mathcal{T}}^{n}_{h,\mathrm{act}}}\mu
h^{2}_{\bar{T}}\|v\|^{2}_{H^{2}(\bar{T}\cap\Omega(t))}.$ (8)

The proof of this stability result follows the ideas in [30] in the case in
which the deformation map $\boldsymbol{\varphi}_{h}^{n}(t)$ is equal to
$\boldsymbol{D}(t)$ on $\Gamma(t)$. Otherwise, the analysis is more technical
and would require to use ideas similar to the ones in [32, 35]. The analysis
therein assumes a constant deformation map and a level-set description of the
domain. In our case, the deformation map can be of higher-order (time-
dependent) and the domain is represented by its boundary. We will leave this
analysis for future work.

In the variational form, the Dirichlet boundary conditions are imposed weakly
with the Nistche method, with a penalty term $\beta_{h}$ that depends on the
spatial cell size $h$. The initial value is also imposed weakly with DG in
time, where the jump is given by $[\hat{u}(t^{n})-\hat{u}^{n-1}(t^{n})]$. The
evaluation of the solution of the previous time slab $\hat{u}^{n-1}(t^{n})$ in
$\Omega^{n}$ requires special attention since it is computed using a different
discretization. Further discussion is provided in Sect. 3.2 and Sect. 4.

The derivatives are projections of the space-time gradient into space and
time, $\partial^{n}_{t}=(\boldsymbol{\nabla}^{n})_{t},$ and
$\boldsymbol{\nabla}^{n}_{x}=(\boldsymbol{\nabla}^{n})_{x},$ resp. They are
obtained by transporting the space-time gradient to the deformed domain
$\boldsymbol{\nabla}^{n}=\mathbf{F}^{-T}\hat{\boldsymbol{\nabla}},$ where
$\mathbf{F}=\boldsymbol{\nabla}\boldsymbol{\varphi}_{h}^{n}$. The operators
$(~{})_{x}$ and $(~{})_{t}$ represent the projection to the space and time
directions, resp. The boundary normal is also transported to the deformed
domain. It is computed as follows,

$\mathbf{n}^{n}_{x}=\left(\frac{\mathbf{F}^{-T}\mathbf{n}_{\partial\hat{Q}^{n}}}{\|\mathbf{F}^{-T}\mathbf{n}_{\partial\hat{Q}^{n}}\|}\right)_{x},$
(9)

where $\mathbf{n}_{\partial\hat{Q}^{n}}$ is the normal vector in the
undeformed space-time domain.

The integration measures to define the domain change are defined by the
jacobians of the deformed domain. The integration measure change of the space-
time domain $\tilde{Q}^{n}$ is given by $J_{Q^{n}}=J$, where
$J=det(\mathbf{F})$. The initial boundary Jacobian is defined as
$J_{\Omega^{n}}=det(\mathbf{F}_{x}(t^{n}))$. In the given case
$J_{\Omega^{n}}=1$. Note that for small deformations, where
$\min_{\Omega}(J)\approx 1$, a simpler approach can be considered by assuming
a deformed initial state, e.g.,
$\Omega^{n}=\boldsymbol{\varphi}^{n-1}_{h}(t^{n})(\Omega^{n-1})$. In this
situation, the evaluation $\hat{u}_{h}^{n-1}(t^{n})$ does not require the
change of reference space described in Sect. 3.2. This approach is equivalent
to a time slab with more than one cell in time and continuous FE spaces in
time could be considered within this macro-cell.

The pullback of the area differential form to the reference domain is
expressed as

$da_{\partial\tilde{Q}^{n}}=J\sqrt{\mathbf{n}^{T}_{\partial\hat{Q}^{n}}\cdot\mathbf{C}^{-1}\cdot\mathbf{n}_{\partial\hat{Q}^{n}}}~{}da_{\partial\hat{Q}^{n}}$
(10)

where $\mathbf{C}=\mathbf{F}^{T}\mathbf{F}$ (see, e.g., [36]).

The space-time gradients are transported to the deformed domain,
$\boldsymbol{\nabla}^{n}=\mathbf{F}^{-T}\hat{\boldsymbol{\nabla}}$. The map
gradient, $\mathbf{F}=\boldsymbol{\nabla}\boldsymbol{\varphi}_{h}^{n}$,
contains the terms as follows:

$\mathbf{F}=\left[\begin{matrix}\mathbf{F}_{x}&0\\\
\partial_{t}\boldsymbol{\varphi}_{h,x}^{n}&1\end{matrix}\right]=\left[\begin{matrix}\mathbf{F}_{x}&0\\\
\mathbf{w}^{T}&1\end{matrix}\right],$ (11)

where $\mathbf{F}_{x}$ is the space gradient and
$\mathbf{w}^{T}=\partial_{t}\boldsymbol{\varphi}_{h,x}^{n}$ the deformation
velocity. Then, the inverse gradient is given by

$\mathbf{F}^{-T}=\left[\begin{matrix}\mathbf{F}_{x}^{-1}&0\\\
-\mathbf{w}^{T}\mathbf{F}_{x}^{-1}&1\end{matrix}\right].$ (12)

Then, by decomposing the space-time gradient into space and time derivatives,
we obtain

$\left[\begin{matrix}\boldsymbol{\nabla}^{n}_{x}\\\
\partial^{n}_{t}\end{matrix}\right]=\mathbf{F}^{-T}\left[\begin{matrix}\hat{\boldsymbol{\nabla}}_{x}\\\
\hat{\partial}_{t}\end{matrix}\right]=\left[\begin{matrix}\mathbf{F}^{-1}_{x}\hat{\boldsymbol{\nabla}}_{x}\\\
\hat{\partial}_{t}-\mathbf{w}^{T}\mathbf{F}^{-1}_{x}\hat{\boldsymbol{\nabla}}_{x}\end{matrix}\right],$
(13)

which already recovers the derivative terms used in ALE formulations.

###  3.2. Inter-slab integration

In the formulation (5)-(6), a DG method is used in time, where the initial
value at the time slab is imposed weakly through an inter-time slab jump
$[\hat{u}^{n}(t^{n})-\hat{u}^{n-1}(t^{n})]$. However, integrating this jump is
not straightforward, since $\hat{u}^{n}(t^{n})$ and $\hat{u}^{n-1}(t^{n})$ in
(6) are expressed in different discrete spaces (and meshes).

To address this evaluation, an intermediate unfitted discretization
$\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ is introduced.
$\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ is a partition of $\Omega^{n}$ that
results of intersecting the embedded discretization
$\mathcal{T}_{\mathrm{cut}}^{n}$ of $\Omega^{n}$ and
$\bar{\mathcal{T}}^{n}_{-}=\boldsymbol{\varphi}^{n-1}_{h}(t^{n})(\bar{\mathcal{T}}^{n-1})$,
i.e., $\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ results from intersecting
$\bar{\mathcal{T}}^{n}$, $\bar{\mathcal{T}}^{n}_{-}$ and $\Omega^{n-1}$. Each
cell $K_{\mathrm{int}}\in\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ has an
injective map to ${K}^{n-1}\in\bar{\mathcal{T}}^{n-1}$ and
${K}^{n}\in\bar{\mathcal{T}}^{n}$. A representation of
$\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ is depicted in Fig. 4, and its
construction is detailed in Sect. 4.

(a) $\boldsymbol{\varphi}^{n-1}_{h}(t^{n})(\bar{\mathcal{T}}^{n-1})$

(b) $\bar{\mathcal{T}}_{\mathrm{int}}^{n}$

(c) $\bar{\mathcal{T}}^{n}$

Figure 4. Mesh sequence for solution transfer between time slabs. The solution
obtained in $\bar{\mathcal{T}}^{n-1}$ in (a) is then evaluated at
$\bar{\mathcal{T}}^{n}$ (c). Acting as a bridge, the intersected mesh
$\bar{\mathcal{T}}_{\mathrm{int}}^{n}=\boldsymbol{\varphi}^{n-1}_{h}(t^{n})(\bar{\mathcal{T}}^{n-1})\cap\bar{\mathcal{T}}^{n}\cap\Omega^{n}$
in (b) facilitates the evaluation by providing injective cell maps to both
active meshes.

The integration of the jump in the time slab interface is performed in
$\bar{\mathcal{T}}_{\mathrm{int}}^{n}$. Thus, we need to define the cell maps
from $\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ to $\bar{\mathcal{T}}^{n}$ and
$\bar{\mathcal{T}}^{n+1}$. For $\hat{u}^{n-1}(t^{n})$ we need a cell map from
$\hat{K}_{\mathrm{int}}\in\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ to
$\hat{K}_{1}\in\bar{\mathcal{T}}^{n}$ such that

$\psi^{-}=\left(\phi_{K_{1}}\right)^{-1}\circ\left(\boldsymbol{\varphi}^{n-1}\right)^{-1}\circ\left(\phi_{K_{\mathrm{int}}}\right),$
(14)

where $\phi_{K_{1}}$ maps from the reference space to the undeformed physical
space of $K_{1}\in\bar{\mathcal{T}}^{n}$ and $\phi_{K_{\mathrm{int}}}$ maps
from the reference space to the physical space of $K_{\mathrm{int}}$, see Fig.
5. In the case of $\hat{u}^{n}(t^{n})$ we need a cell map from
$\hat{K}_{\mathrm{int}}\in\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ to
$\hat{K}^{n}\in\bar{\mathcal{T}}^{n}$ such that

$\psi^{+}=\left(\phi_{K_{2}}\right)^{-1}\circ\left(\phi_{K_{\mathrm{int}}}\right).$
(15)

where $\phi_{K_{2}}$ maps from the reference space to the physical space of
$K_{2}\in\bar{\mathcal{T}}^{n}$. Now, we can numerically integrate the jump by
evaluating the following integral,

$\sum_{K_{\mathrm{int}}\in\bar{\mathcal{T}}_{\mathrm{int}}}\sum_{\hat{q}}\left(\hat{u}^{n}(t^{n})\circ\psi^{+}-\hat{u}^{n-1}(t^{n})\circ\psi^{-}\right)(v^{n}(t^{n})\circ\psi^{+}\left(\hat{q}\right))|J_{K_{\mathrm{int}}}|w_{K_{\mathrm{int}}},$
(16)

where $J_{K_{\mathrm{int}}}$, $\hat{q}$ and $w_{K_{\mathrm{int}}}$ are the
Jacobian, quadrature points and quadrature weights of $K_{\mathrm{int}}$,
resp.

$\hat{K}$
---
$\hat{K}_{\mathrm{int}}$
---
$K_{1}$
---
$K_{2}$
---
$\boldsymbol{\varphi}^{n-1}(K_{1})$
---
$K_{\mathrm{int}}$
---
$\phi_{K_{1}}$
---
$\phi_{K_{2}}$
---
$\boldsymbol{\varphi}^{n-1}$
---
$\hat{x}$
---
$\hat{y}$
---
$x$
---
$y$
---
$\phi_{K_{\mathrm{int}}}$
---
Figure 5. Representation of the cell maps used in the time slab interface
integration. The maps $\phi_{K_{1}}$ and $\phi_{K_{2}}$ send the reference
element $\hat{K}$ to the physical space of $K_{1}\in\bar{\mathcal{T}}^{n-1}$
and $K_{2}\in\bar{\mathcal{T}}^{n}$, resp. The map $\phi_{K_{\mathrm{int}}}$
sends the reference element $\hat{K}_{\mathrm{int}}$ to the physical space of
$K_{\mathrm{int}}\in\bar{\mathcal{T}}^{n-1}_{\mathrm{int}}$.

##  4\. Intersection algorithm for time slab transfer

In this section, we will present one of the key novelties of this work. We
will develop the geometrical algorithms required in Sect. 3.2 to evaluate the
inter-time-slab solution. These algorithms are based on the methods presented
in [12]. Specifically, we compute a linear intersection
$\bar{\mathcal{T}}_{\mathrm{int}}^{n}$ of $\mathcal{T}^{n}$,
$\bar{\mathcal{T}}^{n}_{-}=\boldsymbol{\varphi}^{n-1}_{h}(t^{n})(\bar{\mathcal{T}}^{n-1})$
and $\Omega^{n}=\mathrm{int}(\mathcal{B}_{h}(t^{n}))$. It is important to note
that, to ensure the linearity of the intersections, we need to avoid bilinear
terms in $\boldsymbol{\varphi}^{n-1}_{h}$ using a $P_{k}$ space, e.g.,
decomposing $K\in\bar{\mathcal{T}}^{n}$ into simplices.

To simplify the exposition of the intersection algorithm Alg. 1, we redefine
$\mathcal{T}=\bar{\mathcal{T}}^{n}$,
${\mathcal{T}}_{-}=\bar{\mathcal{T}}_{-}^{n}$ and
$\mathcal{B}=\mathcal{B}_{h}(t^{n})$. Within the cell loop of Alg. 1, we first
consider the cells close to $K\in\mathcal{T}$ by restricting $\mathcal{B}$ and
$\mathcal{T}_{-}$ accordingly (line 3 and line 4). These restriction queries
are computed during a preprocessing stage before entering the loop. Next, in
line 5, we compute the intersection of $K$ with the interior of $\mathcal{B}$
using the algorithms in the loop-body of Algorithm 10 in [12].

These algorithms assume that $\mathcal{B}$ is a linear polytope, which is non-
convex in general. Thus, its domain interior $\mathrm{int}(\mathcal{B})$ is
bounded by the set of planar faces of $\mathcal{B}$. The intersection
$K\cap\mathrm{int}(\mathcal{B})$ requires a convex decomposition of
$\mathcal{B}$ and $K$ before intersecting the half-spaces defined by the
planar faces (see [12]).

Finally, we intersect each polytope in $\mathcal{T}_{\mathrm{cut}}^{K}$ by the
subset of $\mathcal{T}_{-}$ around $K$ (line 8). These intersections are
performed employing convex linear clipping algorithms [37] described in
Algorithm 2 in [12]. Alternatively, in line 10, we intersect the cells
$K\in\mathcal{T}$ within the domain $\mathrm{int}(\mathcal{B})$ bounded by
$\mathcal{B}$ that are not intersected by the domain boundary $\mathcal{B}$.
The information about the cells inside $\mathrm{int}(\mathcal{B})$ is obtained
from the propagation through the untouched cells (see [12]). The returned
triangulation $\mathcal{T}_{\mathrm{int}}$ not only covers the cells cut cells
intersected by the boundary $\mathcal{B}$ but the entire domain enclosed by
$\mathcal{B}$. This process guarantees that each cell
$K_{\mathrm{int}}\in{\mathcal{T}}_{\mathrm{int}}$ has an injective map to
$K\in\mathcal{T}$ and $K_{-}\in\mathcal{T}_{-}$. The cells in
${\mathcal{T}}_{\mathrm{int}}$ are general polytopes that cannot use standard
quadrature rules. In order to numerically integrate in these cells, one can
perform a decomposition of the cells into simplices. Alternatively, one can
reduce the dimension of the integrals with Stokes theorem and moment-fitting
methods; see more details in [13] and references therein.

Algorithm 1 $\mathcal{T}\cap\mathcal{T}_{-}\cap\mathrm{int}(\mathcal{B})$

1:
$\mathcal{T}_{\mathrm{cut}}\leftarrow\emptyset,\quad\mathcal{T}_{\mathrm{in}}\leftarrow\emptyset$

2: for $K\in\mathcal{T}$ do

3: $\mathcal{B}^{K}\leftarrow\mathtt{restrict}(\mathcal{B},K)$

4: $\mathcal{T}_{-}^{K}\leftarrow\mathtt{restrict}(\mathcal{T}_{-},K)$

5: $\mathcal{T}_{\mathrm{cut}}^{K}\leftarrow
K\cap\mathrm{int}(\mathcal{B}^{K})$

6: for $K_{-}\in\mathcal{T}_{-}^{K}$ do

7: if $K\cap\mathcal{B}\neq\emptyset$ then

8: $\mathcal{T}_{\mathrm{cut}}^{K,K_{-}}\leftarrow\\{K_{\mathrm{cut}}\cap
K_{-}:K_{\mathrm{cut}}\in\mathcal{T}_{\mathrm{cut}}^{K}\\};\quad\mathcal{T}_{\mathrm{cut}}\leftarrow\mathcal{T}_{\mathrm{cut}}\cup\mathcal{T}_{\mathrm{cut}}^{K,K_{-}}$

9: else if $K\subset\mathrm{int}(\mathcal{B})$ then

10: $\mathcal{T}_{\mathrm{in}}\leftarrow\mathcal{T}_{\mathrm{in}}\cup(K\cap
K_{-})$

11: end if

12: end for

13: end for

14: return
$\mathcal{T}_{\mathrm{int}}\leftarrow\mathcal{T}_{\mathrm{cut}}\cap\mathcal{T}_{\mathrm{in}}$

We emphasize that the intersection algorithm is as robust as the core cut
algorithm in line 5. Furthermore, it is important to note that, for evaluation
purposes, accurate tolerance management is not required in the intersection
process of line 8. While the algorithm is described for the core cut algorithm
with an exact embedded discretization of explicit domain representation [30],
Alg. 1 is general enough to be used with other unfitted discretizations.

##  5\. Numerical experiments

###  5.1. Objectives

In the experiments of this section, we aim to demonstrate the effectiveness of
the presented methods. In particular, we examine the following aspects within
our formulation:

  * •

We evaluate the $hp$-convergence on both 2D and 3D spatial domains, compared
to the results and analysis in [30].

  * •

We explore numerical stability concerning the cut location and the
approximation degrees.

  * •

We assess the applicability of our formulation to complex 2D and 3D moving
domains derived from STL models.

It is important to note that these experiments focus on comparing our
formulation with the one presented in [30]. Recall that [30] uses space-time
embedded discretizations on implicit geometries determined by level sets, and
the results are presented exclusively in 2D+1D domains (as the geometrical
algorithms for 3D+1D domains are significantly more complex). In contrast, we
design a space-time formulation that works on explicit boundary
representations and only require geometrical intersection algorithms in space
only, e.g., 2D and 3D vs. 3D and 4D in [30]. This is possible by pulling back
the problem into an extruded space-time domain using the formulation in Sect.
3.1.

###  5.2. Environment setup

The numerical experiments have been performed on Gadi and Titani
supercomputers. Gadi is a high-end supercomputer at the NCI (Australia) with
4962 nodes, 3074 of them powered by a 2 x 24 core Intel Xeon Platinum 8274
(Cascade Lake) at 3.2 GHz and 192 GB RAM. Titani is located at UPC (Spain)
with 6 nodes, 5 of them powered by a 2 x 12 core Intel Xenon E5-2650L v3 at
1.8 GHz and 256 GB RAM. The algorithms presented in this work have been
implemented in the Julia programming language [38]. The unfitted FE
computations have been performed using the Julia FE library Gridap.jl [39]
version 0.17.17 and the extension package for unfitted methods
GridapEmbedded.jl version 0.8.1 [40]. STLCutters.jl version 0.1.6 [41] has
been used to compute intersection computations on STL geometries. To mitigate
excessive computational costs, the condition numbers are computed in the
1-norm using cond() Julia function.

###  5.3. Space-time convergence tests

To demonstrate optimal convergence rates we solve a simple PDE with a
manufactured solution out of the FE space. Inspired by [30], we solve the
Poisson equation (3) with the following manufactured solution

$u(x,t)=\sin\left(\frac{\pi\alpha
t}{T}\right)\prod_{i=1}^{D}\sin\left(\frac{\pi x_{i}}{L_{i}}\right),$ (17)

where $\\{L_{1},...,L_{D},T\\}$ are the cartesian dimensions of the space-time
domain and the time parameter is set to $\alpha=0.5$ for the experiments shown
in Fig. 6 and Fig. 7. Furthermore, in the equation (3) we set the diffusion
term as $\mu=1$. The space domain is a $n$-cube with a $n$-cubic hole in the
center. The hole is described by the STL of a cube. The lengths at the domain
sides are $L=3$ while the lengths of the hole sides are $l=1$. The time domain
has size $T=1$. The hole is linearly translated in the $x$ direction. The
displacement map is described as $\boldsymbol{D}(x,t)=(0.2t,0,0)$. In all the
cases, the domain discretization is a regular Cartesian grid, with the same
number of elements $n$ in each direction, both space and time.

For implementation reasons, the two-dimensional examples are computed with a
three-dimensional STL. Thus, we build a pseudo-two-dimensional domain that has
only one cell in the $z$ direction. The space-time domain dimensions are
$\\{L,L,\frac{L}{n},T\\}$ while the number of elements in each direction is
$\\{n,n,1,n\\}$. The number of elements per direction in convergence tests is
$n=2^{i},i=3,...,6$ in two-dimensional runs, while $n=2^{i},i=3,...,5$ for
three-dimensional runs.

We analyze the condition number of the system to be inverted in each time-
slab. Since the system matrix is nonsymmetric, as suggested in [30], a
preconditioner for DG in time [42] can be considered. The effectiveness of
this preconditioner depends on the condition numbers the mass and stiffness
matrices, which are defined as follows,

$\mathbf{M}_{ab}=\int_{\hat{Q}^{n}}\mathcal{E}^{n}(\mathbf{\Phi}^{a})\mathcal{E}^{n}(\mathbf{\Phi}^{b})|J_{\Omega^{n}}|{d}\hat{\boldsymbol{x}}{d}t,\quad\mathbf{A}_{ab}=a_{h}\left(\mathcal{E}^{n}(\mathbf{\Phi}^{a}),\mathcal{E}^{n}(\mathbf{\Phi}^{b})\right).$
(18)

The condition numbers presented in Fig. 6(a) and Fig. 6(b) are computed in the
initial time slab of the two-dimensional convergence experiments. We observe
that the condition number of the mass matrix remains nearly constant, while
the condition number of the stiffness matrix scales with
$\mathcal{O}(h^{-2})$. These observations align with the behavior expected for
AgFEM in space-time domains analyzed in [30].

Now, let us analyze the convergence of the error norms. In Fig. 6, we can
observe the accumulated DG error norm convergence with coercivity constant
$c_{\mu}=1$. We observe that using the AgFEM and constant $h/\tau$ the error
converges with $\mathcal{O}(h^{r})$, where the convergence rate $r=\min(p,q)$.
Fig. 7 shows the $L^{2}(\Omega^{n})$ and $H^{1}(\Omega^{n})$ norms at the
final time $t=T$. The convergence rate of $L^{2}(\Omega^{n})$ norm is
$r=\min(p,q)+1$ while the convergence rate of $H^{1}(\Omega^{n})$ norm is
$r=p$. These results are in agreement with the theoretical and analytical
space-time AgFEM results in [23].

(a) Mass matrix

(b) Stiffness matrix

(c) 2D

(d) 3D

Figure 6. Scaling of the condition numbers of the mass and stiffness matrices
in the initial time slab, (a) and (b). Convergence of the space-time DG norm
error ($e=u-u_{h}$) in two and three-dimensional space domains (c) and (d).
Here, $p$ and $q$ represent the space and time approximation order, resp.

(a) 2D

(b) 3D

(c) 2D

(d) 3D

Figure 7. Convergence of the $L^{2}(\Omega^{n})$ and $H^{1}(\Omega^{n})$ norms
of the error ($e=u-u_{h}$) at the final time $t=T$ in two and three-
dimensional space domains. Here, $p$ and $q$ represent the space and time
approximation order, resp.

###  5.4. Moving domains examples

To demonstrate the applicability of the presented algorithms for the
simulation of fluids around moving boundaries, we solve the flow around two
rotating geometries in 2D and 3D. We solve the Navier-Stokes and Stokes
equations resp. with kinematic viscosity $\nu=10^{-2}$. Regardless of the
dimension, we run similar setups. Both share the same discretization order for
pressure $p_{p}=1$, velocity $p_{u}=2$ and time $q$. We utilize the AgFEM
aggregating all cut cells.

We improve the accuracy of both meshes by clustering the cells around the
geometry. See Fig. 8 and Fig. 9 for 2D and 3D, resp. We map each direction
$i\in 1,...,d$ of the Cartesian meshes as follows,

$\displaystyle\phi^{i}_{M}(\hat{x})=\begin{cases}\hat{\boldsymbol{x}}^{i}_{0}\left(\frac{\hat{x}^{i}}{\hat{x}^{i}_{0}}\right)^{\alpha}&\text{if
}\hat{x}^{i}<\hat{x}^{i}_{0},\\\
1-(1-\hat{x}^{i}_{0})\left(\frac{1-\hat{x}^{i}}{1-\hat{x}^{i}_{0}}\right)^{\alpha}&\text{otherwise.}\end{cases}$
(19)

Here, $\hat{x}^{i}=x^{i}-x^{i}_{0}/L^{i}$ is the reference axis of the
direction $x^{i}$ where $x_{o}^{i}$ is the lower value in the $i$ direction,
and $L^{i}$ is the $i$-length. In (19), we use $\hat{x}^{i}_{0}=0.5$ to
determine the region of element condensation and the exponential factor
$\alpha<1$ for the smoothness of the map.

Figure 8. Background spatial mesh $\bar{\mathcal{T}}^{\mathrm{bg}}$ around a
prismatic gear $\mathcal{B}_{0}$ (geometry id 71711 from _Thingi10k_ [43]).
This mesh is the simplex partition of a mapped Cartesian mesh. The uniform
elements are mapped in each direction with $\phi_{M}$ (19) and an exponential
factor $\alpha=0.5$. The Cartesian mesh has $80\times 40$ elements before the
simplex decomposition. The artificial domain $\Omega^{\mathrm{art}}$ is a box
of size $4.8L_{x}\times 2.4L_{y}$. Here, $L_{x}$ and $L_{y}$ represent the
bounding box size of $\mathcal{B}_{0}$. Figure 9. Background spatial mesh
$\bar{\mathcal{T}}^{\mathrm{bg}}$ around a wing $\mathcal{B}_{0}$ (id 65604 in
_Thingi10k_). This simplex mesh comes from a mapped Cartesian mesh. The
coordinates of the mesh are mapped in the $x$ and $y$ directions with
$\phi_{M}$ (19) and an exponential factor $\alpha=0.3$. The Cartesian mesh has
$20\times 20\times 8$ elements. The artificial domain $\Omega^{\mathrm{art}}$
is a box of size $4L_{x}\times 4L_{y}\times 0.8L_{z}$. Here, $L_{x}$, $L_{z}$
and $L_{y}$ represent the bounding box size of $\mathcal{B}_{0}$.

In the 2D example of Fig. 8 and Fig. 10, we solve the Navier-Stokes equations
with $Re=10^{2}$. We define a parabolic inlet flow in the $x$-direction,

$u(\boldsymbol{x})=U_{\mathrm{max}}(4x_{2}-4x_{2}^{2}),$ (20)

where $U_{\mathrm{max}}=(1,0,0)$. We set zero velocity on the $y$-faces and
zero $z$-velocity in the $z$-faces, i.e., slipping conditions. We weakly
impose the displacement velocity $\boldsymbol{\nabla}\boldsymbol{\varphi}^{n}$
with Nistche’s method on the geometry $\mathcal{B}(t)$. The initial geometry
$\mathcal{B}_{0}$ is mapped by
$\boldsymbol{D}=\boldsymbol{D}_{s}\circ\boldsymbol{\phi}_{t}$,

$\boldsymbol{D}_{s}(\boldsymbol{x},t)=\boldsymbol{x}_{0}+\mathbf{R}(t)\cdot(\boldsymbol{x}-\boldsymbol{x}_{0})+\mathbf{A}_{x}\sin(\omega_{x}t),$

where $\mathbf{R}(t)$ is the rotation matrix with angular velocity
$\omega=\pi/2~{}\mathrm{rad/s}$, $\boldsymbol{x}_{0}$ is the center of the
mesh, $\mathbf{A}_{x}=(0,0.2)$ and $\omega_{x}=\pi/2$. The time map is defined
as, $\boldsymbol{\phi}_{t}(\boldsymbol{x},t)=(\boldsymbol{x},\phi_{t}(t))$,

$\phi_{t}=\begin{cases}\frac{t_{a}}{\gamma}\left(\frac{t}{t_{a}}\right)^{\gamma}&\text{if
}t<t_{a},\\\ (t-t_{a})+\frac{t_{a}}{\gamma}&\text{otherwise.}\end{cases}$ (21)

We set $\gamma=2$ and $t_{a}=1/8$. With this time map, the initial velocity of
the geometry is zero.

$0.0$
---
$2.2$
---
$u_{h}$
---

(a) $t=0.0$

(b) $t=0.5$

(c) $t=1.0$

Figure 10.  Representation of the line integral convolution (LIC) filter of
the fluid simulation around an evolving geometry in Fig. 8 at $Re=10^{2}$. The
time step size is $\tau=1/60$. The geometry combines rotation and
displacements. The red tooth serves as a reference for visualizing the
rotation.

In the 3D experiment of Fig. 9 and Fig. 11, we simulate a viscous flow
($Re\sim 0$) with the same boundary conditions as in the 2D experiment.
However, we define a different inlet profile,

$u(\boldsymbol{x})=U_{\mathrm{max}}(4x_{2}-4x_{2}^{2})(4x_{3}-4x_{3}^{2}),$
(22)

that utilizes the 3D domain. The differences in the $z$-direction are clearly
shown in Fig. 11. Equally, the $z$-faces have a slipping condition to hold the
rotation of the wing in the $z$-direction. The initial geometry
$\mathcal{B}^{0}$ is mapped by

$\boldsymbol{D}(\boldsymbol{x},t)=\boldsymbol{x}_{0}+\mathbf{R}_{\theta}(t)(\boldsymbol{x}-\boldsymbol{x}_{0}),$

where $\boldsymbol{x}_{0}$ is the center of the mesh and $\mathbf{R}_{\theta}$
is the rotation matrix of the angle
$\theta(t)=\theta_{\max}\sin(\omega_{\theta}t)$ over the $z$-axis. Here,
$\theta_{\max}=\pi/10$ and $\omega_{\theta}=\pi$. The 3D experiment is
initialized with a non-zero velocity of the geometry.

We slightly optimize the number of time slab transfer operations in these
experiments. Above a deformation threshold, e.g.,
$\max(\boldsymbol{\nabla}\boldsymbol{\varphi}^{n})<0.8$, we maintain the
spatial active mesh $\bar{\mathcal{T}}^{n}$ between time slabs. In these
cases, the initial value evaluation does not require mesh intersections,
reducing computational cost.

$0.0$
---
$1.0$
---
$u_{h}$
---

(a) $t=0.0,\,\hat{z}=0.5$

(b) $t=0.1,\,\hat{z}=0.5$

(c) $t=0.2,\,\hat{z}=0.5$

(d) $t=0.0,\,\hat{z}=0.75$

(e) $t=0.1,\,\hat{z}=0.75$

(f) $t=0.2,\,\hat{z}=0.75$

(g) $t=0.0,\,\hat{z}=1.0$

(h) $t=0.1,\,\hat{z}=1.0$

(i) $t=0.2,\,\hat{z}=1.0$

Figure 11. Representation of the LIC filter of the viscous flow around Fig. 9.
The matrix representation shows several time slabs on different slices. Here,
$\hat{z}=(z-z_{0})/L_{z}$ represents the relative position in $z$ where
$z_{0}$ is the lower the $z$-coordinate and $L_{z}$ the $z$-length. The step
size of the simulation is $\tau=1/120$. Note that the wing is not an extrusion
of an airfoil. The central section ($\hat{z}=0.5$) has a flat trailing edge
(see Fig. 9).

##  6\. Conclusions and future work

In this work, we have introduced a novel space-time formulation for problems
on moving domains with large displacements. The continuous problem at each
time slab is pulled back to a constant-in-time domain (obtained by extrusion
in time of the initial geometry). The resulting problem is discretized using
unfitted FEM. Intersection algorithms between the boundary of the initial
geometry and background mesh are used to compute the required integrals in the
formulation. This way, we eliminate the requirement for four-dimensional
geometric algorithms. The time discretization relies on DG schemes; the jump
terms across time slabs are exactly integrated by computing the intersection
of the reference and push-forward background meshes.

We have validated this method through an $hp$-convergence analysis of a
manufactured solution using AgFEM. We have observed optimal convergence rates
for two-dimensional and three-dimensional space domains in these experiments.
The convergence has been validated by comparing it with another space-time
analysis for unfitted FEM [30]. Additionally, we have observed the expected
scaling of the condition number of the mass and stiffness matrices.
Furthermore, we have demonstrated the practical applications of this
methodology to the simulation of incompressible flows around rotating
geometries in two examples, one with a two-dimensional geometry and another
with a three-dimensional geometry.

Future work involves the numerical analysis of the method in the case in which
the deformation maps do not much the boundary displacement by designing a
high-order extension of the work in [32] and the extension of this method to
distributed memory machines [44]. The geometrical component of this extension
will be highly scalable, given that the presented algorithms are defined cell-
wise and thus embarrasingly parallel, as they are in [12]. Distributed
computations combined with background mesh refinement, e.g., using octree
meshes, will allow us to solve larger real-world applications. Additional
developments include applying this method to FSI simulations. This extension
will develop the full potential of the method in dynamic interface-coupling
multiphysics simulations.

## Acknowledgments

This research was partially funded by the Australian Government through the
Australian Research Council (project numbers DP210103092 and DP220103160). We
acknowledge Grant PID2021-123611OB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by ERDF “A way of making Europe”. P. A.
Martorell acknowledges the support received from Universitat Politècnica de
Catalunya and Santander Bank through an FPI fellowship (FPI-UPC 2019). This
work was also supported by computational resources provided by the Australian
Government through NCI under the National Computational Merit Allocation
Scheme.

## References

  * Karypis et al. [1997] G. Karypis, K. Schloegel, and V. Kumar.  ParMETIS: Parallel graph partitioning and sparse matrix ordering library.  Technical report, Department of Computer Science and Engineering, University of Minnesota, 1997.
  * Burstedde et al. [2011] 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, Jan. 2011.  doi:10.1137/100791634.
  * Burstedde and Holke [2016] C. Burstedde and J. Holke.  A tetrahedral space-filling curve for nonconforming adaptive meshes.  _SIAM Journal on Scientific Computing_ , 38(5):C471–C503, Jan. 2016.  doi:10.1137/15m1040049.
  * Burman and Fernández [2014] 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, sep 2014.  doi:10.1016/j.cma.2014.07.007.
  * Formaggia et al. [2021] L. Formaggia, F. Gatti, and S. Zonca.  An XFEM/DG approach for fluid-structure interaction problems with contact.  _Applications of Mathematics_ , 66(2):183–211, jan 2021.  doi:10.21136/am.2021.0310-19.
  * Schott et al. [2019] B. Schott, C. Ager, and W. A. Wall.  Monolithic cut finite element–based approaches for fluid-structure interaction.  _International Journal for Numerical Methods in Engineering_ , 119(8):757–796, apr 2019.  doi:10.1002/nme.6072.
  * Dekker et al. [2019] R. Dekker, F. Meer, J. Maljaars, and L. Sluys.  A cohesive XFEM model for simulating fatigue crack growth under mixed-mode loading and overloading.  _International Journal for Numerical Methods in Engineering_ , 118(10):561–577, feb 2019.  doi:10.1002/nme.6026.
  * Giovanardi et al. [2017] B. Giovanardi, L. Formaggia, A. Scotti, and P. Zunino.  Unfitted FEM for modelling the interaction of multiple fractures in a poroelastic medium.  In _Lecture Notes in Computational Science and Engineering_ , pages 331–352. Springer International Publishing, 2017.  doi:10.1007/978-3-319-71431-8_11.
  * Carraturo et al. [2020] M. Carraturo, J. Jomo, S. Kollmannsberger, A. Reali, F. Auricchio, and E. Rank.  Modeling and experimental validation of an immersed thermo-mechanical part-scale analysis for laser powder bed fusion processes.  _Additive Manufacturing_ , 36:101498, dec 2020.  doi:10.1016/j.addma.2020.101498.
  * Neiva et al. [2020] E. Neiva, M. Chiumenti, M. Cervera, E. Salsi, G. Piscopo, S. Badia, A. F. Martín, Z. Chen, C. Lee, and C. Davies.  Numerical modelling of heat transfer and experimental validation in powder-bed fusion with the virtual domain approximation.  _Finite Elements in Analysis and Design_ , 168:103343, jan 2020.  doi:10.1016/j.finel.2019.103343.
  * Badia et al. [2021] S. Badia, J. Hampton, and J. Principe.  Embedded multilevel Monte Carlo for uncertainty quantification in random domains.  _International Journal for Uncertainty Quantification_ , 11(1):119–142, 2021.  doi:10.1615/int.j.uncertaintyquantification.2021032984.
  * Badia et al. [2022] S. Badia, P. A. Martorell, and F. Verdugo.  Geometrical discretisations for unfitted finite elements on explicit boundary representations.  _Journal of Computational Physics_ , 460:111162, jul 2022\.  doi:10.1016/j.jcp.2022.111162.
  * Martorell and Badia [2023] P. A. Martorell and S. Badia.  High order unfitted finite element discretizations for explicit boundary representations.  _arXiv_ , 2023.  doi:10.48550/arXiv.2311.14363.
  * de Prenter et al. [2017] F. de Prenter, C. Verhoosel, G. van Zwieten, and E. van Brummelen.  Condition number analysis and preconditioning of the finite cell method.  _Computer Methods in Applied Mechanics and Engineering_ , 316:297–327, apr 2017.  doi:10.1016/j.cma.2016.07.006.
  * Burman [2010] E. Burman.  Ghost penalty.  _Comptes Rendus Mathematique_ , 348(21-22):1217–1220, nov 2010.  doi:10.1016/j.crma.2010.10.006.
  * Burman et al. [2014] 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, dec 2014.  doi:10.1002/nme.4823.
  * Müller et al. [2016] B. Müller, S. Krämer-Eis, F. Kummer, and M. Oberlack.  A high-order discontinuous galerkin method for compressible flows with immersed boundaries.  _International Journal for Numerical Methods in Engineering_ , 110(1):3–30, nov 2016.  doi:10.1002/nme.5343.
  * Badia et al. [2018a] 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, jul 2018a.  doi:10.1016/j.cma.2018.03.022.
  * Badia et al. [2018b] S. Badia, A. F. Martin, 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, jan 2018b.  doi:10.1137/18m1185624.
  * Verdugo et al. [2019] 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, dec 2019.  doi:10.1016/j.cma.2019.112583.
  * Badia et al. [2021] 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, jan 2021.  doi:10.1137/20m1344512.
  * Neiva and Badia [2021] 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, jul 2021.  doi:10.1016/j.cma.2021.113769.
  * Badia et al. [2022a] S. Badia, E. Neiva, and F. Verdugo.  Robust high-order unfitted finite elements by interpolation-based discrete extension.  _Computers & Mathematics with Applications_, 127:105–126, dec 2022a.  doi:10.1016/j.camwa.2022.09.027.
  * Badia et al. [2022b] S. Badia, E. Neiva, and F. Verdugo.  Linking ghost penalty and aggregated unfitted methods.  _Computer Methods in Applied Mechanics and Engineering_ , 388:114232, jan 2022b.  doi:10.1016/j.cma.2021.114232.
  * Beau et al. [1993] G. L. Beau, S. Ray, S. Aliabadi, and T. Tezduyar.  SUPG finite element computation of compressible flows with the entropy and conservation variables formulations.  _Computer Methods in Applied Mechanics and Engineering_ , 104(3):397–422, may 1993.  doi:10.1016/0045-7825(93)90033-t.
  * Tezduyar et al. [2006] T. E. Tezduyar, S. Sathe, R. Keedy, and K. Stein.  Space–time finite element techniques for computation of fluid–structure interactions.  _Computer Methods in Applied Mechanics and Engineering_ , 195(17-18):2002–2027, mar 2006.  doi:10.1016/j.cma.2004.09.014.
  * Thompson and Pinsky [1996] L. L. Thompson and P. M. Pinsky.  A space-time finite element method for structural acoustics in infinite domains part 1: Formulation, stability and convergence.  _Computer Methods in Applied Mechanics and Engineering_ , 132(3-4):195–227, jun 1996.  doi:10.1016/0045-7825(95)00955-8.
  * Donea et al. [1982] J. Donea, S. Giuliani, and J. Halleux.  An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions.  _Computer Methods in Applied Mechanics and Engineering_ , 33(1-3):689–723, sep 1982.  doi:10.1016/0045-7825(82)90128-1.
  * Nobile and Formaggia [1999] F. Nobile and L. Formaggia.  A stability analysis for the arbitrary lagrangian eulerian formulation with finite elements.  _East-West Journal of Numerical Mathematics_ , 7(2):105–132, 1999.
  * Badia et al. [2023] S. Badia, H. Dilip, and F. Verdugo.  Space-time unfitted finite element methods for time-dependent problems on moving domains.  _Computers & Mathematics with Applications_, 135:60–76, apr 2023.  doi:10.1016/j.camwa.2023.01.032.
  * Heimann et al. [2023] F. Heimann, C. Lehrenfeld, and J. Preuß.  Geometrically higher order unfitted space-time methods for pdes on moving domains.  _SIAM Journal on Scientific Computing_ , 45(2):B139–B165, Mar. 2023.  doi:10.1137/22m1476034.
  * Lehrenfeld and Olshanskii [2019] C. Lehrenfeld and M. Olshanskii.  An eulerian finite element method for pdes in time-dependent domains.  _ESAIM: Mathematical Modelling and Numerical Analysis_ , 53(2):585–614, Mar. 2019.  doi:10.1051/m2an/2018068.
  * de Prenter et al. [2023] F. de Prenter, C. V. Verhoosel, E. H. van Brummelen, M. G. Larson, and S. Badia.  Stability and conditioning of immersed finite element methods: Analysis and remedies.  _Archives of Computational Methods in Engineering_ , 30(6):3617–3656, May 2023.  doi:10.1007/s11831-023-09913-0.
  * Nitsche [1971] 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_ , 36(1):9–15, July 1971.  doi:10.1007/bf02995904.
  * Heimann and Lehrenfeld [2023] F. Heimann and C. Lehrenfeld.  Geometrically higher order unfitted space-time methods for pdes on moving domains: Geometry error analysis.  _arXiv_ , 2023.  doi:10.48550/arXiv.2311.02348.
  * Bonet and Wood [1997] J. Bonet and R. D. Wood.  _Nonlinear continuum mechanics for finite element analysis_.  Cambridge university press, 1997.
  * Sugihara [1994] K. Sugihara.  A robust and consistent algorithm for intersecting convex polyhedra.  _Computer Graphics Forum_ , 13(3):45–54, aug 1994\.  doi:10.1111/1467-8659.1330045.
  * Bezanson et al. [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah.  Julia: A fresh approach to numerical computing.  _SIAM Review_ , 59(1):65–98, jan 2017.  doi:10.1137/141000671.
  * Verdugo and Badia [2022] F. Verdugo and S. Badia.  The software design of gridap: A finite element package based on the julia JIT compiler.  _Computer Physics Communications_ , 276:108341, jul 2022\.  doi:10.1016/j.cpc.2022.108341.
  * Verdugo et al. [2023] F. Verdugo, E. Neiva, and S. Badia.  GridapEmbedded. Version 0.8., Jan. 2023.  Available at https://github.com/gridap/GridapEmbedded.jl.
  * Martorell et al. [2021] P. A. Martorell, S. Badia, and F. Verdugo.  STLCutters.  _Zenodo_ , Sept. 2021.  doi:10.5281/zenodo.5444427.
  * Smears [2016] I. Smears.  Robust and efficient preconditioners for the discontinuous galerkin time-stepping method.  _IMA Journal of Numerical Analysis_ , page drw050, oct 2016.  doi:10.1093/imanum/drw050.
  * Zhou and Jacobson [2016] Q. Zhou and A. Jacobson.  Thingi10K: A Dataset of 10,000 3D-Printing Models.  _arXiv_ , 2016.  doi:10.48550/arXiv.1605.04797.
  * Badia et al. [2022] S. Badia, A. F. Martín, and F. Verdugo.  GridapDistributed: a massively parallel finite element toolbox in julia.  _Journal of Open Source Software_ , 7(74):4157, jun 2022.  doi:10.21105/joss.04157.

  *[FEMs]: finite element method
  *[FSI]: fluid-structure interaction
  *[CAD]: computer-aided design
  *[DG]: discontinuous Galerkin
  *[FE]: finite element
  *[AgFEM]: aggregated finite element method
  *[ALE]: arbitrary Lagrangian Eulerian
  *[FEM]: finite element method
  *[STL]: stereolithography
  *[AMR]: adaptive mesh refinement
  *[DOF]: degree of freedom
  *[DOFs]: degree of freedom
  *[AgFE]: aggregated finite element
  *[CG]: continuous Galerkin
  *[PDEs]: partial differential equation
  *[PDE]: partial differential equation
  *[LIC]: line integral convolution