High order unfitted finite element methods on level set domains using
isoparametric mappings
Christoph Lehrenfeld

arXiv:1509.02762v1 [math.NA] 9 Sep 2015

Institute for Computational and Applied Mathematics, WWU Münster, Münster, Germany

Abstract
We introduce a new class of unfitted finite element methods with high order accurate numerical integration
over curved surfaces and volumes which are only implicitly defined by level set functions. An unfitted
finite element method which is suitable for the case of piecewise planar interfaces is combined with a
parametric mapping of the underlying mesh resulting in an isoparametric unfitted finite element method.
The parametric mapping is constructed in a way such that the quality of the piecewise planar interface
reconstruction is significantly improved allowing for high order accurate computations of (unfitted) domain
and surface integrals. This approach is new. We present the method, discuss implementational aspects and
present numerical examples which demonstrate the quality and potential of this method.
Keywords: numerical integration, level set, unfitted finite element method, isoparametric finite element
method, high order methods, interface problems

1. Introduction
1.1. Motivation
In the recent years unfitted finite element methods have drawn more and more attention. The possibility
to handle complex geometries without the need for complex and time consuming mesh generation is very
appealing. The methodology of unfitted finite element methods, i.e. methods which are able to cope with
interfaces or boundaries which are not aligned with the grid, has been investigated for different problems. For
boundary value problems with non-aligned boundaries methods such as penalty methods [1, 2], the fictitious
domain method [3, 4], the immersed boundary method [5], and other unfitted finite element methods [6]
have been developed. For unfitted interface problems extended finite element methods (XFEM) have been
developed in (among others) [7, 8, 9, 10, 11]. Also partial differential equations on surfaces which are not
meshed have been considered using the trace finite element method (TraceFEM) [12].
In the community of unfitted finite element methods mostly piecewise linear discretizations are considered. One major issue in the design and realization of high order methods is the problem of numerical
integration on domains which are only prescribed implicitly, for instance as the zero level of a scalar function,
the so-called level set function. The use of standard integration rules (which ignore the cut position on cut
elements) is not a good option as the integrand does not provide the necessary smoothness.
The purpose of this paper is the presentation of a new approach which allows for high order numerical
integration on domains prescribed by level set functions. The approach is new, robust and fairly simple
to implement. The method is geometry-based and can be applied to unfitted interface or boundary value
problems as well as to partial differential equations on surfaces.
1.2. Literature
We briefly discuss the state of the art in the literature to put the new approach into context. One
approach that is often used in order to realize numerical integration on implicit domains consists out of two
step: the approximation of the interface with a piecewise planar interface and a tesselation algorithm to
divide the subdomains of a cut element into simple geometries, on which standard quadrature rules can be
Preprint submitted to Comp. Meth. Appl. Mech. Eng.

September 10, 2015

applied. A famous example for quadrilaterals and cubes is the marching cube algorithm [13]. For simplices
a detailed discussion of this approach can be found in (among others) [14, Chapter 5],[15] for triangles and
tetrahedra and in [16, 17] also for 4-prisms and pentatopes (4-simplices). Many simulation codes which
apply unfitted finite element discretizations, e.g. [18, 19, 20, 21, 22] make use of such strategies. In order to
increase the accuracy of this integration one often combines this approach with subdivisions and adaptivity.
Especially on octree-based meshes this can be done very efficiently [23]. However, this tesselation approach
is only second order accurate (w.r.t. the finest subtriangulation) which complicates the realization of high
order methods.
One approach to solve this problem is the tailoring of quadrature points and weights which provide high
order accurate integration rules for implicit domains. Such a construction of points and weights is based
on the idea of fitting integral moments, cf. [24, 25]. Although this results in accurate integration rules it
has two shortcomings. First, the computation of quadrature rules is fairly involved. This aspect is typically
outwayed by the resulting high accuracy. Secondly, the major problem, quadrature weights are in general
not positive. This can lead to stability problems as positivity of mass or stiffness matrices in finite element
formulations can no longer be guaranteed.
For special cases satisfactory answers to the question of high order numerical integration strategies have
been found which allow for an implementation of high order unfitted finite element methods. We mention
a few interesting approaches. For quadrilaterals and hexahedra in [26] a numerical integration algorithm is
presented which can achieve arbitrary high order accuracy and guarantee positivity of integration weights
at the same time. The approach is based on the idea of interpreting the interface as a piecewise graph over
a hyperplanes. In [27] an unfitted boundary value problem is considered. Instead of aiming at a high order
accurate geometrical approximation of the boundary a correction in the imposition of boundary values is
applied in order to recover a high order method. For the discretization of partial differential equations on
surfaces, in [28] a parametric mapping of the interface from a piecewise planar representation to a high
order representation based on quasi-normal fields is applied. Although the approach can not be carried over
straight-forwardly to the case where also the integration on sub-domains is important, that approach has
important similarities to the approach presented in this paper.
In the literature of extended finite element methods (XFEM), there exist other approaches which are
based on a tesselation approach and aim at improving the resulting piecewise planar approximations by
means of a parametric mapping of the sub-trianguation [29, 30]. These approaches are typically technically
involved, especially in more than two dimensions. Moreover, robustness of these approaches is not clear.
The approach presented in this paper is similar to the approaches in [28, 29, 30] in that it is also based
on a piecewise planar interface which is significantly improved using a parametric mapping. The important
difference is, that we consider a parametric mapping of the underlying mesh rather than the sub-triangulation
or the interface. According to the mapping of the mesh the use of isoparametric finite element is natural.
At this point, we would also like to mention the paper [31] where the construction of a mesh deformation,
which is also used in combination with isoparametric finite elements, is specifically designed to align a given
mesh to a given interface position. The goal in that paper is to obtain a mesh that is conforming to a given
interface without changing the mesh topology while keeping the quality of the mesh. Our approach is in a
similar virtue.
1.3. The concept
In this paper we present a strategy for the numerical approximation of domains implicitly defined by an
approximate signed distance function φ, the level set function. The approach is based on the assumption
that a numerical strategy to deal with interfaces stemming from a piecewise linear level set function exactly
is available.
R
In the context of unfitted finite element methods integrals of the form S f dx have to be computed for
S being an interface or a domain which is only implicitly defined through the level set function φ. Consider
the case S = Ωi = Ω ∩ {φ R< 0} for a polygonal domain Ω which has a suitable triangulation T. As an exact
evaluation of the integral Ωi f dx is in general practically not feasible we consider the approximation of Ωi
with the domain
Ωi,h := Ω ∩ {Ih φ < 0}
2

where Ih is the continuous piecewise linear interpolation of φ. To Ωi,h an explicit representation can be
found on which quadrature rules can be applied element by element:
Z
Z
XX
f dx ≈
f dx ≈
ωi f (xi )
(1.1)
Ωi

Ωi,h

T ∈T

i

Here ωi und xi are the integration weights and points which depend on the cut configuration on T . The
accuracy of this approach is limited by the approximation quality of Ωi,h which is only second order accurate.
By a transformation of the mesh which is represented by an explicit mapping Ψh , parametrized by a
finite element function, the piecewise planar interface is mapped approximately onto the zero level set
of a corresponding high order accurate level set function, cf. Figure 1.1.

Ψh

Figure 1.1: Main idea: The piecewise planar interface on the initial mesh is mapped close to the desired interface via a mesh
transformation Ψh .

According to this mapping we have that Ψh (Ωi,h ) is a high order accurate approximation to Ωi which,
by construction, still has an explicit representation. The integral can hence be approximated as
Z
Z
XX
f dx ≈
f dx ≈
ωi |det(Ψh (xi ))| f (Ψh (xi )),
(1.2)
Ωi

Ψh (Ωi,h )

T ∈T

i

where the integration weights and points are the same as in (1.1). In contrast to (1.1) the accuracy of the
quadrature in (1.2) is no longer bounded to second order accuracy but essentially depends on Ψh . The finite
element spaces (used in a discretization of an (unfitted) PDE problem) have to be adapted correspondingly
which renders a resulting method an isoparametric finite element method: Let Wh be the finite element
space corresponding to the piecewise planar interface approximation with Γ1 := {Ih φ = 0}. Then, after
transformation with Ψh , the appropriate isoparametric finite element space is
Whk := {ϕh ◦ Ψ−1
h |ϕh ∈ Wh }.
Note that this implies that shape functions are not necessarily polynomials on the mapped domain.
From a computational point of view only two things change compared to the case of a piecewise planar
interface. First, a suitable (approximated) mapping Ψh has to be constructed. A major contribution of this
paper is the discussion of this. Secondly, the deformation has to be considered in the definition of the shape
functions and the quadrature of deformed elements. The treatment of the latter aspect is well-known and
is not discussed here.
Owing to the combination of a tesselation algorithm for a piecewise planar interface and a parametric
mapping, we obtain the four most important features of this approach:
• An explicit high order accurate geometrical approximation of the exact interface.
• Guaranteed positiveness of quadrature weights on the interface and in the sub-domains.
3

• Unchanged cut topology compared to the piecewise planar case.
• Easy integration into existing codes as the approach builds on a piecewise planar interface.
Together, this renders the method highly accurate, robust, efficient and fairly simple to implement. We
again note, that the improved geometry approximation is based on a parametric mapping of the initial
mesh, s.t. a combination with isoparametric finite elements is crucial.
1.4. Content and structure of the paper
The main purpose of this study is the presentation of the new approach which allows for an explicite
high order geometry approximation of domains which are implicitly defined as level sets. We focus on the
discussion of an efficient construction of a suitable mesh deformation and related implementational aspects.
Although this discussion goes beyond simple concepts in that it addresses important details of the method,
we do not aim at a thorough error analysis, yet. This is topic of ongoing research and will be published in
a different paper.
The paper is structured in the following way. The starting point for the method is the ability to
deal with interfaces which are piecewise planar. As we mainly consider simplicial meshes this piecewise
planar interface typically corresponds to a (continuous) piecewise linear level set function. We consider a
corresponding decomposition strategy as given and do not discuss it here. Instead, we refer to the literature
instead, cf. [14, Chapter 5] or [17, Chapter 4].
The crucial component of the new approach is the construction of an explicit mapping Ψh which is suitable
for implementation. The construction of Ψh is discussed in section 2 and consists of several steps. We first
characterize the desired properties of the transformation Ψh for a general case. Afterwards we restrict to
an important special case with simplicial elements and a given level set function which is an approximate
signed distance function. For this case we present an explicit construction of a suitable mapping.
Numerical examples demonstrating the quality of the explicit geometrical representation obtained by
this new approach are presented in section 3 and section 4. While the examples in section 3 focus on the
accuracy of the geometry approximation, in section 4 a high order unfitted (isoparametric) finite element
formulation for an interface problem is considered which proves the practical use of the method.
2. Construction of the mesh transformation Ψh
We start with the formulation of properties that we demand from a suitable transformation Ψh . Let
Γ1 be a piecewise planar interface and [Vhk ]d the space of continuous vector-valued piecewise polynomial
functions of degree k. We seek for a transformation Ψh ∈ [Vhk ]d , s.t.
Ψh = argmin dist(Ψh (Γ1 ), Γ)

or at least

Ψh ∈V ∗

dist(Ψh (Γ1 ), Γ) ≤ chk+1

(2.1)

with V ∗ ⊂ [Vhk ]d the subset of admissible transformations which fulfil the following constraints:
• Homeomorphy:

0
d
Ψh , Ψ−1
h ∈ [C (Ω)]

(2.2a)

• Shape regularity: Consider an arbitrary element T in the triangulation T. Let Φh be the transformation
from the reference element T̂ to T = Φh (T̂ ). Then the mapped element T = Ψh (T ) = (Ψh ◦ Φh )(T̂ )
should also be shape regular. Moreover, the shape regularity should be comparable to the shape
regularity of T . We therefore define the transformation
Θh := Ψh ◦ Φh : T̂ → T
and demand
κ(∇Θh ) ≤ C κ(∇Φh ).
for a given constant C > 1.
4

(2.2b)

• Locality:

Ψh (x) 6= x only in the vicinity of the interface Γ1 .

(2.2c)

The last constraint is desirable for efficiency reasons but not crucial.
Note, that due to the fact that Γ1 is already a second order accurate approximation to the interface, the
transformation Ψh is essentially only a local high order correction. This is in contrast to the transformations
tailored in [31]. Further note, that Ψh is obviously not uniquely defined with the requirements formulated
in (2.1) and (2.2a)-(2.2c). Different choices are possible. For the case of a simplicial mesh we introduce a
special choice for Ψh in the following.
In order to construct a suitable transformation Ψh the above given characterization is too general to
be of direct practical use. In the remainder of this study we restrict to simplicial meshes and assume the
interface to be described by an approximate signed distance function φ which is a continuous piecewise
polynomial function. Further, the piecewise planar approximation Γ1 is assumed to be the zero level of the
piecewise linear interpolant Ih φ. These restrictions allow for an explicit construction of a function Ψh which
is of practical use.
We derive this in several steps. In section 2.1 we introduce notation and assumptions. In section 2.2 we
give an ideal mapping Ψ which fulfils Ψ(Γ1 ) = Γk := {φ = 0} (and hence (2.1)) and is defined pointwise.
Together with a suitable modification (section 2.3) and an approximation of this function in [Vhk ]d (section
2.4) we implement conditions (2.1), (2.2a) and (2.2c). In cases where the interface is well resolved by
the piecewise linear level set function Ih φ the transformation Ψh is close to the identity and the shape
regularity of the transformed mesh is ensured. Nevertheless, for practical applications a mechanism to
ensure robustness even if the interface is not well resolved becomes necessary. The aspect of controlling the
shape regularity as in (2.2b) is implemented with a limitation step which is presented and discussed in detail
in section 2.6.
In section 2.7 we summarize the algorithm to determine Ψh and the most important properties of the
method from a computational point of view.
2.1. Notation and assumptions
We introduce some basic notation and assumptions. Ω is a polygonal domain in Rd , d = 2, 3. It is
decomposed into a shape regular partition T of Ω consisting of elements T which are simplices. Inside Ω an
internal interface Γ (Γ ∩ ∂Ω = ∅) separates two domains Ω1 and Ω2 . The interface and the subdomains Γ,
Ω1 , Ω2 are implicitely described by a level set function φex with Γ = {x| φex (x) = 0}, Ω1 = {x| φex (x) < 0}
and Ω2 = {x| φex (x) > 0}. The level set function is a signed distance function:

dist(x, Γ), x ∈ Ω1 ,
ex
φ (x) =
−dist(x, Γ), x ∈ Ω2 .
We make the following assumptions on the smoothness of the interface Γ and the resolution of Γ by the
triangulation T.
Assumption 1. We assume that the interface Γ is C m -smooth, m ≥ 2. Then there exists a δ > 0, such
that φex ∈ H m (Ωδ ) with Ωδ := {x| kx − yk2 ≤ δ for a y ∈ Γ}.
Assumption 2. Let h be the maximum mesh size of the triangulation T. Then we assume that there holds
h < δ with δ as in assumption 1.
Assumption 2 can always be satisfied for sufficiently small mesh sizes. Let Vhk denote the space of continuous
piecewise polynomials of degree k. The level set function φex is approximated by a function φ ∈ Vhk ,
1 ≤ k ≤ m − 1. We think of φ as a projection of φex under a suitable projection operator Πφ : L2 (Ω) → Vhk .
In practice one often only knows φ and not φex . We add an assumptions on the approximation error induced
by the projection Πφ .
Assumption 3. With φ = Πφ φex and Γk := {φ = 0} we assume that there holds
5

kφ − φex kL2 (Ωδ ) ≤ chk+1 |Ωδ |,
ex

k∇φk∞,Ωδ ≤ ck∇φ k∞,Ωδ ,

(2.3a)
(2.3c)

kφ − φex k∞,Ωδ ≤ chk+1 ,
k+1

dist(Γk , Γ) ≤ ch

.

(2.3b)
(2.3d)

with constants c only depending on k and Γ.
Note that assumption 3 is justified for Πφ the standard L2 projection on Ω or the projection discussed in
section 2.5.
The standard nodal interpolation of φ into the space of piecewise linears Vh1 is denoted by Ih φ. The zero
level of that level set function Γ1 := {Ih φ = 0} is piecewise planar which allows for an explicit representation
of the interface and the sub-domains. Γ1 is a second order accurate approximation to Γ. Let
TΓ := {T ∈ T| T ∩ Γ1 6= ∅} and ΩΓ := {x ∈ T | T ∈ TΓ }
be the set and region of cut elements, respectively. Due to assumption 2 the interface is resolved in the
sense that ΩΓ ⊂ Ωδ . With the assumed smoothness (m ≥ 2 in assumption 1) there also holds:
kIh φ − φkL2 (ΩΓ ) ≤ ch2 kφex kH 2 (ΩΓ ) .

(2.4)

We further introduce the following notation for the extension of cut elements by its neighbors:
TΓ,+ := {T ∈ T| ∂T ∩ ΩΓ 6= ∅} and ΩΓ,+ := {x ∈ T | T ∈ TΓ,+ }.
2.2. Ideal transformation Ψ
In this section we define an ideal transformation Ψ which is defined pointwise and is not a finite element
function. We assume the ideal case that φ = φex and that the assumptions 1 and 2 are valid. This specifically
implies that ∇φ is continuous in the set of cut elements ΩΓ and Γ = Γk . At the end of the section we discuss
the case φ 6= φex .
We ask for Ψ(Γ1 ) = Γk = Γ. This only determines the transformation Ψ on Γ1 . Hence, we use a condition
which describes a reasonable extension on ΩΓ . Later, in section 2.4, the continuous transition to Ψh (x) = x
away from the interface is implemented within the approximation in a finite element space.
Definition of the mapping Ψ:. We choose a transformation Ψ which maps iso levels of the piecewise linear
approximation of the level set function {Ih φ = c} onto the corresponding iso levels of the level set function
{φ = c}, i.e. we ask for
Ih φ = φ ◦ Ψ ∀ x ∈ ΩΓ .
(2.5)
Such a transformation is sketched in Figure 2.1. Note that (2.5) implies Ψ(Γ1 ) = Γ. The mapping can also

Ψ

Figure 2.1: Ideal transformation: Piecewise linear level sets are mapped on exact level sets.

be formulated as the following problem: For every point x ∈ ΩΓ with (linearly interpolated) level set value
c = Ih φ(x), find a point x∗ in Ω such that φ(x∗ ) = c.
Note that (2.5) does not define Ψ uniquely, yet. In order to determine a unique transformation, for every
point x in ΩΓ we search for the corresponding point x∗ ∈ {φ(x∗ ) = c} with smallest minimal distance to x.
Due to x ∈ Ωδ and assumption 1 this point is unique. With the help of the signed distance function φ we
6

can characterize x∗ = x + rs using the search direction s = ∇φ and the distance r. By specifying the search
direction s(x) for each point in ΩΓ , we can define Ψ as
Ψ(x) := x + r(x) s(x), x ∈ ΩΓ with r(x) ∈ R such that

Ih φ(x) = φ(x + r(x) s(x)).

(2.6a)
(2.6b)

With φ = φex , we have that s = ∇φex is continuous in ΩΓ which results in Ψ ∈ [C 0 (ΩΓ )]d . By construction
we have Ψ(xV ) = xV for every vertex xV in the mesh as there holds Ih φ(xV ) = φ(xV ).
The non-ideal case φ 6= φex :. In practice one mostly has φ 6= φex . Then, typically s = ∇φ is not continuous
any more and as a result the same holds for a correspondingly defined Ψ. There are essentially two ways
to deal with this problem. Either continuity of the search directions is restored by means of a projection
Πs of ∇φ into [Vhk ]d . Possible projections are for instance an L2 (ΩΓ )-projection or the projection discussed
in section 2.5. Alternatively, one could stick with the discontinous search direction s = ∇φ accepting the
resulting discontinuous Ψ. This is acceptable as finally, for a realization of the transformation, Ψ is projected
into a continuous finite element space with a projector ΠΨ similar or identical to Πs . In our experience both
approaches give similar results. Nevertheless, in the remainder of the paper we consider the use s = Πs ∇φ
as it simplifies the discussion of the method. We also mention that it is not necessary to have a high order
accurate approximation of ∇φ for the search direction. In the numerical examples in section 3 and section
4 a continuous piecewise linear search direction (using the projection discussed in section 2.5) is used.
Another important aspect of the non-ideal case φ 6= φex is the fact, that the transformation Ψ is typically
not smooth (only Lipschitz continuous) within each element T as φ typically has discontinuities in all
derivatives across element interfaces. A typical situation leading to a non-smooth Ψ is sketched in Figure
2.2. In section 2.3 we propose a modification which (alongside its other benefits) results in a transformation
which is smooth within all elements.

T1

T2

Figure 2.2: Sketch of a situation where the ideal transformation Ψ is not smooth within a given element (T2 ). Points on the
piecewise planar level set (green) are mapped on corresponding points on the corresponding level set of φ (blue). Ψ inherits
the missing regularity of φ with a shift in locations: kinks across element interfaces in φ result in kinks of Ψ within elements.

T1

T2

Evaluation of Ψ:. At each point where the transformation Ψ has to be evaluated, the solution of problem
(2.6b) for a fixed x has to be implemented. Let T be the simplex, such that x ∈ T . The search direction s(x)
and the value (Ih φ)(x) can easily be evaluated at the given point. This information is available elementlocal. In order to find the corresponding projected point x∗ , it seems natural to try the following Newton
iteration:
(Ih φ)(x) − φ(xk )
x0 = x, xk+1 = xk −
· s(x)
(2.7)
∇φ(xk ) · s(x)

This procedure has two disadvantages. First, φ may have kinks at element interfaces (as it is only piecewise
polynomial) which may complicate the solution of the non-linear problem. This problem can be tackled
by replacing the Newton-iteration by more sophisticated methods. More importantly, for x0 ∈ T it can
/ T , such that evaluations of φ need to be carried out in a neighborhood of element T .
happen that xk ∈
The evaluation of finite element functions on neighbor elements is often troublesome, esspecially if one is
concerned with an efficient parallelization. We propose a modification which solves both problems in the
next section.
7

2.3. Localized (piecewise smooth) transformation ΨT
We modify the previously discussed construction of the deformation in order to render an implementation
less complicated and more efficient. We achieve this with a localization of the computations on elements.
That means that we avoid evaluating finite element functions from neighbor elements. To this end we first
introduce a discontinuous (across element interfaces) transformation ΨT in section 2.3 which is afterwards
projected into the space of continuous finite element functions with a projection such as the one discussed
in section 2.5.
We define a modified transformation ΨT by its element-wise contributions ΨT and only consider elements
T ∈ TΓ . On each element T the level set φ is a polynomial of order k. Let ET be the extension operator so
that ET (φ) denotes the polynomial on Rd the restriction of which on T coincides again with φ. As φ is an
accurate approximation of a smooth function, ET (φ) is a good approximation to φ also in the neighborhood
of T . Replacing φ with ET (φ) in (2.6b), results in the definition of ΨT :
ΨT (x) := x + rT (x) s(x), x ∈ T with rT (x) ∈ R such that

Ih φ(x) = (ET (φ))(x + rT (x) s(x)).

(2.8a)
(2.8b)

As the function to evaluate in (2.8b) is now a polynomial on Rd the Newton iteration
xk+1 = xk −

(Ih φ)(x) − (ET (φ))(xk )
· s(x)
∇(ET (φ))(xk ) · s(x)

(2.9)

is much simpler to evaluate and should converge even faster. For moderate order k even exact root finding
algorithms could be applied. Another advantage of this approach is that we can guarantee that the resulting
transformation is smooth within each element, cf. Figure 2.3.

T0

T

Figure 2.3: Sketch of the same configuration as in Figure 2.2, this time for the localized transformation ΨT . Points on the
piecewise planar level set (green) are mapped on corresponding points on the corresponding locally extended (dotted) level set
of φ (blue) resulting in a smooth function ΨT within each element T .

We briefly explain why we expect this to be a valid approximation of the original problem. To this end,
we assume that the interface is resolved. In that case ET (φ) only has to be evaluated in a small neighborhood
of T , denoted as ω(T ). In that neighborhood a simple triangle inequality gives

T0

T

kET (φ) − φk∞,ω(T ) ≤ kET (φ) − φex k∞,ω(T ) + kφ − φex k∞,ω(T )

(2.10)

where the latter term can be bounded in O(hk+1 ) due to assumption 3 and the first term on the right hand
side can also be bounded in O(hk+1 ) with a Bramble-Hilbert argument.
ΨT is discontinuous across element boundaries. However, the jumps over element boundaries should be
in the order of O(hk+1 ). Again by a projection ΠΨ we make sure that the final discrete transformation Ψh
is a continuous finite element function. A simple and efficient choice for ΠΨ is discussed in section 2.5.
2.4. Approximation with finite elements
Until now we have discussed how to define and evaluate Ψ in ΩΓ . Our goal, however, is the construction
of a finite element function Ψh ∈ [Vhk ]d which fulfils condition (2.1), but furthermore guarantees (2.2a) and
(2.2c).
8

We write Ψh ∈ [Vhk ]d as Ψh = x + dh with dh ∈ [Vhk ]d the function describing the deformation of the
mesh. Now we can split
NΓ
NΩ
X
X
Γ
dh = dΓh + dΩ
=
u
ϕ
(x)
+
vj ϕ Ω
i
h
i
i (x)
i=1

j=1

with ui ∈ Rd the degrees of freedom corresponding to the basis functions ϕΓi which are located at the
Γ
interface, supp(ϕΓi )∩ΩΓ 6= ∅, and vj the remaining degrees of freedom (supp(ϕΩ
j )∩Ω = ∅). By construction
Γ
Ω
Ω
Γ
we then have Ψh |ΩΓ = (x + dh )|ΩΓ and dh |ΩΓ = 0, i.e. dh has no contribution in Ω .
This splitting decomposes the problem again into two parts. First, dΓh has to be chosen such that Ψh is a
Γ,+
good approximation to Ψ in ΩΓ . Secondly, dΩ
. The second
h has to be chosen such that Ψh (x) = x in Ω \ Ω
Ω
problem immediately suggests to choose dh = 0. Hence, (2.2a) and (2.2c) are fulfilled and we only have to
consider the approximation in (2.1). On that account we apply a projection ΠΨ of ΨT (or Ψ) into [Vhk ]d .
This projection could be a nodal interpolation (only if Ψ is used), an L2 (ΩΓ )-projection or the projection
in section 2.5 and we expect the following assumption to be valid.
Assumption 4. Consider Ψ∗ ∈ {Ψ, ΨT }. With Ψh = ΠΨ Ψ∗ there holds
kΨh − Ψ∗ kL2 (ΩΓ ) ≤ chk+1 kΨ∗ kH k+1 (TΓ ) ≤ chk+1 |ΩΓ |

(2.11)

h→0

where the last inequality holds due to Ψ∗ (x) −→ x. Further there holds
kΨh − Ψ∗ k∞,ΩΓ ≤ chk+1 .

(2.12)

A rigorous justification of this assumption (for both, Ψ and ΨT ) requires further analysis which is topic of
ongoing research.
Consider Ψh = ΠΨ Ψ. Defining Γh := Ψh (Γ1 ) we expect the following error estimate to hold:
dist(Γh , Γ) ≤ kφex k∞,Γh = kφex ◦ Ψh k∞,Γ1 = kφex ◦ Ψh − Ih φk∞,Γ1

≤ kφex ◦Ψh −φ ◦Ψh k∞,Γ1 + kφ ◦Ψ−φ ◦Ψh k∞,Γ1 +kφ ◦Ψ − Ih φk∞,Γ1
{z
}
|
≤ kφex − φk∞,ΩΓ + k∇φk∞,ΩΓ kΨ − Ψh k∞,Γ1 ≤ chk+1 =0
{z
}
|
{z
}
| {z } |
≤ck∇φex k∞,ΩΓ

≤chk+1

(2.13)

≤chk+1

where in the last step we used assumptions 3 and 4. For the case Ψh = ΠΨ ΨT the estimate still applies
(replace Ψ with ΨT ) except for the term kφ ◦ΨT − Ih φk∞,Γ1 in (2.13) which is not exactly zero but can be
estimated as follows:
kφ ◦ΨT − Ih φk∞,Γ1 ≤ max

T ∈TΓ



kET (φ) ◦ΨT − Ih φk∞,Γ1 ∩T +kET (φ) ◦ΨT − φ ◦ΨT k∞,T
|
{z
}
=0

(2.10)

≤ max kET (φ) − φk∞,ω(T ) ≤ chk+1
T ∈TΓ

2.5. An Oswald-type projection
At several occasions throughout this paper we mentioned a generic projection operator Π : [L2 (Ω∗ )]n →
k n
[Vh ] with Ω∗ = Ω or Ω∗ = ΩΓ and n ∈ {1, d}, for instance Πφ , Πs and ΠΨ . An obvious choice for all these
projection is a simple L2 (Ω∗ ) projection. However, to avoid the solution of global linear systems it is worth
mentioning a simple and more efficient alternative which we also used in the numerical examples. For ease
of presentation we only consider the case where we seek a projection Π : L2 (ΩΓ ) → Vhk , i.e. the scalar case
on the interface region. The translation to Ω and/or vector-valued functions is then obvious.
9

The Oswald-type interpolation operator consists of two steps: A projection into a discontinuous finite
element space and an averaging into Vhk . Let Vhk,disc := {u|T ∈ P k (T )| T ∈ TΓ } be the space of piecewise
polynomials which are discontinuous across element boundaries. Due to the missing continuity-restriction,
the standard L2 -projection Πdisc : L2 (Ω) → Vhk,disc can be computed in an element-by-element fashion and
can be implemented very efficiently. Next, we define a simple averaging operator Πav : Vhk,disc → Vhk which
is a high order version of the Oswald interpolation operator [32]. Let ωi := {T ∈ TΓ | supp(ϕi ) ∩ T 6= ∅}
be the set of elements where the ith basis function ϕi of Vhk is supported in ΩΓ . Further, let ui,T be the
coefficients such that
X X
ui,T ϕi |T , u ∈ Vhk,disc .
u=
i

av

We then define Π (u) :=

P

i

T ∈ωi

wi ϕi

with wi =

1 X
ui,T if supp(ϕi ) ∩ ΩΓ 6= ∅ and wi = 0 otherwise.
#ωi
T ∈ωi

Together, we obtain the projection operator Πo : L2 (Ω) → Vhk with Πo = Πav ◦ Πdisc . We note that in an
implementation of Πo the finite element space Vhk,disc does not need to be constructed explicitly.
2.6. Shape regularity
So far we assumed that the interface is smooth and well resolved. An approach which aims to be of
practical use should however be capable of dealing with non-smooth or underresolved interfaces. In this
section we characterize a simple sufficient condition to ensure shape regularity. Based on this we discuss two
different cases. First, we assume the nice case, where the interface is well-resolved and smooth and show
that shape regularity is not an issue. Secondly, we discuss the case where the interface is not resolved and
propose a simple strategy to ensure shape regularity.
Before we discuss the details of both cases, we introduce transformations between the reference simplex,
the simplex before and after the transformation with Ψh and an intermediate curved simplex.
Simplex transformations:. A planar simplex can be characterized as a reference simplex under an affine
transformation. The mesh obtained after the transformation with Ψh can be interpreted as a collection
of simplices under the concatenation of two transformations: an affine one and another one parametrized
by a finite element function. This characterization facilitates the investigation of the shape regularity of
the transformed elements. We therefore introduce notation for the (curved and planar) simplices and the
transformations relating them. In figure 2.4 the simplices and transformations are also depicted.
By T̂ and x̂ we denote the reference simplex and its coordinates. The affine mapping Φh transforms the
reference simplex to the planar simplex T = Φh (T̂ ).
x = Φ(x̂) = A · x̂ + x0

(2.14a)

The final (curved) element T is obtained after the finite element transformation of T by Ψh . The shape
functions (corresponding to Ψh ) are typically defined w.r.t. to the reference simplex T̂ , s.t. we have
y = Ψh (x) = x +

X
i

di (ϕi ◦ Φ−1
h )(x)

x=Φh (x̂)

=

x0 + A · x̂ +

X
i

di · ϕi (x̂).

(2.14b)

Thus, we have T = Θh (T̂ ) with Θh = Ψh ◦ Φh . We can also characterize Θh as a mapping which first
applies the non-planar deformation and the affine mapping afterwards. We therefore introduce the mapping
Ψ̂h , s.t.
Θh = (Ψh ◦ Φh ) = (Φh ◦ Ψ̂h )
(2.14c)
with
ŷ = Ψ̂h (x̂) = x̂ +

X
i

dˆi · ϕi (x̂),
10

dˆi = A−1 · di .

(2.14d)

Ψ̂h (x̂)

Φh (ŷ)

Φh (x̂)

T̂ , ŷ

Ψh (x)

T ,y

T̂, x̂
T, x

Figure 2.4: Reference and transformed domains before and after curving and before and after affine transformation to the
physical domain.

Let x̂V and xV be the coordinates of the vertices of T̂ and T . Then Ψ̂h (x̂V ) = x̂V and Ψh (xV ) = xV .
From the previous observations and we can conclude that
X
∇Θ = ∇Φh +
di ∇ϕi = ∇Φh · ∇Ψ̂h , ∇Φh = A.
(2.15)
i

Hence, if the transformation Ψ̂h : T̂ → T̂ (curving of the reference element) and the transformation Φh :
T → T (affine transformation) are well-behaved, the same also follows for Θ : T̂ → T .
As a measure to quantify the shape regularity we consider the relative spectral condition number of the
Jacobian. We have the simple estimate
κ(∇Θ) ≤ κ(∇Φh ) · κ(∇Ψ̂h ).

(2.16)

We may assume that the initial mesh is shape regular such that κ(∇Φh ) ≤ C for reasonable moderate
number C. In order to fulfil (2.2b) it remains to ask for the impact of the non-linear transformation, i.e.
κ(∇Ψ̂h ).
We have ∇Ψ̂h = I + ∇dˆh with dˆh (x̂) the deformation of T̂ w.r.t. the reference simplex T̂ . As dˆh is a
polynomial on the reference triangle, we can bound its gradient
k∇dˆh k∞,T̂ ≤ ck kdˆh k∞,T̂

(2.17)

with kqk∞,T̂ = maxx∈T̂ kq(x)k2 for a vector-valued function q and kQk∞,T̂ = maxx∈T̂ kQ(x)kF for a matrixvalued function Q. The constant ck only depends on the polynomial degree. Due to the choice of norms, it
is sufficient to have (2.17) for a scalar-valued polynomial.
A simple sufficient condition for the shape regularity, formulated as κ(∇Ψ̂h ) ≤ C for a chosen constant
C > 1, can be derived as follows. With elementary calculations we have
k∇Ψ̂h k2

k(∇Ψ̂h )−1 k2

≤ 1 + k∇dˆh k2

≤ (1 − k∇dˆh k2 )−1

≤ 1 + k∇dˆh kF

≤ (1 − k∇dˆh kF )−1
11

≤ 1 + k∇dˆh k∞,T̂

≤ (1 − k∇dˆh k∞,T̂ )−1 .

With (2.17) we can easily conclude that
kdˆh k∞,T̂ ≤

C −1
ck (C + 1)

(2.18)

is a sufficient condition for
κ(∇Ψ̂h ) ≤ C.
The resolved case:. Due to kIh φ − φkL2 (ΩΓ ) ≤ ch2 there holds dist(Γ1 , Γk ) ≤ ch2 and we can conclude that
kdˆh k∞,T̂ ≤ ch. This implies that the deformation gets arbitrary small which implies (2.18) for sufficiently
small h. We can thus conclude, that the shape regularity is easily provided for sufficiently small mesh sizes,
i.e. for a sufficiently fine resolution of the interface by the mesh.
The under-resolved case:. The considered construction of the transformation may lead to undesirable deformation, e.g. self-intersecting geometries and arbitrary small angles if the interface is not sufficiently well
resolve. This may lead to a blow up in κ(∇Ψ̂h ) or even a singular ∇Ψ̂h at some points.
To avoid these situations, we introduce a “barrier step” which garantuees the quality of the resulting
mesh. The idea is simple. We replace the deformation d = Ψ − x (or d = ΨT − x) with a limited deformation
d¯ = min{γh · s/ksk, d}.

(2.19)

On a fixed element T , we then have dˆh = (∇Φh )−1 ΠΨ d¯ and hence
kdˆh k∞ ≤ cT cπ γ

(2.20)

with cπ the L∞ stability constant of the projection ΠΨ and cT = hk(∇Φh )−1 k∞,T a constant depending on
the shape regularity of T and quasi-uniformity of TΓ .
For sufficiently small parameter γ we can then ensure condition (2.18) and get shape regularity independent of the resolution of the interface. The price, however, is a reduction to an essentially only second order
accurate approximation if the limitation step has to be used (d¯ 6= d).
In cases where the interface is resolved the limitation is not needed which leads to d = d. In the
worst case, where the interface is not resolve, one essentially recovers the quality of the piecewise linear
approximation. In both cases the deformed elements are shape regular which results in a robust method.
2.7. Summary of computational aspects
To clarify on the (low) computational complexity of the resulting scheme we briefly summarize the
algorithmic structure for the computation of Ψh . We restrict to the modified version of the algorithm, cf.
d
section 2.3, which is also used in the
Pnumerical examples. To this end we determine the coefficients di ∈ R
of the representation Ψh (x) = x + i di ϕi (x), di ∈ Rd where ϕi are the basis functions of Vhk :
1. Set coefficients and counter (for later averaging) to zero:
di = 0 and ni = 0, i = 1, .., dim(Vhk ).
2. Compute element contributions (Πdisc ΨT ):
Loop over elements T ∈ TΓ :
(a) Determine the Ψh |T ∈ [P k (T )]d by solution of the local L2 (T ) problem:
Z
Z
Ψh |T · v dx =
ΨT · v dx ∀ v ∈ [P k (T )]d
T

T

For each integration point xj this involves:
• The evaluation of ΨT (xj ):
This requires the solution of (2.8b) by the Newton iteration (2.9).
• The limitation step:
If the deformation ΨT (xj )−xj is too large, it is restricted according to (2.19).
12

(b) Add element contribution to global coefficient vectors (for all d.o.f. i of T ):
• di ← di + dTi where dTi are the basis coefficients corresponding to Ψh |T
• ni ← ni + 1
3. Average the element contributions (Πav ):
di ← di /ni if ni > 0

In this overview we neglected the computation of the search direction s = Πs ∇φ. The corresponding
algorithm goes the same lines with the only difference that the evalution of the right hand side functional
in the (element-)local L2 projection is much simpler.
The resulting deformation Ψh has the following properties:
• Only elements in the neighborhood of Γ are deformed. In Ω \ ΩΓ,+ the mesh stays unchanged.
• Where the resolution of the interface is sufficiently high, the interface can be resolved with high order
accuracy.
• In any case shape regularity can be guaranteed.
3. Numerical examples I: Geometry approximation
Preliminaries
In the following we consider complex geometries described by level set functions and their approximation
with the approach proposed before. For the barrier parameter we choose γ = 0.1 (independent of the degree
k). To solve the one-dimensional nonlinear problem in (2.6b) we choose a (relative) tolerance of 10−14 which
resulted in only of 2 − 3 iterations for each point.
For the quasi-normal field s we use a continuous piecewise linear field. In all examples we use the
transformation Ψh = Πo (ΨT ), cf. section 2.3. The implementation is included in an add-on library ngsxfem
to the finite element library NGSolve [33].
We are interested in the geometrical error between the discrete interface Γh = Ψh (Γ1 ) (with an explicit
description) and the ideal interface Γ = {φex (x) = 0} (with an implicit description). Note that the error
consists of two approximations: The approximation of the level set function φ ≈ φex and the approximation
of Γk = {φ = 0} with Ψh (Γ1 ). The error we are interested in is the maximum distance between the discrete
and the exact interface, dist(Γh , Γ).
3.1. A two-dimensional example: The flower-shape
The geometry in this example is inspired by the one in [34]. Due to its shape we call it the “flower-shape”
example. The domain is a square domain Ω = [−1, 1]2 and the interface is given as
Γ := {(x1 , x2 )| x1 = R(θ) cos(θ), x2 = R(θ) sin(θ), θ ∈ [0, 2π]}

(3.1)

with
R(θ) = r0 + 0.1 sin(ωθ), r0 = 0.5, ω = 8.
p
Γ is the zero level of the level set function φ(x) = x21 + x22 − R(θ) with θ(x) = arctan(x1 /x2 ). The inner
domain is Ω1 := {x ∈ Ω, φ(x) < 0}. A sketch of the geometrical configuration is displayed in Figure 3.1 (left)
together with the initial unstructured mesh which consists of 230 triangles. Starting from that mesh, local
refinements around the interface are successively performed 12 times. The final mesh consists of roughly
1.5 million triangles. We consider different polynomial degrees k, where we always use the same polynomial
degree for the approximation of φ and the deformation Ψh .
In this example φ is not a signed distance function, s.t. dist(Γh , Γ) 6= kφk∞,Γh . However, there holds
k∇φk ≥ 1, s.t. we have dist(Γh , Γ) ≤ kφk∞,Γh . In Figure 3.1 the convergence of kφk∞,Γh is depicted which
is an upper bound for the distance.
If the mesh size h in the vicinity of the interface is sufficiently small the geometrical error decreases
with the (optimal) order k + 1. For coarse grids, where the resolution of the interface is not sufficient the
transformation has to be limited in order to garantuee shape regularity. On the corresponding grids we can
not expect to observe a high order convergence which we also observe on the levels 0 − 2.
13

10−2

kφk∞,Γh

10−4
10−6
k
k
k
k
k
k

10−8

10−10
10−12
0

=1
=2
=3
=4
=6
=8
2

4

6

8

10

12

refinements (towards interface)
Figure 3.1: Mesh and convergence behavior of geometrical error for the flower example.

3.2. A three-dimensional example: The gyroid
In this example we consider the cube Ω = [−1, 1]3 with an interface prescribed by the zero level of the
following level set function:
φ(x, y, z) = cos(πx1 ) sin(πx2 ) + cos(πx2 ) sin(πx3 ) + cos(πx3 ) sin(πx1 )

(3.2)

A sketch of the interface and the initial unstructured mesh, consisting of 577 tetrahedra, is depicted in
Figure 3.2. In a neighborhood of the interface Γ there holds k∇φk ≥ 1, s.t. we again have that the
computed quantity kφk∞,Γh is an upper bound for dist(Γh , Γ). We consider seven adaptive refinements
only towards the interface. The final mesh consists of around 40 million tetrahedra. Note that this time
the interface is not contained in Ω such that Ψh (Ω) 6= Ω. However, for the approximation of Γ and the
evaluation of the error this is not important as φ is trivially extended to R3 .
100
10−2

kφk∞,Γh

10−4
10−6
k=1
k=2
k=3
k=4
O(hk+1 )

10−8
10−10
2

3

4

5

6

7

refinements (towards interface)
Figure 3.2: Mesh and convergence behavior of geometrical error for the gyroid example.

The results in Figure 3.2 reveal that the asymptotic behavior of the method is as good as in the previous
two-dimensional example, the order of convergence O(hk+1 ) is clearly visible. However, the pre-asymptotic
14

regime is significantly larger. This is due to the more difficult geometry. We observe that until refinement
level 4 the error drops only with O(h2 ). However, once the interface is resolved by the piecewise linear
interface Γ1 , the limitation of the deformation is no longer necessary and the high order resolution of the
interface kicks in. On the finest level, the accuracy between different polynomial degrees differs dramatically.
An increase of the polynomial degree by one order results in 2 orders of magnitue decrease in the error.
3.3. Summary of numerical examples in this section
These example represent configurations with realistic and challenging complexity. We conclude that the
method is highly accurate if the interface is smooth and well resolved. If the interface is not well resolved,
we at least obtain a robust method with the “standard” second order accuracy.
4. Numerical examples II: A high order unfitted isoparametric finite element method for an
interface problem
4.1. Elliptic interface problem: The disk
Problem description:. We consider the problem from [17, section 2.5.1.4], where an elliptic interface problem
of the following form is considered.
−div (α∇u) = f

in Ωi ,
on Γ,

(4.1b)

[[βu]] = 0

on Γ,

(4.1c)

on ∂Ω.

(4.1d)

[[α∇u · n]] = 0

u = gD

i = 1, 2,

(4.1a)

Here αi , βi , i = 1, 2 are domain-wise constants, (α1 , α2 ) = (2, 1), (β1 , β2 ) = (1, 3/2). We consider a circular
interface Γ with radius R = 0.6 as interface Γ in the square domain [−1, 1]2 . The data gD and f is chosen
such that the solution to (4.1) is

α2 U (r(x)) + β2 , x ∈ Ω1
u(x) =
α1 U (r(x)) + β1 , x ∈ Ω2
q
πr2
with U (r) = cos( 2 ) and r(x) = x21 + x22
2R
The Nitsche-XFEM discretization:. We consider an isoparametric version of the discretization in [17, section
2], where a combination of an extended finite element (XFE) space combined with a Nitsche formulation
for the unfitted interface is applied. The discretization is compactly described next. For a more thorough
discussion of the underlying (planar) discretization we refer to [17, chapter 2].
Let Ri be the restriction operator to domain Ωi,h with Ω1,h := {Ih φ < 0} and Ω2,h := {Ih φ > 0}, such
that (Ri u)|Ωj = δi,j · u|Ωi with δi,j the Kronecker delta symbol. The extended finite element space is then
defined as VhΓ = R1 Vhk ⊕ R2 Vhk . Basis functions from VhΓ are continuous within the sub-domain Ωi,h but
may be discontinuous across Γ1 . According to the transformation Ψh we define the isoparametric finite
Γ
element space VhΓ = {vh ◦ Ψ−1
h |vh ∈ Vh } with correspondingly mapped shape functions.
The finite element space is not conforming w.r.t. the interface condition (4.1c). To incorporate the
interface condition we consider a Nitsche-type discretization:
Find uh ∈ VhΓ such that
a(uh , vh ) + Nh (uh , vh ) = f (vh ) for all vh ∈ VhΓ
(4.2)
with Nh (u, v) := Nhc (u, v) + Nhc (v, u) + Nhs (u, v) and
X Z
a(u, v) :=
αi ∇u∇v dx,
i=1,2

Nhc (u, v)

Z
:=
Γh

Z
f (v) :=

Ψh (Ωi,h )

f v dx
Ω

Nhs (u, v)

{{−α∇u ·n}}[[βv]] ds,
15

Z
:=
Γh

ᾱλk 2
[[βu]][[βv]] ds,
h

with Γh = Ψh (Γ1 ), ᾱ = 12 (α1 +α2 ) and λ a stabilization parameter that only depends on the shape regularity
of the mesh.
The terms a(uh , vh ) and f (vh ) ensure consistency inside the subdomains, Nhc (uh , vh ) ensures consistency
w.r.t. the interface condition (4.1b). Nhc (vh , uh ) render the discrete l.h.s. operator symmetric while keeping
consistency and Ns (uh , vh ) ensures stability in a consistent way.
An important detail in this discretization is the choice of the element-wise defined averaging operators
{{·}}. We define the averaging weights κi , s.t. {{u}} := κ1 u|Ω1 + κ2 u|Ω2 . We consider the “Heaviside” choice
where κi = 1 if |Ti | > 21 |T | and κi = 0 if |Ti | < 12 |T |. Here |Ti | refers to the cut configuration on the
undeformed mesh. This choice in the averaging renders the scheme in (4.2) stable for arbitrary polynomial
degrees k, independent of the cut position of Γ. The interested reader is referred to the error analysis in [17,
section 2.3], specifically remark 2.3.1 therein.
This test case has been considered for k = 1, 2, 3 in [17, section 2.5.1.4] with a piecewise planar approximation of the interface (on a eight times adaptively refined subtriangulation). We now consider k = 1, 2, 3, 4
in VhΓ in combination with a level set function and a transformation Ψh of the same polynomial degree.
Numerical setup:. We consider the mesh displayed in Figure 1.1 as our initial mesh and successively apply
uniform refinements of the mesh and measure the error in the L2 norm, the broken H 1 (semi-) norm as well
as the error in the interface condition (4.1c) in the L2 (Γh )-norm. To solve the arising linear systems we use
a sparse direct solver. We note that the arising linear systems are extremely ill-conditioned. We observed
that in some cases, for k > 2, the linear systems were so ill-conditioned that the initial residual of the linear
system could only be reduced by a factor of 10−6 .
Remark 4.1. Efficient linear solvers for this kind of problem are rarely addressed in the literature, especially
for high order unfitted finite elements. The only robust strategy (without additional stabilization) for linear
systems arising from Nitsche-XFEM discretizations, that we know of, which is based on a rigorous analysis,
is presented in [35]. In that paper however only the case k = 1 is considered. The design and analysis
of linear solvers for linear systems arising from Nitsche-XFEM discretizations remains a challenging task
which requires further attention.
Numerical results:. We note that no barrier step (cf. section 2.6) has been necessary after one refinement.
Hence, we expect to be in the asymptotic regime w.r.t. the geometry error such that kφk∞,Γh ≤ O(hk+1 )
holds. We observe the (optimal) convergence rates for the volume error that are predicted in the error analysis
in [17, section 2.3]. Note that error sources due to geometry approximation and numerical integration have
not been considered in that error analysis. The error in the interface condition [[βu]] = 0 converges with
O(hk+1 ) which is even half an order better than predicted.
ku − uh kH 1 (Ω1 ∪Ω2 )

ku − uh kL2 (Ω)

1

10

10

10−1

10−4

10−5

−2

10

k=1
k=2
k=3
k=4
O(hk+1 )

−7

10−9

k=1
k=2
k=3
k=4
O(hk+1 )

10−2

100

10−3

10

k[[βuh ]]kL2 (Γh )

0

0

2

k=1
k=2
k=3
k=4
O(hk )

10−4
4

refinements (uniform)

6

0

10−6
10−8

2

4

refinements (uniform)

6

0

2

refinements (uniform)

Figure 4.1: Error convergence in different norms for the circle example.

16

4

6

Conclusion
A new unfitted finite element method with a high order geometrical approximation of level set domains
is presented and discussed in detail. The method is efficient and easy to implement. Numerical examples
reveal its potential in handling complex geometries robust and highly accurate.
As an example for the potential of the method, the method has been combined with a rather standard
Nitsche-XFEM discretization for an unfitted interface problem and optimal order of convergence has been
observed. The potential of the method goes beyond this specific example. As the method is geometry-based
it allows to consider high order methods for a large range of unfitted finite element methods. However,
many open questions remain. Stable discretizations and suitable linear solvers for high order unfitted finite
element methods are difficult to develop and analyze.
The presented method requires further attention. A rigorous error analysis of the presented method is
missing, yet. We will address this in a forthcoming paper. Further, extensions to non-simplex meshes and
time-dependent problems (with moving interfaces) are interesting topics for future studies.
Acknowledgements
The author would like to express his appreciation to Arnold Reusken for his valuable and constructive
feedback on a first draft of this manuscript. The author also greatly appreciates the fruitful discussion with
Joachim Schöberl on the topic of mesh transformations and shape regularity in the context of this study.
References
[1] I. Babuška, The finite element method with penalty, Math. Comp. 27 (122) (1973) 221–228.
[2] J. W. Barrett, C. M. Elliott, Finite element approximation of the Dirichlet problem using the boundary penalty method,
Numer. Math. 49 (4) (1986) 343–366.
[3] R. Glowinski, T.-W. Pan, J. Periaux, A fictitious domain method for Dirichlet problem and applications, Comput. Meth.
Appl. Mech. Eng. 111 (3–4) (1994) 283–303.
[4] E. Burman, P. Hansbo, Fictitious domain finite element methods using cut elements: II. a stabilized Nitsche method,
Applied Numerical Mathematics 62 (4) (2012) 328–341.
[5] C. S. Peskin, D. M. McQueen, A three-dimensional computational method for blood flow in the heart I. immersed elastic
fibers in a viscous incompressible fluid, J. Comput. Phys. 81 (2) (1989) 372–405.
[6] P. Bastian, C. Engwer, An unfitted finite element method using discontinuous Galerkin, International journal for numerical
methods in engineering 79 (12) (2009) 1557–1576.
[7] T.-P. Fries, T. Belytschko, The extended/generalized finite element method: an overview of the method and its applications, International Journal for Numerical Methods in Engineering 84 (3) (2010) 253–304.
[8] A. Hansbo, P. Hansbo, An unfitted finite element method, based on nitsches method, for elliptic interface problems,
Computer methods in applied mechanics and engineering 191 (47) (2002) 5537–5552.
[9] S. Groß, V. Reichelt, A. Reusken, A finite element based level set method for two-phase incompressible flows, Comput.
Visual. Sci. 9 (2006) 239–257.
[10] R. Massjung, An unfitted discontinuous Galerkin method applied to elliptic interface problems, SIAM J. Numer. Anal.
50 (6) (2012) 3134–3162.
[11] R. Becker, E. Burman, P. Hansbo, A Nitsche extended finite element method for incompressible elasticity with discontinuous modulus of elasticity, Computer Methods in Applied Mechanics and Engineering 198 (4144) (2009) 3352 – 3360.
[12] M. A. Olshanskii, A. Reusken, J. Grande, A finite element method for elliptic equations on surfaces, SIAM journal on
numerical analysis 47 (5) (2009) 3339–3358.
[13] W. E. Lorensen, H. E. Cline, Marching cubes: A high resolution 3d surface construction algorithm, in: ACM SIGGRAPH
Computer Graphics, Vol. 21, ACM, 1987, pp. 163–169.
[14] T. A. Nærland, Geometry decomposition algorithms for the Nitsche method on unfitted geometries, Master’s thesis,
University of Oslo (2014).
[15] U. M. Mayer, A. Gerstenberger, W. A. Wall, Interface handling for three-dimensional higher-order XFEM-computations
in fluid–structure interaction, International Journal for Numerical Methods in Engineering 79 (7) (2009) 846–869.
[16] C. Lehrenfeld, The Nitsche XFEM-DG space-time method and its implementation in three space dimensions, SIAM
Journal on Scientific Computing 37 (1) (2015) A245–A270.
[17] C. Lehrenfeld, On a space-time extended finite element method for the solution of a class of two-phase mass transport
problems, Ph.D. thesis, RWTH Aachen (February 2015).
[18] S. Gross, et al., DROPS package for simulation of two-phase flows (2015).
URL http://www.igpm.rwth-aachen.de/DROPS

17

[19] C. Engwer, F. Heimann, Dune-UDG: a cut-cell framework for unfitted discontinuous Galerkin methods, in: Advances in
DUNE, Springer, 2012, pp. 89–100.
[20] E. Burman, S. Claus, P. Hansbo, M. G. Larson, A. Massing, CutFEM: Discretizing geometry and partial differential
equations, International Journal for Numerical Methods in Engineering.
URL http://dx.doi.org/10.1002/nme.4823
[21] Y. Renard, J. Pommier, GetFEM++, an open-source finite element library (2014).
URL http://home.gna.org/getfem
[22] T. Carraro, S. Wetterauer, On the implementation of the eXtended finite element method (XFEM) for interface problems,
arXiv preprint arXiv:1507.04238.
[23] A. Y. Chernyshenko, M. A. Olshanskii, An adaptive octree finite element method for pdes posed on surfaces, Computer
Methods in Applied Mechanics and Engineering 291 (2015) 146–172.
[24] B. Müller, F. Kummer, M. Oberlack, Highly accurate surface and volume integration on implicit domains by means of
moment-fitting, International Journal for Numerical Methods in Engineering 96 (8) (2013) 512–528.
[25] Y. Sudhakar, W. A. Wall, Quadrature schemes for arbitrary convex/concave volumes and integration of weak form in
enriched partition of unity methods, Computer Methods in Applied Mechanics and Engineering 258 (2013) 39–54.
[26] R. Saye, High-order quadrature method for implicitly defined surfaces and volumes in hyperrectangles, SIAM Journal on
Scientific Computing 37 (2) (2015) A993–A1019.
[27] E. Burman, P. Hansbo, M. G. Larson, A cut finite element method with boundary value correction, arXiv preprint
arXiv:1507.03096.
[28] J. Grande, A. Reusken, A higher order finite element method for partial differential equations on surfaces, Tech. Rep. 403,
Institut für Geometrie und Praktische Mathematik, RWTH Aachen (2014).
[29] K. W. Cheng, T.-P. Fries, Higher-order XFEM for curved strong and weak discontinuities, International Journal for
Numerical Methods in Engineering 82 (5) (2010) 564–590.
[30] K. Dréau, N. Chevaugeon, N. Moës, Studied X-FEM enrichment to handle material interfaces with higher order finite
element, Computer Methods in Applied Mechanics and Engineering 199 (29) (2010) 1922–1936.
[31] S. Basting, M. Weismann, A hybrid level set – front tracking finite element approach for fluid–structure interaction and
two-phase flow applications, Journal of Computational Physics 255 (2013) 228 – 244.
[32] P. Oswald, On a BPX-preconditioner for P1 elements, Computing 51 (1993) 125–133.
[33] J. Schöberl, C++11 implementation of finite elements in NGSolve, Tech. Rep. ASC-2014-30, Institute for Analysis and
Scientific Computing (September 2014).
[34] Z. Li, A fast iterative algorithm for elliptic interface problems, SIAM Journal on Numerical Analysis 35 (1995) 230–254.
[35] C. Lehrenfeld, A. Reusken, Optimal preconditioners for Nitsche-XFEM discretizations of interface problems, arXiv preprint
arXiv:1408.2940.

18

