THE NITSCHE XFEM-DG SPACE-TIME METHOD AND ITS
IMPLEMENTATION IN THREE SPACE DIMENSIONS

arXiv:1408.2941v1 [math.NA] 13 Aug 2014

CHRISTOPH LEHRENFELD

∗

Abstract. In the recent paper [C. Lehrenfeld, A. Reusken, SIAM J. Num. Anal., 51 (2013)] a
new finite element discretization method for a class of two-phase mass transport problems is presented
and analyzed. The transport problem describes mass transport in a domain with an evolving interface.
Across the evolving interface a jump condition has to be satisfies. The discretization in that paper is
a space-time approach which combines a discontinuous Galerkin (DG) technique (in time) with an
extended finite element method (XFEM). Using the Nitsche method the jump condition is enforced in
a weak sense. While the emphasis in that paper was on the analysis and one dimensional numerical
experiments the main contribution of this paper is the discussion of implementation aspects for the
spatially three dimensional case. As the space-time interface is typically given only implicitly as the
zero-level of a level-set function, we construct a piecewise planar approximation of the space-time
interface. This discrete interface is used to divide the space-time domain into its subdomains. An
important component within this decomposition is a new method for dividing four-dimensional prisms
intersected by a piecewise planar space-time interface into simplices. Such a subdivision algorithm is
necessary for numerical integration on the subdomains as well as on the space-time interface. These
numerical integrations are needed in the implementation of the Nitsche XFEM-DG method in three
space dimensions. Corresponding numerical studies are presented and discussed.
Key words. finite elements, evolving surface, parabolic PDE, space-time, two-phase flow,
computations, XFEM, 4D quadrature, discontinuous galerkin
AMS subject classifications. 65D30, 65M30

1. Introduction. We consider two immiscible incompressible fluids which are
contained in a convex polygonal domain Ω ⊂ Rd , d = 2, 3. We assume (for simplicity)
that one fluid is completely surrounded by the other. A typical configuration of this
sort is a droplet moving in a flow field. The domain containing the surrounding fluid is
Ω2 , whereas the inner phase is Ω1 = Ω \ Ω2 . Both domains are separated by a smooth
interface Γ := Ω̄1 ∩ Ω̄2 . The hydrodynamic evolution of the fluids is typically modeled
by the incompressible Navier-Stokes equations. Suitable interface conditions have to
be added to account for continuity and momentum balance which includes surface
tension forces. In [4, 11, 19, 12] one can find a more detailed discussion on the model.
For the mass transport problem under consideration in this paper we assume that the
velocity field w as well as the evolution of the domains Ω1 (t), Ω1 (t) (and the interface
Γ(t)) are known and there holds div w = 0 and VΓ = w · n on Γ(t) with VΓ the velocity
of the interface motion in normal direction and n the unit normal of ∂Ω1 (t). The latter
assumption states that the relative convection velocity across the moving interface
Γ(t) is zero.
We consider the transport of a dissolved species u in such a two-phase flow configuration.
A standard model consists of convection-diffusion equations within both domains,
where Fick’s law is applied as the diffusion model:
∂u
+ w · ∇u − div(α∇u) = f in Ωm (t), m = 1, 2, t ∈ [0, T ],
(1.1a)
∂t
with f a source term and α = α(x, t) the diffusion coefficient which is assumed to be
piecewise constant,
α = αm (t) > 0

in Ωm (t).

∗ Institut für Geometrie und Praktische Mathematik, RWTH Aachen, D-52056 Aachen, Germany;
email: lehrenfeld@igpm.rwth-aachen.de

1

and different in the subdomains α1 =
6 α2 . The solutions in the subdomains are coupled
through interface conditions
[α∇u · n]Γ = 0

at Γ(t),

(1.1b)

[βu]Γ = 0

at Γ(t),

(1.1c)

where [v]Γ denotes the jump across the interface Γ of a (sufficiently smooth) quantity
v, [v]Γ = (v1 )|Γ − (v2 )|Γ , where vm = v|Ωm is the restriction of v to Ωm . The first
condition (1.1b) is derived from the conservation of mass principle. The condition in
(1.1c) is called Henry condition, c.f. [26, 25, 5, 4] and accounts for difference in the
solubilities of the species in the different fluids. The coefficient β = β(x, t) is piecewise
constant:
β = βm (t) > 0

in Ωm (t).

As the solubilities do not coincide, in general we have β1 6= β2 which leads to a jump
discontinuity of u across the interface. We assume that the ratio between β1 and β2 is
bounded by a constant c > 0:
1
β1 ≤ β2 ≤ cβ1 .
c
The fact that the solution has a discontinuity across the (moving) interface is the
demanding aspect of the problem.
To close the model we add initial and boundary conditions:
u(·, 0) = u0 in Ωm (0), m = 1, 2,
u(·, t) = 0

on ∂Ω, t ∈ [0, T ].

(1.1d)
(1.1e)

For the treatment of the problem the case of a non-stationary interface is significantly
different from the case of a stationary interface. In this paper we discuss several
aspects of the finite element discretization method for the non-stationary interface
proposed in [23].
The method in [23] is based on a Nitsche-XFEM method. The Nitsche-XFEM method
has been introduced and studied for different applications in [13, 15, 14, 16, 2, 8].
For such transport problems as in (1.1a)-(1.1e) the Nitsche-XFEM method has been
investigated in [13, 24, 22]. In [13] the Nitsche-XFEM method is analyzed for the
simpler situation of a stationary problem without a discontinuity, i.e. β1 = β2 and
without convection, i.e. w = 0. This method is extended to the instationary case with a
stationary interface, a jump discontinuity across the interface and a convection velocity
w 6= 0 in [24]. The problem of stabilizing the method in [24] w.r.t. a dominating
convection term is discussed in [22]. In all these papers, the Nitsche-XFEM method
is analyzed only for the case of a stationary interface. In [29] a parabolic problem
with a moving interface is discretized with a backward Euler/Nitsche-XFEM method
and analyzed. This paper however only covers the case without discontinuities, i.e.
β1 = β2 and does not consider convection terms.
In [12, Chapter 10] it is shown that a well-posed weak formulation of the parabolic
problem (1.1a)-(1.1e) is most naturally formulated in terms of a space-time variational
formulation. A combination of the Nitsche-XFEM method with a suitable space-time
discontinuous Galerkin method, the space-time XFEM-DG method is presented and
analyzed in [23]. The discretization has a (proven) second order accuracy in space and
time for the transport problem in (1.1a)-(1.1e).
2

An important characteristic of the Nitsche-XFEM method is the fact that the triangulation is not fitted to the interface. This also holds true for the space-time XFEM-DG
method which renders it an Eulerian method.
The idea of combining space-time and extended finite elements has already been
presented in [9] and applied to spatially one dimensional conservation equations
without any analysis.
In implementations of the spatially three dimensional case one has to deal with four
dimensional geometries. In the literature space-time discretizations which explicitly
construct four dimensional triangulations are rare. In [3] a decomposition of four
dimensional prisms into simplices is applied to achieve real space-time adaptivity within
one time slab. More global four dimensional simplex triangulations are considered in
[21, 20].
The main new contributions of this paper are the following two:
1. A strategy to decompose the space-time domain which is intersected by an
implicitly defined space-time interface into simplices and
2. numerical studies for the Nitsche XFEM-DG method in three space dimensions.
As in [23], we restrict to the diffusion dominated case, i.e., no stabilization w.r.t.
convection is needed.
structure and content of this paper. In section 2 we present the discretization
method that is considered in the analysis in [23]. This method is often not feasible
in practice, due to the fact that it is assumed that volume integrals close to the
space-time interface and surface integrals over the space-time interface are evaluated
exactly. In practice a strategy for approximating these integrals will be necessary.
Implementation aspects are discussed in section 3. Especially a strategy to construct
a piecewise planar space-time interface approximation and corresponding polygonal
approximations of the space-time subdomains is presented. A crucial point in this
strategy is the decomposition of prisms (intersected by a piecewise planar interface)
into (uncut) simplices. In the spatially three dimensional case this involves four
dimensional geometries. A main contribution of this paper is a solution algorithm
for this problem that is presented separately in section 5. In section 4 numerical
experiments in three spatial dimensions are presented and discussed.
2. The Nitsche XFEM-DG discretization. In this section we present the
discretization method. We introduce some notation. The space-time cylinder Q =
Ω × (0, T ] ⊂ Rd+1 is partitioned into N time slabs Qn = Ω × In with In := (tn−1 , tn ]
and 0 = t0 < t1 < . . . < tN = T . The length of each time interval In is ∆tn = tn −tn−1 .
For simplicity we assume ∆t = ∆tn = const. The discretization allows to solve for
the discrete solution time slab by time slab which leads to a computational structure
comparable to other (implicit) time-stepping schemes. In the following we will consider
the problem on one time slab Qn , we will therefore omit the subscript n.
Within each time slab we assume that the triangulation of the spatial domain Ω is a
shape regular decomposition into simplices T = {T }. The corresponding characteristic
mesh size is denoted by h. Corresponding to this triangulation of Ω, we have a
canonical triangulation of Q into d + 1-dimensional prisms. This triangulation is
denoted by T ∗ = {QT } where for each prism we have QT = T × In for a corresponding
d-dimensional simplex T . Note, that for different time slabs the triangulation is allowed
to change. Further, the triangulation is not fitted to the interface Γ(t) (cf. figure 2.1).
We introduce the space-time finite element space used for both trial and test functions
in the discretization. Based on Vh , the spatial finite element space of continuous,
piecewise linear functions (with zero boundary values on ∂Ω), the space-time finite
3

t

Γ∗

tn

Qn2
tn−1

Qn−1
1

x

t
tn

Qn2

Γ∗

Qn−1
1

tn−1

x

Fig. 2.1. Sketch of the space-time domains Qn
m . Note that within each time slab the triangulation
has a tensor product structure T × [tn−1 , tn ]. The triangulation is allowed to change between the
time slabs.

element space of continuous piecewise bilinear (linear in time, linear in space) functions
is given by
W := { v : Q → R | v(x, t) = φ0 (x) + tφ1 (x), φ0 , φ1 ∈ Vh }.

(2.1)

S
The space-time interface Γ∗ := t∈In Γ(t) separates the time slab into its subdomains
S
Qm := t∈In Ωm (t), m = 1, 2. With Rm the restriction operator on L2 (Q), such that
Rm v = v|Qm , we introduce the space-time XFEM space as
W Γ∗ := R1 W ⊕ R2 W.

(2.2)

To enforce the Henry condition (1.1c) we use a space-time variant of the Nitsche method.
To do this we need suitable jumps and averages across Γ∗ . Within a prism-element
QT , the averaging operator is defined as
{v}Γ∗ (x, t) := κ1 (v1 )|Γ∗ + κ2 (v2 )|Γ∗ ,

(x, t) ∈ Γ∗

with QTm := Qm ∩ QT and the averaging weights
κm|QT :=

|QTm |
.
|QT |

This is the space-time version of the averaging choice proposed in [13].
We use a similar notation for the jump operators:
[v]Γ∗ (x, t) = (v2 )|Γ∗ − (v1 )|Γ∗ ,

(x, t) ∈ Γ∗ .

Remark 2.1. The choice for κm |QT is slightly different from the one in [23]. As κm |QT
is time-independent within one element QT it is better suited for implementation.
Under mild assumptions on the (approximated) space-time interface the theoretical
results in [23] carry over with only minor adjustments.
Further, for the weak enforcement of continuity in time (using the discontinuous
n−1
Galerkin (DG) method) we use the definition v+
(·) := lim↓0 v(·, tn−1 + ).
4

To transform the iterated integrals which appear in the Nitsche formulation to an
integral on the space-time interface we use the following transformation formula:
Z

tn

Z

Z

f (s) 1 + (w · nΓ )2

f (s, t) ds dt =
tn−1

Γ(t)

− 21

Z
ds =:

Γ∗

f (s)ν(s) ds,
Γ∗

− 1
with ν(s) = 1 + (w · nΓ )2 2 and w the (interface) velocity. Note that both the
surface measure on Γ(t) as well as on Γ∗ are denoted with ds. Under the assumption
that the space-time interface is sufficiently smooth, there holds for a constant c0 > 0
c0 ≤ ν(s) ≤ 1

for all s ∈ Γ∗ .

Below in the Nitsche bilinear form we use a weighting with ν(s).
With the introduced notation we formulate the discrete variational formulation on the
space-time slab Q. We define bilinear forms a(·, ·), b(·, ·) and NΓ∗ (·, ·) corresponding
to the partial differential equation, the weak enforcement of continuity in time (DG)
and the Nitsche terms to weakly impose the interface conditions, respectively. Further
we define the linear forms f (·) and c(·; ·) corresponding to the right hand side source
term f and the initial conditions un−1 :
2 Z
X

∂um
a(u, v) =
+ w · ∇um βm vm + αm βm ∇um · ∇vm dx,
∂t
m=1 Qm
Z
Z
NΓ∗ (u, v) = −
ν(s){α∇u · n}Γ∗ [βv]Γ∗ ds −
ν(s){α∇v · n}Γ∗ [βu]Γ∗ ds
Γ∗
Γ∗
Z
+λh−1
ν(s)[βu]Γ∗ [βv]Γ∗ ds,
T
Γ∗

b(u, v) =

2 Z
X
m=1

f (v) =

2
X
m=1

c(ū; v) =

Ωm (tn−1 )

Z
βm f v dx,
Qm

2 Z
X
m=1

n−1 n−1
βm u+
v+ dx

βm un−1
Ωm (tn−1 )

n−1
v+
dx.

Note that both the volume measure on Ωm (t) as well as on Qm are denoted with dx.
Further, the (space-time) volume integrals are weighted with the factor β. In the
definition of c(·, ·) the notation ū should emphasize that un−1 is given (initial) data
for the computations on the current time slab Q, e.g. for n = 1, un−1 is the initial
data u0 of (1.1d). In NΓ∗ (·, ·) the parameter λ ≥ 0 has to be chosen sufficiently large.
These bi- and linear forms are well-defined on the space-time XFEM space W Γ∗ and
the problem within one time slab reads as: Find u ∈ W Γ∗ , s.t.
B(u, v) := a(u, v) + b(u, v) + N Γ∗ (u, v) = f (v) + c(ū, v)

∀ v ∈ W Γ∗

(2.3)

Note that the solution of (2.3) is the result of one time step within a time stepping
scheme. In [23] it has been shown, that the presented method is a stable discretization
scheme and a second order (in time and space) error bound for the L2 (Ω(T ))-norm
has been proven. We recall the main result from [23].
5

S
[0,T ]
[0,T ]
[0,T ]
Theorem 2.1. Let Qm = t∈[0,T ] Ωm (t), m = 1, 2 and u ∈ H 2 (Q1
∪ Q2 )
denote the solution of (1.1a)-(1.1e). Further, let U be the discrete solution of (2.3)
on every timeslab Qn , n = 1, .., N . Under certain assumptions on the regularity of the
solution of the homogeneous backward problem to (1.1a)-(1.1e) (cf. [23] for precise
statement) the following holds:
k(u − U )(·, T )kL2 (Ω) ≤ c(h2 + ∆t2 )kukH 2 (Q[0,T ] ∪Q[0,T ] ) .
1

2

Remark 2.2 (Linear systems). Linear systems arising from a Nitsche-XFEM discretization may become ill-conditioned depending on the position of the interface.
For the stationary case this has been discussed in [6, 7]. This is also true for the
discretization considered in this paper. However in our experience a remedy to this
problem is a simple diagonal scaling. In the numerical examples in section 4 we solved
the arising linear systems with a Jacobi-preconditioned GMRES method. In all these
examples the iteration counts stayed within reasonable bounds. We do not know
of any literature that rigorously discusses this issue without introducing additional
stabilizations even for the stationary case. This is topic of ongoing research.
Remark 2.3 (Computational efficiency). A concern with the discretization method
might be that the number of degrees of freedom is approximately twice as much as for
a comparable first-order (in time) method. However the fact that we have a (proven)
second order convergence in time justifies, in our opinion, the additional computational
overhead. A crucial point for the discussion of the computational efficiency is the
question of how to solve the arising linear systems efficiently. At this point we do
not have a satisfying answer. In the implementation used for the numerical examples
in section 4 we use standard iterative solvers and preconditioners (see also remark
2.2) which are not adapted to the space-time context of the method. The question of
how to design more sophisticated solvers which exploit the space-time structure of the
problem is the topic of ongoing research. A promising approach has been presented
in [1] for the case without an interface. It is however not clear if these ideas can be
applied in the context of XFEM methods.
3. Implementation aspects. In order to implement the method presented in
section 2 one has to compute the matrices and vectors representing the bi- and linear
forms. In section 3.1 we explain which types of integrals need to be calculated for that.
In practice the space-time interface Γ∗ is typically not given explicitly, but implicitly,
e.g. as the zero-level of a level-set function. In our application the level-set function
is a piecewise quadratic function in space, and typically only given at discrete time
levels tn . In order to apply quadrature on the space-time objects, we want to have an
explicit representation of the space-time interface. As this is practically not feasible,
we construct an appropriate approximation Γ∗,h which has an explicit representation.
Such an approximation is discussed in section 3.2. A different approach is presented
in the recent paper [18] where quadrature rules for implicitly given domains by means
of moment-fitting are derived.
The special tensor product structure of each time slab is reflected in the construction
and representation of finite element shape functions. A short discussion on the shape
functions, especially w.r.t. XFEM can be found in section 3.3.
Once, the (space-time) geometries and finite element shape functions are defined
one needs suitable quadrature rules for (regular) prisms and simplices. As these are
discussed in many standard references (e.g. [27]) for the case d + 1 ≤ 3, but only rarely
in the 3 + 1-dimensional case, section 3.4 and section 3.5 address this issue. In section
6

3.4 tensor product quadrature rules for prisms are applied for elements that are not
intersected. The more involved situation when elements are cut by Γ∗,h make use of
the decomposition of an intersected d + 1-prism into d + 1-simplices. This is discussed
in 3.5.
3.1. Integral types. Every bi- and linear form P
in (2.3) has a natural decomposition into its element contributions, e.g. a(u, v) = QT ∈T ∗ aT (u, v). We consider
the task of computing the element contributions of the (bi-)linear forms a(·, ·), b(·, ·),
NΓ∗ (·, ·), f (·) and c(ū; ·).
As we need to calculate (approximations of) integrals of different kinds, we categorize
these integrals before we discuss their numerical treatment. We distinguish those
integrals in terms of the sets S we are integrating on. The cases are denoted as case
(m,n,o) where m is the dimension of S, n is the co-dimension of S and o ∈ {c, n}
describes if the set S is intersected by the space-time interface Γ∗ (o=c) or not (o=n).
We recall the notation for a prism QT = T ×In . Accordingly we define QTm = QT ∩Qm .
3.1.1. a(·, ·), f (·): d + 1-dimensional measure, co-dimension 0. Integrals
appearing on each element for a(·, ·) are integrals on d + 1-dimensional objects like
Z

Z

tn

Z
f dx , with Tm (t) = T ∩ Ωm (t)

f dx =
QT
m

tn−1

Tm (t)

We distinguish two different situations: The prism QTm is not intersected by the
(approximated) interface, i.e. the prism is completely in one phase and thus the
volume to integrate on is the prism itself. We consider this as case (d+1,0,n) where
numerical integration can exploit the tensor product structure. If on the other hand
the prism QTm is intersected by the (approximated) interface, the geometry Qm
T is much
more difficult to handle. In that case d+1-dimensional quadrature on subsimplices
has to be applied. This is denoted by case (d+1,0,c).
3.1.2. b(·, ·), c(û, ·): d-dimensional measure, co-dimension 0. The integrals
in the element contributions of b(·, ·) and c(ū; ·) have the form
Z
f dx
Tm (tn−1 )

and thus are d-dimensional measures. Also here, we distinguish the case of a one phase
element (i.e. an element which is not intersected), denoted by case (d,0,n) and the
case of an intersected element, case (d,0,c).
3.1.3. N Γ∗ (·, ·): d-dimensional measure, co-dimension 1. For the spacetime integrals stemming from the Nitsche stabilization bilinear form N Γ∗ (·, ·) on each
element we get terms like
Z
ν(s) · ds.
ΓT
∗

where ΓT∗ = Γ∗ ∩ QT . Some terms also depend on the normal direction nΓ . These
integrals only appear on elements that are intersected. The measure is d-dimensional
on the manifold Γ∗ with co-dimension 1. This case is denoted as case (d,1,c).
7

3.2. Approximation of the space-time interface Γ∗ and the space-time
volumes Qm . For the integration of space-time volumes and the space-time interface
a discrete approximation of Γ∗ and Qm which is feasible for numerical integration has
to be found. This is relevant for the integration cases (d+1,0,c) and (d,1,c).
Note that case (d,0,c) is an integral on a d-dimensional simplex T . To deal with
those we follow the strategies (for d ≤ 3) discussed in [12, Section 7.3]. In order to
get approximations Γ∗,h and Qm,h to the space-time interface and volumes (case
(d+1,0,c)) we proceed similarly as in lower (d) dimensions.
We consider the prism QT with a characteristic spatial length h of T (e.g. diameter)
and the time step size ∆t = tn − tn−1 . We apply regular subdivisions in time and
space. Each edge of T is divided into ms parts of equal length and the time interval
is divided into mt parts (see Figure 3.1). We get mt · mds smaller prisms {Qi } with
spatial resolution h/ms and temporal resolution ∆t/mt .

→

mt = 1, ms = 1

→

mt = 2, ms = 1

mt = 2, ms = 2

Fig. 3.1. An intersected prism in d + 1 dimensions, with d = 2. The original prism (left), the
prisms after uniform subdivision in time (middle) and after subdivision in space and time (right).

Each (smaller) prism Qi is subdivided into d + 1 (d+1)-simplices {Pj } (cf. section
5.2). On Pj the level-set function is interpolated as a linear function in space-time (by
simply evaluating the vertex values only). As the level-set function is now represented
as a linear function on each simplex, the according approximation of the zero-level of
the level-set function is piecewise planar.
3.3. Finite elements. We briefly explain the construction of the basis functions
for W Γ∗ . Let J W = {(i, k), i ∈ {1, .., NV }, k = 1, 2} with NV the number of vertices
of the spatial mesh denote the index set corresponding to basis functions qi,k in W .
Based on J W we can define the index subset of basis functions in W which have to be
“enriched”. These are all basis functions which have a non-zero trace on the space-time
interface:

JΓ∗ := { (i, k) ∈ J W | measd Γ∗ ∩ supp(qi,k ) > 0}.
For every index (i, k) ∈ JΓ∗ we add the following basis function
Γ∗
qi,k
(x, t) := qi,k (x, t) · (HΓ∗ (x, t) − HΓ∗ (xi , tn−2+k )),

(j, k) ∈ JΓ∗

(3.1)

where HΓ∗ is the characteristic function of Q2 with HΓ∗ (x, t) = 1 if (x, t) ∈ Q2 and
Γ∗
Γ∗
zero otherwise. It is ensured that qi,1
(xi , tn−1 ) = 0, qi,2
(xi , tn ) = 0 holds for all
vertices xi . Using this construction it is obvious that the XFEM space W Γ∗ can
8

be decomposed into the standard space W and additional basis functions which are
discontinuous across the space-time interface:
 Γ∗
W Γ∗ = W ⊕ span qi,k
| (i, k) ∈ JΓ∗ .
Let the local index set of the local standard basis functions on one prism QT be
J T = {(i, k), i = 1, .., 4, k = 1, 2}. The according shape functions are
qi,k (x, t) := φi (x)ψk (t),

(i, k) ∈ J T

with φi (x) = λi where λj denotes the barycentric coordinate corresponding to vertex
1
1
(tn − t), ψ2 (t) = ∆t
(t − tn−1 ).
j inside of tetrahedron T and ψ1 (t) = ∆t
3.4. Numerical integration on non-intersected (space-time) volumes
(case (d+1,0,n)). Whenever a (tetrahedral) element T is not intersected by the
interface Γ(t) for the whole time interval (tn−1 , tn ] the volume integrals of a(·, ·) and
f (·) act on (the complete) prismatic element QT . Consider for example the diffusion part on QT . One corresponding matrix
entry for I = (i, k), J = (j, l) ∈ J W is
R
T
T
T
GI,J := g (qi,k , qj,l ) with g (u, v) := QT αm βm ∇u∇v dx dt and can be computed
using iterated integrals:
Z

T

tn

Z

g (qi,k , qj,l ) =

αm βm (∇φi )ψk (∇φj )ψl dx dt
tn−1 T
Z tn

=

Z

Z
αm βm ∇φi ∇φj dx = ∆t Mk,l

ψk ψl dt
tn−1

T

αm βm ∇φi ∇φj dx
T

where {Mk,l }k,l=1,2 is the one-dimensional mass matrix

M=

1/3
1/6

1/6
1/3


.

Thus, in this case the quadrature problem is reduced to a standard problem in d
dimensions. For the convection or r.h.s. term the velocity w and the source term f are
not necessarily separable (in terms of t and x). Nevertheless applying interpolation in
time for f and w, e.g.
f∆t = f (tn−1 , x) · ψ1 (t) + f (tn , x) · ψ2 (t) =

2
X

f (tk , x) · ψk (t)

k=1

or tensor product quadrature rules, e.g.
Z

T

tn

Z
βm (∇φi )ψk (w · ∇φj )ψl dx dt

h (qi,k , qj,l ) :=
tn−1 T

≈

Nip
X

Z
βm (∇φi )(w(ti ) · ∇φj )ψl dx dt

ωi ψk (ti )ψl (ti )
T

i=0

where Nip , ti and ωi are the information of the (1D) quadrature rule, reduces the
complexity to a d-dimensional quadrature problem.
9

3.5. Numerical integration on intersected space-time volumes and the
space-time interface (case (d+1,0,c) and case (d,1,c)). If QT is intersected
some simplices Pj within the decomposition QT = {Pj } (see section 3.2) are intersected
by a planar approximation of the interface. Using the simplex and the (hyper-) plane
(k)
one can find a decomposition of Pj into simplices {Pj } which are no longer intersected
S (k)
and form a decomposition of Pj , Pj = k Pj . Furthermore the plane intersecting
one simplex Pj can also be decomposed into uncut d-dimensional simplices. As this
decomposition is neither obvious nor standard for the case d = 3 a solution strategy is
presented in detail in section 5. Thus one can achieve an explicit decomposition of QTm
into uncut (d+1)-dimensional simplices and of Γ∗h into d-dimensional simplices. Once
this decomposition is determined the integration can be applied simplex by simplex.
Quadrature rules of high order for simplices can be found in standard references (see
e.g. [27]) if the dimension of the simplex is d + 1 ≤ 3. For d = 3, i.e. the simplex is
four-dimensional (a pentatope) this is no longer standard. In the literature only a few
integration rules can be found (see eg. [3] and [27]). In appendix A we quote lower
order rules and strategies to generate higher order ones. Further in appendix B we
comment on the computation of the weighting factor ν(s), s ∈ Γ∗ .
3.6. Numerical integration on intersected space volumes. For the treatment of the cases (d,0,n) and (d,0,c) we refer to [12, Section 7.3].
4. Numerical examples. In this section we present results of numerical experiments. Different from the experiments in [23], where only spatially one-dimensional
situations have been considered, we consider spatially three-dimensional cases to assess
the convergence behavior of the method. In this setting we restrict ourselves to
piecewise linear (in time and space) finite element approximations. In all examples we
consider the L2 (Ω(T ))-error for different space and time resolutions. The time interval
[0, T ] is divided into nt time slabs of equal size. Accordingly the time step size is
∆t = nTt . The spatial domain is always a cuboid which is either a cube or divided into
a small number of cubes. The cubes are divided into (ns )d smaller cubes which are
then divided into tetrahedra. The error behaviour is investigated w.r.t. refinements,
i.e. series of ns and nt . For all computations we used a third order rule from [3]
for the numerical integration on the (sub-)pentatopes. In appendix A we review on
integration rules for pentatopes.
4.1. Moving plane, quasi-1D. This case is the three-dimensional counterpart
to the one-dimensional test case in [23]. The domain is the cube Ω = [0, 2]3 . The “inner”
phase is contained in the domain Ω1 (t) = {x ∈ Ω : |x1 − q(x2 , x3 ) − r(t)| ≤ D/2},
where q : [0, 2]2 → [0, 2] is the graph describing the shape of the domain Ω1 and
r : [0, T ] → R the function describing the time-dependent shift of the interface in
x1 -direction. D = 23 is the width of the domain in x1 -direction. The complementary
domain is Ω2 (t) = Ω \ Ω1 (t).
The velocity field is given as w = ( ∂r
∂t (t), 0, 0). As boundary conditions we apply
periodicity in all directions, u(xi = 0) = u(xi = 2), i = 1, 2, 3. This renders the problem
essentially one-dimensional if q(x2 , x3 ) = const.
We prescribe the r.h.s. source term f , such that the solution is given by
u(x, t) = sin(kπt) · U m (x1 − q(x2 , x3 ) − r(t)),

x ∈ Ωm (t),

m = 1, 2

with
U 1 (y) = ay + by 3

and U 2 (y) = sin(πy)
10

(4.1)

q

U 1 /U 2

x1

y

w
Γ

Γ

Fig. 4.1. Sketch of the geometrical setup of the numerical example in section 4.1 (left) , Ω1 is
blue and Ω2 is green. The time-independent part of the solution U 1 /U 2 (see (4.1)) is also sketched
(right).

where a and b are chosen such that the interface conditions hold.
The diffusivities are (α1 , α2 ) = (1, 2) and the Henry weights (β1 , β2 ) = (1.5, 1), resulting
in a ≈ 1.02728 and b ≈ 6.34294. The problem is considered in the time interval [0, T ]
with T = 1.
4.1.1. Planar (in space and time) interface. We choose q(x2 , x3 ) = 1 and
r(t) = 0.25t, hence the space-time interface is planar. Thus the proposed method for
the approximation of the space-time interface is exact for every mt , ms ≥ 1, where mt
and ms are the number of subdivisions in time and each space direction, respectively,
cf. section 3.2. We choose ms = mt = 1.
100

10−1

10−1

10−2

10−2
nx = 8
nx = 16
nx = 32
nx = 64
nx = 128
order 1,2,3

10−3

10−4

1

2

4

nt = 4
nt = 8
nt = 16
nt = 32
nt = 64
order 1,2

10−3

8
nt

16

32

64

8

16

32
nx

64

128

Fig. 4.2. Convergence in L2 (Ω(T ))-norm w.r.t. refinements in time (left) and space (right) for
test case in section 4.1.1.

In Table 4.1 and Figure 4.2 we give the resulting error kuh − ukL2 (Ω(T )) for different
mesh and time step sizes. Similar to the results in [23] we observe a third order
convergence w.r.t. the time step size ∆t and a second order convergence w.r.t. the
mesh size.
Remark 4.1. In [28, Theorem 12.7] for the corresponding DG-FEM method applied
to the standard heat equation an error bound with third order convergence w.r.t. ∆t
has been derived. Note however that the analysis does not carry over for the case of
the Nitsche-XFEM discretization.
4.1.2. Nonlinear moving interface. We consider q(x2 , x3 ) = 78 + 14 x22 (2 − x2 )2
1
and r(t) = 4π
sin(2πt) which leads to a space-time interface which is no longer
11

nt \ ns
1
2
4
8
16
32
64
eocs

8
0.6826
0.2809
0.0624
0.0755
0.0831
0.0844
0.0846

16
0.7005
0.3138
0.0603
0.0165
0.0233
0.0247
0.0249
1.76

32
0.72760
0.32771
0.07139
0.00840
0.00524
0.00646
0.00667
1.90

64
0.73301
0.33260
0.07557
0.01132
0.00103
0.00147
0.00168
1.99

128
0.742317
0.334165
0.076649
0.012328
0.001585
0.000261
0.000403
2.06

eoct
1.15
2.12
2.64
2.96
2.6
-0.63

Table 4.1
kuh − ukL2 (Ω(T )) for different refinements in time and space for the test case in section 4.1.1.
The last column shows the estimated order of convergence w.r.t. time (eoct ) using the finest spatial
resolution (ns = 128), the last row shows the estimated order of convergence w.r.t. space (eocs )
using the finest temporal resolution (nt = 64).

planar. The geometrical approximation of the space-time interface in this paper is
piecewise planar, i.e. the maximum distance between Γ∗ and its approximation Γ∗,h
converges with second order w.r.t. increasing nt · mt , ns · ms . In Table 4.2 the error
kuh − ukL2 (Ω(T )) on a fixed (fine) spatial grid with resolution 64×64×64 for different
numbers of time steps is listed. In order to investigate the impact of the approximation
of Γ∗ we performed the computation with different numbers of subdivisions ms , mt .
The results, shown in Table 4.2, indicate an error bound behaviour of the form
nt
ms = 1, mt = 1
eoct
ms = 1, mt = 2
eoct
ms = 1, mt = 4
eoct
ms = 1, mt = 8
eoct
ms = 4, mt = 8
eoct

1
2.50
2.49
0.520
0.491
0.491

2
2.89
-0.21
0.817
1.61
0.481
0.11
0.412
0.25
0.412
0.25

4
0.547
2.40
0.168
2.28
0.0985
2.29
0.0910
2.18
0.0909
2.18

8
0.137
2.00
0.0374
2.17
0.0167
2.56
0.0143
2.67
0.0142
2.68

16
0.0342
2.00
0.00878
2.09
0.00284
2.56
0.00189
2.92
0.00179
2.99

32
0.00879
1.96
0.00301
1.54
0.00219
0.37
0.00212
-0.17
0.00207
-0.20

64
0.00317
1.49
0.00241
0.32
0.00236
-0.11
0.00236
-0.15

Table 4.2
Error kuh − ukL2 (Ω(T )) for different temporal refinements and quadrature subdivisions on a
regular 64×64×64 tetrahedral mesh for the test case in section 4.1.2.
2

kuh − ukL2 (Ω(T )) ≤ C1 ∆t3 + C2 (∆t/mt ) + C3 (h)
where C1 is independent of the approximation of Γ∗ . C2 is directly related to the
interface approximation errors. If the interface approximation is exact (as in the last
section) C2 is zero. The function C3 (h) describes the spatial error due to the method
and the piecewise linear interface approximation for the numerical integration. It
is thereby the part of the error that can not be reduced by refinements in time. In
this examples C3 (h) ≈ 0.002. Furthermore, we observe that in this example C3 is
essentially independent of ms .
p
For mt sufficiently large, i.e. mt > C2 /(C1 ∆t) and h sufficiently small, the first term
dominates the error, such that one observes a third order in time convergence. This
does not hold if mt is too small. Especially for ms = mt = 1, the error is converging
with (only) second order, due to a dominating interface approximation error.
Remark 4.2. In order to investigate the additional effort within one time step due
to additional XFEM unknowns, we consider the ratio between the maximum number
12

of extended (XFEM) unknowns and standard (space-time) finite element unknowns
within one time slab. In Figure 4.3 a sketch of the corresponding situation is shown. If
the interface is well resolved, the number of unknowns close to the interface increases
by a factor of 2d−1 for one uniform (spatial) refinement whereas the overall number
of unknowns increases with 2d . Thus the ratio decreases linearly with the spatial
resolution. In Table 4.3 the corresponding numbers for this test case are given which
are in agreement with the expected behaviour.
t

t

Γ∗

standard d.o.f.

XFEM d.o.f.

x
t

Γ∗

standard d.o.f.

XFEM d.o.f.

x

Γ∗

standard d.o.f.

XFEM d.o.f.

x

Γ∗

standard d.o.f.

XFEM d.o.f.

x

t

Fig. 4.3. Sketch of different spatial and temporal resolutions for a spatially one-dimensional
problem. Standard and XFEM degrees of freedom are marked. For nt = 2, ns = 4 (upper left) the
highest ratio between XFEM unknowns and standard unknowns within one time slab is 80%, for
nt = 8, ns = 4 (lower left) the ratio is also 80%, for nt = 2, ns = 16 (upper right) the ratio is 59 %
and for nt = 8, ns = 16 (lower right) the ratio is 41%.

Std. unkn.
nt \ ns
1
2
4
8
16
32
64

736
656
656
656
640
640
640

1024
8
(72%)
(64%)
(64%)
(64%)
(63%)
(63%)
(60%)

3648
3008
3008
2880
2624
2496
2496

8192
16
(45%)
(37%)
(37%)
(35%)
(32%)
(30%)
(30%)

19328
14336
14336
13056
11136
10368
9856

65536
32
(30%)
(22%)
(22%)
(20%)
(17%)
(16%)
(15%)

524288
64
119808 (23%)
78848 (15%)
78848 (15%)
67072 (13%)
52992 (10%)
45824 (9%)
41472 (8%)

4194304
128
813056(19%)
479232(11%)
479232(11%)
383488 (9%)
275456 (7%)
198656 (5%)
182784 (4%)

Table 4.3
Number of standard (space-time) unknowns (first row) and maximal number of additional XFEM
unknowns for one time slab for different spatial and temporal resolutions for test case in section
4.1.2. In brackets the ratio between XFEM and standard unknowns is added.

Remark 4.3. To decrease the (space-time) interface approximation error one can
either choose smaller time steps or a larger refinement factor mt for the construction
of Γ∗,h . The computation with a fixed ∆t = ∆t̃ and mt = m̃t > 1 is cheaper than a
computation with ∆t = ∆t̃/m̃t and mt = 1. In Figure 4.4 a sketch of both strategies is
13

shown. For mt > 1 additional effort due to the decomposition strategy and quadrature
within one time step is required. However if the interface is resolved, this is only
required for a small number of elements (cf. Remark 4.2). The number of time steps
and thereby the number of linear systems that have to be solved however is reduced
by a factor of m̃t . Note that the solution of linear systems is typically the most time
consuming part.
t

Γ∗

t

Γ∗,h

Γ∗

Γ∗,h

tn

tn
tn−1
tn−1
standard d.o.f.

XFEM d.o.f.

x

subtriangulation for approximating Γ∗

standard d.o.f.

x

XFEM d.o.f.

subtriangulation for approximating Γ∗

Fig. 4.4. Sketch of geometry approximation and degrees of freedom for the case mt = 1 and
∆t = ∆t̃/4 (left) and for the case mt = 4 and ∆t = ∆t̃ (right). Note that the approximation quality
of the piecewise planar interface is the same.

4.2. Rotational symmetric solution on a moving sphere. In this example
we consider a rotational symmetric solution for a stationary sphere. The sphere
is then translated with a time-dependent velocity. The time interval is [0, T ] with
T = 0.5 and the domain is the cube Ω = [0, 2]3 . One phase is contained in the domain
U 1 /U 2

y

w
y
Γ

Γ

Fig. 4.5. Sketch of geometrical setup (left) for test case in section 4.2, Ω1 is blue and Ω2 is
green. And sketch of the time-independent part of the solution U 1 /U 2 (see (4.2)) (right).

Ω1 (t) = {x ∈ Ω : kx − (p0 + r(t) · e1 )k ≤ R}, where p0 is the center of the initial sphere
and r(t) the motion of the interface in x1 -direction, e1 is the corresponding unit vector.
R = 31 is the radius of the sphere. The complementary domain is Ω2 (t) = Ω \ Ω1 (t).
The velocity field w is given as w = ( ∂r
∂t (t), 0, 0). As boundary conditions we apply
suitable Dirichlet boundary conditions everywhere.
We prescribe these boundary conditions and the r.h.s. source term f , such that the
solution is given by
u(x, t) = sin(kπt) · U m (kx − (p0 + r(t) · e1 )k),

x ∈ Ωm (t),

m = 1, 2

with
U 1 (y) = a + by 2

and
14

U 2 (y) = cos(πy),

(4.2)

where a and b are chosen s.t. the interface conditions hold. The diffusivities are
(α1 , α2 ) = (10, 20) and the Henry weights (β1 , β2 ) = (2, 1) resulting in a ≈ 1.1569
and b ≈ −8.1621. The problem is considered in the time interval [0, T ] with T = 0.5.
1
We choose p0 = (0.5, 1, 1)T and r(t) = 4π
sin(2πt). For the approximation of the
space-time interface we consider ms = mt = 1. We observe an error behaviour which
is of (at least) second order in time and space (O(h2 + ∆t2 )) (see Figure 4.6). For the
finest spatial resolution (ns = 64) we observe an order around 2.5 for the convergence in
time. In contrast to the previous test cases the spatial error dominates the overall error
already for coarse temporal resolutions. We expect that for finer spatial resolutions
and better geometry approximations (mt > 1, ms ≥ 1) one could retain the third order
convergence in time.
nx = 8
nx = 16
nx = 32
nx = 64
order 1
order 2
order 3

10−1

nt \ns
1
2
4
8
16
32
64

10−2

1

2

4

8
nt

16

32

8
0.185
0.180
0.201
0.208
0.209
0.209
0.209

16
0.1929
0.0520
0.0440
0.0509
0.0527
0.0530
0.0530

32
0.2202
0.0607
0.0113
0.0113
0.0131
0.0135
0.0136

64
0.22794
0.06760
0.01409
0.00268
0.00293
0.00332
0.00340

64

Fig. 4.6. Convergence in L2 (Ω(T ))-norm w.r.t. refinements in time and space as a plot and as
tabulated values for test case in section 4.2.

It is also relevant to study the accuracy of the method w.r.t. the interface condition.
Therefore, in Table 4.4 we consider the error
T

Z

1
2

! 12

Z

kν [βuh ]kL2 (Γ∗,h ) =

2

[βuh ] ds dt
0

Γh (t)

under space and time refinement and also observe an O(h2 + ∆t2 ) behaviour.
nt \ ns
2
4
8
16
32
64
eocs

8
0.0495
0.0430
0.0417
0.0414
0.0413
0.0413

16
0.00700
0.00384
0.00253
0.00205
0.00190
0.00186
4.48

32
0.0198
0.00567
0.00164
0.000716
0.000523
0.000477
1.96

64
0.0587
0.0227
0.00517
0.00117
0.000275
0.000131
1.86

eoct
1.37
2.13
2.14
2.10
1.07

Table 4.4
1
Interface error kν 2 [βuh ]kL2 (Γ∗,h ) for different refinements in time and space for the test case
in section 4.2. The last column shows the estimated order of convergence w.r.t. time (eoct ) using
the finest spatial resolution (ns = 64), the last row shows the estimated order of convergence w.r.t.
space (eocs ) using the finest temporal resolution (nt = 64).

4.3. Deforming bubble in a vortex. As a last example we consider a more
complex configuration. We consider an ellipsoidal bubble which is deforming under a
vortex velocity field. The domain is Ω = [0, 2] × [0, 2] × [0, 1] and the velocity field is
15

given as:
w(x, t) = w(x) = q(r(x)) · (x2 − c2 , c1 − x1 , 0)T

(4.3)

where (c1 , c2 , x3 ) with c1 = c2 = 1 describes the rotation axes of the vortex and
r(x) :=

p

(x1 − c1 )2 + (x2 − c2 )2 .

The angular velocity varies with a changing distance to the rotation axes. The
dependency is described by q(r) which is chosen as follows:

2 510 2
− r)3
π  1 +4 2 313 r (0.9
2
q(r) =
·
10 (x − 0.8) (x − 1)2
10 
0

r≤ 0.9
0.9 <r< 1
1 ≤r

(4.4)

The time interval is [0, T ] with T = 20 which corresponds to one full rotation of the
bubble. Due to q(r) 6= const the bubble deforms during that process. In Figure 4.7
the interface at different time levels is shown. Initially we set

t=0

t=5

t = 10

t = 15

t = 20

Fig. 4.7. Interface position at several times for the test case in section 4.3, view from the top.


Ω1 (0) = 32 (x1 − b01 )2 + 62 (x2 − b02 )2 + 32 (x3 − b03 )2 ≤ 1
with (b01 , b02 , b03 ) = (0.5, 1, 0.5), the center of the initial bubble. The complementary
initial domain is Ω2 (0) = Ω \ Ω1 (0). The interface evolves along the characteristics
given by ẋ(t) = w(x, t) which can be calculated explicitely. The diffusivities are
(α1 , α2 ) = (8 · 10−3 , 4 · 10−3 ) and the Henry weights are (β1 , β2 ) = (1, 1.75). As initial
condition we choose u1 (x, 0) = 1 and u2 (x, 0) = 0. Note that u(·, 0) does not fulfill
√ the
interface condition (1.1c) which leads to a parabolic boundary layer of size O( αt).
For this setup we do not know the exact solution. We therefore computed a reference
solution uref on a spatial fine mesh (96 × 96 × 48) with a characteristic mesh size of
h = 1/96 and a temporal resolution of 1024 timesteps for the whole time interval, i.e.
∆tref ≈ 2 · 10−2 . This solution is denoted by uref . We only investigate the temporal
convergence by varying the number of time steps nt . For the approximation of the
space-time interface we consider ms = 2, mt = 1. A contour plot of the solution and
the convergence table is given in Figure 4.8. We observe a convergence order in time
of nearly two. The range in which we observe this order of convergence is limited by
nt < 256. This is due to several effects. At first the reference solution is different from
the exact solution. Secondly due to the dependence of the finite element space and
the approximated space-time interface on the time step size the impact of the spatial
discretization errors can be different for different time step sizes. Hence u − uref does
contain spatial errors which do not vanish for ∆t → 0.
16

nt
32
64
128
256

ku − uref (·, T )kL2 (Ω)
1.22569 · 10−3
3.32611 · 10−4
9.65349 · 10−5
4.08119 · 10−5

eoct
1.88
1.78
1.24

Fig. 4.8. Numerical solution and convergence table of the test case in section 4.3. The contour
plot (left) shows the numerical solution (ns = 96, nt = 1024) of the test case in section 4.3 at T = 20
in the cutting plane z = 0.5. The table shows the convergence in time w.r.t. a reference solution.

5. A strategy to decompose intersected 4-prisms into pentatopes. In
this section we introduce a decomposition strategy that allows for a decomposition of
four dimensional prisms into pentatopes as needed in section 3.5. This approach is
new.
Firstly, we introduce the definitions of relevant four dimensional geometries in section
5.1. The decomposition of a 4-prism into four pentatopes is presented in section
5.2. This is already needed to construct (via interpolation of the level-set function)
the piecewise planar space-time interface in section 3.2. In section 5.4 a strategy is
presented that allows us to decompose a pentatope which is intersected by a hyperplane
(representing an approximation of the space-time interface) into pentatopes which are
not intersected. Figure 5.1 sketches the algorithmic structure of the decomposition
strategy. In this algorithm we need a particular geometrical object, that we call
hypertriangle, which can be decomposed into six pentatopes following the decomposition
rule in section 5.3.
Remark 5.1 (Small angles). The resulting pentatopes/tetrahedra in this decomposition
can have arbitrary small angles. Note that this does not lead to stability problems as
we are using the decomposition only for numerical integration.
5.1. Definition of simple geometries in four dimensions. By ei ∈ Rn we
denote the i-th unit vector with (ei )j = δi,j for i = 1, .., n and e0 := 0.
Definition 5.1 (4-simplex / pentatope). Let xi ∈ R4 for i = 1, .., 5 and di,j := xi −xj .
Iff the vectors di,1 for i = 2, .., 5 are linearly independent, we call the convex hull
P = conv({xi }i=1,..,5 ) the 4-simplex or pentatope.
Remark 5.2 (reference pentatope). Every pentatope P can be represented as an
affine transformation applied to the reference pentatope P̂ = conv({ei }i=0,..4 ). The
transformation has the form
Φ : P̂ → P, (x̂1 , x̂2 , x̂3 , x̂4 ) →

5
X

λ̂i xi ,

i=1

where λ̂i (x̂1 , x̂2 , x̂3 , x̂4 ) are the barycentric coordinates of P̂ w.r.t. vertex ei−1 .
Definition 5.2 (4-prism). Let xi ∈ R4 for i = 1, .., 4 and y ∈ R4 . Iff {xi }i=1,..,4
defines a 3-simplex (tetrahedron) T = conv({xi }i=1,..,4 ) and y is linearly independent
17

one cut
4-prism

uncut
pentatope

(sec. 5.2)

one uncut
pentatope

four
pentatopes
case 1
(sec. 5.4.1)

one
4-prism

(sec. 5.2)

four uncut
pentatopes

one hypertriangle

(sec. 5.3)

six uncut
pentatopes

cut
pentatope
(sec. 5.4)

case 2
(sec. 5.4.2)

cut geometries

uncut geometries

uncut pentatopes

Fig. 5.1. Algorithmic structure of the decomposition strategy proposed in section 5.

of {di,1 }i=2,..,4 , with di,j := xi − xj , the set
Q = conv({xi }i=1,..,4 , {xi + y}i=1,..,4 ) = {x + αy, x ∈ conv({xi }i=1,..,4 ), α ∈ [0, 1]}
is called 4-prism.
Remark 5.3 (reference 4-prism). Every 4-prism can be represented as an affine
linear transformation applied to the reference 4-prism Q̂ = conv({ei }i=0,..,3 }, {ei +
e4 }i=0,..,3 }). The transformation has the form
Φ : Q̂ → Q, (x̂1 , x̂2 , x̂3 , x̂4 ) →

4
X

µ̂i xi + x̂4 y,

i=1

where µ̂i (x̂1 , x̂2 , x̂3 ) are the barycentric coordinates of the reference tetrahedron T̂ =
conv({ei }i=0,..3 ). The next geometry is a little bit more complex. It later occurs as
one part of a pentatope cut by a hyperplane.
Definition 5.3 (hypertriangle). We define the reference hypertriangle as
Ĥ := {(x1 , x2 , x3 , x4 ) ∈ R4+ , x1 + x2 ≤ 1, x3 + x4 ≤ 1}
= conv({x̂i,j }i=1,..,3,j=1,..,3 ) = K̂ × K̂
where K̂ ⊂ R2 denotes the reference triangle K̂ = conv({χ1 , χ2 , χ3 }) ⊂ R2 with
χi = ei−1 ∈ R2 , i = 1, 2, 3 and x̂i,j = (χi , χj ) ∈ R4 , i, j = 1, 2, 3. Now, let xi,j ∈
R4 , i, j = 1, 2, 3. The convex hull H = conv({xi,j }i,j=1,2,3 ) is called a hypertriangle
18

x3

x3

x3

x1
x2

x1

x2

x1

x2

x4

x4

x4

Fig. 5.2. Sketch of reference geometries. Reference 4-prism Q̂ (left), reference pentatope P̂
(center) and reference hypertriangle Ĥ (right). The dotted line in the left picture are parallel to the
x4 -axes and connect the tetrahedra at x4 = 0 and x4 = 1.

iff there exists a transformation
Φ : Ĥ → R4 , (x̂1 , x̂2 , x̂3 , x̂4 ) →

3 X
3
X

ρ̂i (x̂1 , x̂2 )ρ̂j (x̂3 , x̂4 )xi,j ,

i=1 j=1

where ρ̂i (x̂1 , x̂2 ) is the barycentric coordinate of the reference triangle K̂corresponding
to the vertex χi such that there holds Φ(Ĥ) = H.
Remark 5.4 (Sketches). The sketches in Figures 5.2, 5.3 and 5.4 in this section
show two dimensional parallel projections of four dimensional objects. Straight lines
in the sketch represent a line (in four dimensions) between two vertices. Note that
the preimage of a point in the two dimensional sketch of the parallel projection is a
two-dimensional set.
5.2. Decomposition of a 4-prism into four pentatopes. We consider an
arbitrary prism element QT = T × In with a tetrahedral element T and a time
interval In . For each QT there exists a linear transformation Φ mapping from the
reference 4-prism Q̂ to QT which is of the form Φ(x̂, t̂) = (Φx (x̂), Φt (t̂))T with the
time transformation Φt (t̂) = t̂ · tn + (1 − t̂) · tn−1 and the space transformation Φx (x̂)
mapping from the reference tetrahedron T̂ to T .
It is sufficient to consider the decomposition of the reference 4-prism Q̂ into four
pentatopes as applying Φ to each pentatope of this decomposition results in a valid
decomposition of Q into four pentatopes. With xi := ei−1 and yi := ei−1 + e4 for
i = 1, .., 4 for the reference 4-prism there holds Q̂ = conv({xi }i=1,..,4 , {yi }i=1,..,4 ).
We decompose Q̂ into four pentatopes P̂1 , P̂2 , P̂3 , P̂4 , which are defined as follows:
P̂1 := conv({x1 , x2 , x3 , x4 , y4 }), P̂2 := conv({x1 , x2 , x3 , y3 , y4 })
P̂3 := conv({x1 , x2 , y2 , y3 , y4 }), P̂4 := conv({x1 , y1 , y2 , y3 , y4 })
A sketch of those can be found in Figure 5.3. The pentatopes are disjoint (except for
a part with measure zero) and sum up to the reference prism. To see this, we give the
following characterization of the pentatopes P̂i in terms of constrained sets and their
19

pentatope

P̂1 :

P̂2 :

pentatope

x3

pentatope

x3

x2
x1

x1

P̂3 :

x3

x2
x1

x4

x4

P̂4 :

pentatope

x3

x2
x1

x4

x2

x4

Fig. 5.3. Sketch of pentatopes P̂1 , P̂2 , P̂3 , P̂4 which form a valid decomposition of the reference
4-prism Q̂.

partial sums B̂i =

Si

j=1

P̂j :

P̂1 = {x ∈ Q̂, x3 ≥ x4 },
P̂2 = {x ∈ Q̂, x3 ≤ x4 , x3 + x2 ≥ x4 },

B̂2 = {x ∈ Q̂, x3 + x2 ≥ x4 },

P̂3 = {x ∈ Q̂, x1 + x2 + x3 ≥ x4 , x3 + x2 ≤ x4 },

B̂3 = {x ∈ Q̂, x1 + x2 + x3 ≥ x4 },

P̂4 = {x ∈ Q̂, x1 + x2 + x3 ≤ x4 },

B̂4 = {x ∈ Q̂}.

Note further that the measure of all pentatopes are the same, i.e. meas4 (P̂i ) = 1/24.
5.3. Decomposing the reference hypertriangle. Let ui = x̂1,i , vi = x̂2,i ,
w = x̂3,i , i = 1, .., 3 with x̂i,j as in definition 5.3. We decompose Ĥ into six pentatopes
which are defined as follows:
i

D̂u = conv({u1 , u2 , u3 , v2 , w3 }), D̂v = conv({u1 , v1 , v2 , v3 , w3 })
D̂w = conv({u1 , v2 , w1 , w2 , w3 }), D̂1 = conv({u1 , v1 , v2 , w1 , w3 })
D̂2 = conv({u1 , u2 , v2 , w2 , w3 }), D̂3 = conv({u1 , u3 , v2 , v3 , w3 })
Note that there is a simple structure behind this decomposition. We define the
“diagonal triangle” as K̂diag = conv(u1 , v2 , w3 ). To the three vertices of K̂diag we
add the missing vertices (underlined) of one of the following six triangles
K̂u = conv({u1 , u2 , u3 }), K̂v = conv({v1 , v2 , v3 }), K̂w = conv({w1 , w2 , w3 }),
K̂1 = conv({u1 , v1 , w1 }), K̂2 = conv({u2 , v2 , w2 }), K̂3 = conv({u3 , v3 , w3 }).
A sketch of those pentatopes is given in Figure 5.4. Also here, one can show that the
pentatopes are disjoint (except for a part with measure zero), and sum up to Ĥ. To
this end we divide the hypertriangle according to three binary decisions and define
Ĥi,j,k := Ĥ ∩ {(−1)i x2 ≤ (−1)i x4 } ∩ {(−1)j x1 ≤ (−1)j x3 }
∩ {(−1)k (x1 + x2 ) ≤ (−1)k (x3 + x4 )},

i, j, k = 0, 1

Note that Ĥ1,1,0 and Ĥ0,0,1 are sets with measure zero. All other sets can be identified
with a pentatope from the decomposition:
Ĥ1,1,1 = D̂u , Ĥ1,0,0 = D̂v , Ĥ0,1,0 = D̂w , Ĥ0,0,0 = D̂1 , Ĥ0,1,1 = D̂2 , Ĥ1,0,1 = D̂3 .
20

pentatope

D̂u :

pentatope

x2

x1

pentatope

x2

x1

D̂1 :

pentatope

x4

D̂2 :

pentatope

x3

x2

x1

x2

x1

x4

x3

D̂w :

x3

x4

pentatope

D̂v :

x3

x3

D̂3 :

x3

x2

x1

x4

x4

x2

x1

x4

Fig. 5.4. Sketch of pentatopes D̂u , D̂v , D̂w , D̂1 , D̂2 , D̂3 which form a valid decomposition of
the reference hypertriangle Ĥ. The edges of K̂diag are highlighted in red, whereas the triangles
K̂u , K̂v , K̂w , K̂1 , K̂2 , K̂3 are filled with the corresponding color. The triangle corresponding to each
pentatope is highlighted especially.

5.4. Decomposition of a pentatope intersected by the space-time interface. We assume that the space-time interface is approximated in a piecewise planar
fashion, s.t. within each pentatope the space-time interface is a (hyper-)plane. This
plane divides a pentatope into two parts. Note that due to the pentatope being a
convex set each of the two parts will still be convex. We now consider a pentatope
P which is cut by the plane G = {x ∈ R4 : x · nG = c} which represents the local
approximation of the space-time interface. Each vertex v is marked corresponding
to one of the two halfspaces. Vertices with v · nG < c are marked with a plus (+),
all others with a minus (-). Note that this classification includes the cases where the
space-time interface hits vertices (v · nG = c). We thus can only have two non-trivial
situations:
Case 1: One vertex has a sign that is different from all the others or
Case 2: Two vertices have a sign that is different from the other three vertices.
In the following we will consider these cases separately and construct a decomposition
of the parts into pentatopes. Without loss of generality we assume that the vertices in
the smaller group of vertices are those marked with a plus (+).
21

5.4.1. Case 1: Decomposition into one pentatope and one 4-prism.
We consider the case where one vertex of a pentatope, say x5 , is marked with
a plus (+). All other vertices (x1 , x2 , x3 , x4 ) are marked with a minus (-). The
cutting points of the hyperplane G with the edges are b1 := x1 x5 ∩ G, b2 :=
x2 x5 ∩ G, b3 := x3 x5 ∩ G, b4 := x4 x5 ∩ G. The geometry containing the separated
vertex is the pentatope P + := conv({b1 , b2 , b3 , b4 , x5 }) while the remainder is
Q− := conv({x1 , x2 , x3 , x4 , b1 , b2 , b3 , b4 }). Consider the mapping
Φ : Q̂ → Q− , (x̂1 , x̂2 , x̂3 , x̂4 ) →

4
X

µi (x̂1 , x̂2 , x̂3 )(x̂4 bi + (1 − x̂4 )xi )

i=1

with µ̂(x̂1 , x̂2 , x̂3 ) the barycentric coordinates of the reference tetrahedron T̂ . The
decomposition of the reference 4-prism Q̂ into the four pentatopes P̂i , i = 1, .., 4 as
described in section 5.2 can be used as a triangulation of Q̂. Let Φh be the (pentatope-)
piecewise linear interpolation of Φ at the vertices of this triangulation. One can show
that Φh is an isomorphism between Q̂ and Q− and furthermore for each pentatope P̂i
the image Φh (P̂i ) is again a pentatope. Thus the decomposition rule for the reference
4-prism can also be applied here and we get a valid decomposition by taking the four
pentatopes
P1 = Φh (P̂1 ) = conv({x1 , x2 , x3 , x4 , b4 }), P2 = Φh (P̂2 ) = conv({x1 , x2 , x3 , b3 , b4 })
P3 = Φh (P̂3 ) = conv({x1 , x2 , b2 , b3 , b4 }), P4 = Φh (P̂4 ) = conv({x1 , b1 , b2 , b3 , b4 }).
Decomposition of the space-time interface into tetrahedra for case 1. The triangulation of the interface is trivially obtained with the tetrahedron
P ∩ G = I = conv({b1 , b2 , b3 , b4 }).
5.4.2. Case 2: Decomposition into one 4-prism and one hypertriangle.
Let us consider the case where two vertices of a pentatope are marked with a plus (+),
these are (w.l.o.g.) vertices x4 and x5 . All other vertices (x1 , x2 , x3 ) are marked with
a minus (-). The cutting points of the hyperplane G with the edges are c1 := x1 x4 ∩ G,
c2 := x2 x4 ∩ G, c3 := x3 x4 ∩ G, d1 := x1 x5 ∩ G, d2 := x2 x5 ∩ G, d3 := x3 x5 ∩ G. Thus
we have to decompose the two parts H− and Q+ into pentatopes with
H− := conv({x1 , x2 , x3 , c1 , c2 , c3 , d1 , d2 , d3 }),
Q+ := conv({c1 , c2 , c3 , d1 , d2 , d3 , x4 , x5 }).
Let us start with the decomposition of H− . Consider the mapping
Φ : Ĥ → H− , (x̂1 , x̂2 , x̂3 , x̂4 ) →

3 X
3
X

ρi (x̂1 , x̂2 )ρj (x̂3 , x̂4 )qi,j

i=1 j=1

with qi,1 = xi , qi,2 = ci and qi,3 = di where ρi (x̂1 , x̂2 ) are the barycentric coordinates
of the reference triangle K̂ ⊂ R2 . Following section 5.3, we have a triangulation
of Ĥ into pentatopes {D̂i }. One can show that the (pentatope-) piecewise linear
interpolation Φh of Φ is an isomorphism between Ĥ and H− and each image Φh (D̂i )
is again a pentatope. Therefore we can apply the decomposition of the reference
22

hypertriangle Ĥ into pentatopes to get the six pentatopes
Du = Φh (D̂u ) = conv({x1 , x2 , x3 , c2 , d3 }),
Dv = Φh (D̂v ) = conv({x1 , c1 , c2 , c3 , d3 }),
Dw = Φh (D̂w ) = conv({x1 , c2 , d1 , d2 , d3 }),
D1 = Φh (D̂1 ) = conv({x1 , c1 , c2 , d1 , d3 }),
D2 = Φh (D̂2 ) = conv({x1 , x2 , c2 , d2 , d3 }),
D3 = Φh (D̂3 ) = conv({x1 , x3 , c2 , c3 , d3 }.
We now turn over to Q+ . For notational convenience define c4 := x4 and d4 := x5 .
Thus Q+ = conv({c1 , c2 , c3 , c4 , d1 , d2 , d3 , d4 }). Now the structure is similar to the
−
situation for QS
in Case 1 and we can apply the same procedure and get a valid
decomposition Pi = Q+ with
1 2 3
4
5
Q
1 2 3
3
5
P1 = Φ Q
h (P̂1 ) = conv({c , c , c , x , x }), P2 = Φh (P̂2 ) = conv({c , c , c , d , x })
1 2
2
3
5
Q
1
1
2
3
5
P3 = Φ Q
h (P̂3 ) = conv({c , c , d , d , x }), P4 = Φh (P̂4 ) = conv({c , d , d , d , x })

with ΦQ
h the corresponding piecewise linear transformation for the 4-prism.
Decomposition of the space-time interface into tetrahedra for case 2. With similar techniques as done for the four dimensional volume, we can proceed with the
triangulation of the interface which is isomorph to a 3-prism resulting in tetrahedra Ii :
I1 = conv({c1 , c2 , c3 , d3 }), I2 = conv({c1 , c2 , d2 , d3 }), I3 = conv({c1 , d1 , d2 , d3 }).
Conclusion. The method presented and analysed in [23] is applied to a spatially
three dimensional mass transport problem with a moving discontinuity. A strategy for
constructing a piecewise planar approximate space-time interface and corresponding
polygonal space-time domains is presented. For the integration on the polygonal
space-time domains a method which decomposes the domain into a set of simplices
is presented. Furthermore the construction of the tensor product space-time finite
elements and the (non-standard) issue of quadrature rules in four space dimensions is
addressed. Results of numerical experiments reveal at least a second order convergence
in time.
Acknowledgement. The author wants to thank Arnold Reusken and the referees
for their comments that greatly improved the manuscript. Financial support from
the German Science Foundation (DFG) within the Priority Program (SPP) 1506
“Transport Processes at Fluidic Interfaces” is gratefully acknowledged.
Appendix A. Quadrature rule(s) on pentatopes. We shortly review lower
order quadrature rules on pentatopes and discuss how to achieve higher order rules
using tetrahedron rules, Duffy transformation and 1D Gauss-Jacobi integration rules.
Integration rules are given for the reference pentatope P̂.
R
R
A.1. First order rule. There holds P̂ 1 dx̂ = 1/24 and P̂ q(x̂) dx̂ = 1/120 for
q(x̂) ∈ {x̂1 , x̂2 , x̂3 , x̂4 }, s.t. the following rule is obviously exact for all polynomials up
to degree one:
4
X
1
I (f ) = 1/120
f (ei )
i=0

23

A.2. Third order rule. A third order rule, taken from [3], is as follows:
I 3 (f ) = 1/120

5
X

f (si ) with λj (si ) = α for i 6= j and λj (si ) = β for i = j

i=1

where λj is the barycentric coordinate of vertex j in the reference pentatope and the
coefficients are α = 0.118350341907227374 and β = 0.526598632371090503.
A.3. Higher order rules using the Duffy transformation. A more general
approach to derive integration rules for pentatopes is based on the Duffy transformation
3
4
[10].
Let ŷ =
R
R (ŷ1 , ŷ2 , ŷ3 ) ∈ R , and x̂ = (ŷ, t) ∈ R . The problem to compute
f (x̂) dx̂ = P̂ f (ŷ, t) d(ŷ, t) can be transformed using the transformation (ŷ, t) →
P̂
(1/(1 − t)ŷ, t) = (ỹ, t) (see also Figure A.1 for a sketch):
Z

Z

1

Z

1−t

1−t−ŷ1

Z

Z

1−t−ŷ1 −ŷ2

f (ŷ, t) d(ŷ, t) =

f (ŷ, t) dŷ3 dŷ2 dŷ1 dt

P̂

0

0

0

0

[ỹ = 1/(1 − t)ŷ]

Z

1
3

Z

1

1−ỹZ
1 1−ỹ1 −ỹ2

Z

f ((1 − t)ỹ, t) dỹ3 dỹ2 dỹ1 dt

(1 − t)

=
0

Z
=
0

0

0
1

(1 − t)3

Z

0

f˜(ỹ, t) dỹ dt =

T̂

Z

1

(1 − t)3 g̃(t) dt

0

R
with f˜(ỹ, t) = f ((1 − t)ỹ, t) and g̃(t) = T̂ f˜(ỹ, t) dỹ . In this form one can apply a
one-dimensional integration rule of the form
Z

1

(1 − t)3 g̃(t) dt ≈

0

N
X

ωi g̃(ti )

k=0

where ωi and ti are weights and points of the corresponding quadrature rule. In
order to approximate g̃(ti ) at every integration point ti a standard 3D quadrature
rule can be applied. Let’s assume this 3D quadrature rule has order q accuracy. The
highest order for the pentatope rule at lowest costs is achieved if a Gauss-Jacobi rule
(corresponding to the weight (1 − t)3 ) of order q is used for the numerical integration
w.r.t. t. The resulting quadrature rule is positive, but not symmetric. In principle
also the quadrature rule for the tetrahedron can be derived from lower dimensional
quadrature rules applying the idea recursively. This generic procedure generates
quadrature rules which have slightly more points than symmetric Gauss rules. For
that reason, we use symmetric Gauss rules for the tetrahedron.
Appendix B. Computing the weighting factor ν. The weighting factor ν(s) in
the Nitsche XFEM-DG method can be computed using the space-time normal nΓ∗ of
the space-time interface. One can show that there holds
1

ν(s) = (1 + (w · nΓ )2 )− 2 = k(n1 , .., nd )T k,

s ∈ Γ∗

with nΓ∗ = (n1 , .., nd+1 )T the space-time normal at the interface.
As we use a piecewise planar approximation of the space-time interface consisting of
d-simplices in d + 1 dimensions we have to compute a normal to the d-simplex. It is
known that for d = 2 one can use the standard cross-product to compute the normal.
In the next section we quote a generalized cross-product which allows to do the same
if d = 3.
24

t
t̃= t,
1
ỹ= 1−t
y

1

t̃
1

t= t̃,
y=(1 − t̃)ỹ
1 y
0 ≤ t ≤ 1, 0 ≤ y ≤ 1 − t

1 ỹ
0 ≤ t̃ ≤ 1, 0 ≤ ỹ ≤ 1

Fig. A.1. Sketch of the Duffy transformation for d = 1

B.1. Computing normals to tetrahedra in 4 dimensions. In [17] a generalization of the cross-product is given. Given three vectors u1 , u2 , u3 ∈ R4 one can
compute the cross-product v = X(u1 , u2 , u3 ), s.t.
• X(u1 , u2 , u3 ) = 0 iff u1 , u2 , u3 are linear dependent.
• Iff u1 , u2 , u3 are linear independent then for v = X(u1 , u2 , u3 ), there holds:
v ⊥ ui , i = 1, .., 3.
• αX(u1 , u2 , u3 ) = X(αu1 , u2 , u3 ) = X(u1 , αu2 , u3 ) = X(u1 , u2 , αu3 ), α ∈ R
• X(u1 , u2 , u3 ) = sign(π)X(uπ(1) , uπ(2) , uπ(3) ), where π is a permutation, i.e.
changing the order of the arguments switches the sign.
This cross-product can be used to compute normals to tetrahedra. The computation
is given below:
Given u, v, w ∈ R4 . Compute z = X(u, v, w) ∈ R4 as follows:
a1,2 = u1 · v2 − u2 · v1 ,
a1,3 = u1 · v3 − u3 · v1 ,
z1 = w2 · a3,4 − w3 · a2,4 + w4 · a2,3
a1,4 = u1 · v4 − u4 · v1 ,
z2 = −w1 · a3,4 + w3 · a1,4 − w4 · a1,3
a2,3 = u2 · v3 − u3 · v2 ,
z3 = w1 · a2,4 − w2 · a1,4 + w4 · a1,2
a2,4 = u2 · v4 − u4 · v2 ,
z4 = −w1 · a2,3 + w2 · a1,3 − w3 · a1,2
a3,4 = u3 · v4 − u4 · v3 ,
REFERENCES
[1] S. Basting and S. Weller, Efficient preconditioning of variational time discretization methods
for parabolic partial differential equations, Submitted to Numer. Math.. Preprint available.,
(2013).
[2] R. Becker, E. Burman, and P. Hansbo, A Nitsche extended finite element method for
incompressible elasticity with discontinuous modulus of elasticity, Comput. Methods Appl.
Mech. Engrg., 198 (2009), pp. 3352–3360.
[3] M. Behr, Simplex space-time meshes in finite element simulations, Inter. J. Numer. Meth.
Fluids, 57 (2008), pp. 1421–1434.
[4] D. Bothe, M. Koebe, K. Wielage, J. Prüss, and H.-J. Warnecke, Direct numerical
simulation of mass transfer between rising gas bubbles and water, in Bubbly Flows: Analysis,
Modelling and Calculation, M. Sommerfeld, ed., Heat and Mass Transfer, Springer, 2004.
[5] D. Bothe, M. Koebe, K. Wielage, and H.-J. Warnecke, VOF-simulations of mass transfer
from single bubbles and bubble chains rising in aqueous solutions, in Proceedings 2003
ASME joint U.S.-European Fluids Eng. Conf., Honolulu, 2003, ASME. FEDSM2003-45155.
[6] E. Burman and P. Hansbo, Fictitious domain finite element methods using cut elements: Ii. a
stabilized nitsche method, Applied Numerical Mathematics, 62 (2012), pp. 328 – 341. Third
Chilean Workshop on Numerical Analysis of Partial Differential Equations (WONAPDE
2010).
[7] E. Burman and P. Zunino, Numerical approximation of large contrast problems with the
unfitted nitsche method, in Frontiers in Numerical Analysis - Durham 2010, J. Blowey
and M. Jensen, eds., vol. 85 of Lecture Notes in Computational Science and Engineering,
Springer Berlin Heidelberg, 2012, pp. 227–282.
25

[8] A. Chernov and P. Hansbo, An hp-Nitsche’s method for interface problems with nonconforming
unstructured finite element meshes, Lect. Notes Comput. Sci. Eng., 76 (2011), pp. 153–161.
[9] J. Chessa and T. Belytschko, Arbitrary discontinuities in space-time finite elements by level
sets and x-fem, Inter. J. Numer. Meth. Engng., 61 (2004), pp. 2595–2614.
[10] M. G. Duffy, Quadrature over a pyramid or cube of integrands with a singularity at a vertex,
SIAM J. Numer. Anal., 19 (1982), pp. 1260–1262.
[11] S. Groß, V. Reichelt, and A. Reusken, A finite element based level set method for two-phase
incompressible flows, Comp. Visual. Sci., 9 (2006), pp. 239–257.
[12] S. Groß and A. Reusken, Numerical Methods for Two-phase Incompressible Flows, Springer,
Berlin, 2011.
[13] A. Hansbo and P. Hansbo, An unfitted finite element method, based on Nitsche’s method, for
elliptic interface problems, Comput. Methods Appl. Mech. Engrg., 191 (2002), pp. 5537–
5552.
[14]
, A finite element method for the simulation of strong and weak discontinuities in solid
mechanics, Comput. Methods Appl. Mech. Engrg., 193 (2004), pp. 3523–3540.
[15] A. Hansbo, P. Hansbo, and M. Larson, A finite element method on composite grids based on
Nitsche’s method, Math. Model. Numer. Anal., 37 (2003), pp. 495–514.
[16] P. Hansbo, C. Lovadina, I. Perugia, and G. Sangalli, A Lagrange multiplier method for the
finite element solution of elliptic interface problems using non-matching meshes., Numer.
Math., 100 (2005), pp. 91–115.
[17] S. Holasch, Four-space visualization of 4d objects, master’s thesis, Arizona State University,
1991.
[18] B. Müller, F. Kummer, and M. Oberlack, Highly accurate surface and volume integration
on implicit domains by means of moment-fitting, Int. J. Num. Meth. Eng., (2013).
[19] M. Muradoglu and G. Tryggvason, A front-tracking method for computation of interfacial
flows with soluble surfactant, J. Comput. Phys., 227 (2008), pp. 2238–2262.
[20] M. Neumüller, Space-Time Methods, Fast Solvers and Applications, phd thesis, TU Graz,
2013.
[21] M. Neumüller and O. Steinbach, Refinement of flexible spacetime finite element meshes and
discontinuous galerkin methods, Comp. Vis. Sci., 14 (2011), pp. 189–205.
[22] A. Reusken and C. Lehrenfeld, Nitsche-XFEM with streamline diffusion stabilization for a
two-phase mass transport problem, SIAM J. Sci. Comp., 34 (2012), pp. 2740–2759.
[23]
, Analysis of a Nitsche XFEM-DG discretization for a class of two-phase mass transport
problems, SIAM J. Num. Anal., 51 (2013), pp. 958–983.
[24] A. Reusken and T. Nguyen, Nitsche’s method for a transport problem in two-phase incompressible flows, J. Fourier Anal. Appl., 15 (2009), pp. 663–683.
[25] S. Sadhal, P. Ayyaswamy, and J. Chung, Transport Phenomena with Droplets and Bubbles,
Springer, New York, 1997.
[26] J. Slattery, Advanced Transport Phenomena, Cambridge Universtiy Press, Cambridge, 1999.
[27] A. H. Stroud, Approximate calculation of multiple integrals, Math. Comput., 27 (1973),
pp. 437–440.
[28] V. Thomee, Galerkin finite element methods for parabolic problems, Springer, Berlin, 1997.
[29] P. Zunino, Analysis of backward euler/extended finite element discretization of parabolic
problems with moving interfaces, Comput. Methods Appl. Mech. Engrg., 258 (2013),
pp. 152–165.

26

