THE AGGREGATED UNFITTED FINITE ELEMENT METHOD FOR ELLIPTIC
PROBLEMS

arXiv:1709.09122v1 [] 26 Sep 2017

SANTIAGO BADIA, FRANCESC VERDUGO, AND ALBERTO F. MARTÍN
Abstract. Unfitted finite element techniques are valuable tools in different applications where the generation of body-fitted meshes is difficult. However, these techniques are prone to severe ill conditioning
problems that obstruct the efficient use of iterative Krylov methods and, in consequence, hinders the practical usage of unfitted methods for realistic large scale applications. In this work, we present a technique
that addresses such conditioning problems by constructing enhanced finite element spaces based on a cell
aggregation technique. The presented method, called aggregated unfitted finite element method, is easy to
implement, and can be used, in contrast to previous works, in Galerkin approximations of coercive problems with conforming Lagrangian finite element spaces. The mathematical analysis of the new method
states that the condition number of the resulting linear system matrix scales as in standard finite elements
for body-fitted meshes, without being affected by small cut cells, and that the method leads to the optimal finite element convergence order. These theoretical results are confirmed with 2D and 3D numerical
experiments.

Keywords: unfitted finite elements; embedded boundary methods; ill-conditioning.
Contents
1. Introduction
2. Embedded boundary setup and cell aggregation
3. Aggregated unfitted Lagrangian finite element spaces
4. Approximation of elliptic problems
5. Numerical analysis
5.1. Stability of the coordinate vector extension matrix
5.2. Mass matrix condition number
5.3. Inverse inequality
5.4. Coercivity and Nitsche’s coefficient
5.5. Well-posedness of the unfitted finite element (FE) problem
5.6. Error estimates
6. Numerical experiments
6.1. Setup
6.2. Moving domain experiment
6.3. Convergence test
7. Conclusions
References

1
3
5
7
9
9
10
10
11
11
13
14
14
15
16
17
19

1. Introduction
Unfitted FE techniques are specially appealing when the generation of body-fitted meshes is difficult.
They are helpful in a number of contexts including multi-phase and multi-physics applications with moving
interfaces (e.g., fracture mechanics, fluid-structure interaction [1], or free surface flows), or in situations
in which one wants to avoid the generation of body-fitted meshes to simplify as far as possible the preprocessing steps (e.g., shape or topology optimization frameworks, medical simulations based on CT-scan
Date: September 27, 2017.
SB gratefully acknowledges the support received from the Catalan Government through the ICREA Acadèmia Research
Program. E-mails: sbadia@cimne.upc.edu (SB), fverdugo@cimne.upc.edu (FV), amartin@cimne.upc.edu (AM).
1

AGGREGATED UNFITTED FINITE ELEMENT METHOD

2

data, or parallel large-scale simulations). In addition, the huge success of isogeometrical analysis (splinebased discretization) and the severe limitations of this approach in complex 3D geometries will probably
increase the interest of unfitted methods in the near future [2]. Unfitted FE methods have been named in
different ways. When designed for capturing interfaces, they are usually denoted as eXtended FE methods
(XFEM) [3], whereas they are usually denoted as embedded (or immersed) boundary methods, when
the motivation is to simulate a problem using a (usually simple Cartesian) background mesh (see, e.g.,
[4, 5, 6]).
Yet useful, unfitted FE methods have known drawbacks. They pose problems to numerical integration,
imposition of Dirichlet boundary conditions, and lead to ill conditioning problems. Whereas different
techniques have been proposed in the literature to address the issues related with numerical integration
(see, e.g., [7]) and the imposition of Dirichlet boundary conditions (see, e.g., [8]), the conditioning problems
are one of the main showstoppers still today for the successful use of this type of methods in realistic large
scale applications. For most of the unfitted FE techniques, the condition number of the discrete linear
system does not only depend on the characteristic element size of the background mesh, but also on
the characteristic size of the cut cells, which can be arbitrary small and have arbitrarily high aspect
ratios. This is an important problem. At large scales, linear systems are solved with iterative Krylov subspace methods [9] in combination with scalable preconditioners. Unfortunately, the well known scalable
preconditioners based on (algebraic) multigrid [10] or multi-level domain decomposition [11] are mainly
designed for body-fitted meshes and cannot readily deal with cut cells. Different preconditioners for
unfitted FE methods have been recently proposed, but they are mainly serial non-scalable algorithms
(see, e.g., [12, 13, 14, 15]). Recently, a robust domain decomposition preconditioner able to deal with
cut cells has been proposed in [16]. Even though this method has proven to be scalable in some complex
3D examples, it is based on heuristic considerations without a complete mathematical analysis and its
application to second (and higher) order FEs is involved. This lack of preconditioners for unfitted FEs
can be addressed with enhanced formulations that provide well-posed discrete systems independently of
the size of the cut cells. Once the conditioning problems related to cut cells are addressed, the application
of standard preconditioners for body-fitted meshes to the unfitted case is strongly simplified, opening the
door to large-scale computations.
The main goal of this work is to develop such an enhanced unfitted FE formulation that fixes the problems associated with cut cells. The goal is to achieve condition numbers that scale only with the element
size of the background mesh in the same way as in standard FE methods for body-fitted meshes. Our
purpose is to implement it in FEMPAR, our in-house large scale FE code [17]. Since FEMPAR is a parallel multi-physics multi-scale code that includes different continuous and discontinuous FE formulations
and several element types, it is crucial for us that the novel formulation fulfills the following additional
properties: 1) It should be general enough to be applied to several problem types, 2) it should deal with
both continuous and discontinuous FE formulations, 3) it should deal with high order interpolations, and
4) it should be easily implemented in an existing parallel FE package.
To our best knowledge, none of the existing unfitted FE formulations fulfill these requirements simultaneously. For instance, one can consider the ghost penalty formulation used in the CutFEM method [4, 18]
However, it leads to a weakly non-consistent algorithm, and it requires to compute high order derivatives
on faces for high order FEs, which are not at our disposal in general FE codes and are expensive to
compute, certainly complicating the implementation of the methods and harming code performance. Alternatively, for finite volume and discontinuous Galerkin (DG) formulations, one can consider the so-called
cell aggregation (or agglomeration) techniques [19, 20]. E.g., for DG formulations, the idea is simple: cells
with the small cut cell problem, i.e., the ratio between the volume of the cell inside the physical domain
and the total cell volume is close to zero, are merged with neighbor full cells forming aggregates. A new
polynomial space is defined in each aggregate that replaces the local FE spaces of all cells merged in it.
This process fixes the conditioning problems, since the support of the newly defined shape functions is at
least the volume of a full cell. Even though this idea is simple and general enough to deal with different
problem types and high order interpolations, the resulting discrete spaces are such that the enforcement
of continuity through appropriate local-to-global degrees of freedom (DOFs) numbering, as in standard
FE codes (see, e.g., [17]), is not possible, limiting their usage to discontinuous Galerkin or finite volume
formulations. Up to our best knowledge, there is no variant of cell agglomeration currently proposed

AGGREGATED UNFITTED FINITE ELEMENT METHOD

3

in the literature producing conforming FE spaces, which could be used for classical continuous Galerkin
formulations. It is the purpose of this work.
In this article, we present an alternative cell aggregation technique that can be used for both continuous
and discontinuous formulations, the aggregated unfitted FE method. We start with the usual (conforming)
Lagrangian FE space that includes cut cells, which is known to lead to conditioning problems. The main
idea is to eliminate from this space all the potentially problematic DOFs by introducing a set of judiciously
defined constraints. These constraints are introduced using information provided by the cell aggregates,
without altering the conformity of the original FE space. Alternatively, the method can be understood
as an extension operator from the interior (well-posed ) FE space that only involves interior cells to a
larger FE space that includes cut cells and covers the whole physical domain. Discontinuous spaces can
also be generated as a particular case of this procedure, which makes the method compatible also with
DG formulations. In contrast to previous works, we also include a detailed mathematical analysis of the
method, in terms of well-posedness, condition number estimates, and a priori error estimates. For elliptic
problems, we mathematically prove that 1) the method leads to condition numbers that are independent
from small cut cells, 2) the condition numbers scale with the size of the background mesh as in the standard
FE method, 3) the penalty parameter of Nitsche’s method required for stability purposes is bounded above,
and 4) the optimal FE convergence order is recovered. These theoretical results are confirmed with 2D
and 3D numerical experiments using the Poisson equation as a model problem.
The outline of the article is as follows. In Section 2, we introduce our embedded boundary setup and
the strategy to build the cell aggregates. In Section 3, we describe the construction of the novel FE spaces
based on the cell aggregates. In Section 4, we introduce our elliptic model problem. The numerical analysis
of the method is carried out in Section 5. Finally, we present a complete set of numerical experiments in
Section 6 and draw some conclusions in Section 7.
2. Embedded boundary setup and cell aggregation
Let Ω ⊂ Rd be an open bounded polygonal domain, with d ∈ {2, 3} the number of spatial dimensions.
For the sake of simplicity and without loss of generality, we consider in the numerical experiments below
.
that the domain boundary is defined as the zero level-set of a given scalar function ψ ls , namely ∂Ω = {x ∈
Rd : ψ ls (x) = 0}.1 We note that the problem geometry could be described using 3D CAD data instead of
level-set functions, by providing techniques to compute the intersection between cell edges and surfaces
(see, e.g., [21]). In any case, the way the geometry is handled does not affect the following exposition. Like
in any other embedded boundary method, we build the computational mesh by introducing an artificial
domain Ωart such that it has a simple geometry that is easy to mesh using Cartesian grids and it includes
the physical domain Ω ⊂ Ωart (see Fig. 1a).

internal cells
cut cells
external cells

(a)

(b)

Figure 1. Embedded boundary setup.
Let us construct a partition of Ωart into cells, represented by Thart , with characteristic cell size h. We
are interested in Thart being a Cartesian mesh into hexahedra for d = 3 or quadrilaterals for d = 2, even
1Analogous assumption have to be made for body-fitted methods.

AGGREGATED UNFITTED FINITE ELEMENT METHOD

4

though unstructured n-simplex background meshes can also be considered. Cells in Thart can be classified as
follows: a cell K ∈ Thart such that K ⊂ Ω is an internal cell ; if K ∩ Ω = ∅, K is an external cell ; otherwise,
K is a cut cell (see Fig. 1b). The set of interior (resp., external and cut) cells is represented with Thin and
its union Ωin ⊂ Ω (resp., (Thext , Ωext ) and (Thcut , Ωcut )). Furthermore, we define the set of active cells as
.
Thact = Thin ∪ Thcut and its union Ωact . In the numerical analysis, we assume that the background mesh is
quasi-uniform (see, e.g., [22, p.107]) to reduce technicalities, and define a characteristic mesh size h. The
maximum element size is denoted with hmax .
We can also consider non-overlapping cell aggregates AK composed of cut cells and one interior cell
K such that the aggregate is connected, using, e.g., the strategy described in Algorithm 2.1. It leads to
another partition Thagg defined by the aggregations of cells in Thact ; interior cells that do not belong to any
aggregate remain the same. By construction of Algorithm 2.1, there is only one interior cell per aggregate,
denoted as the root cell of the aggregate, and every cut cell belongs to one and only one aggregate. For a
cut cell, we define its root cell as the root of the only aggregate that contains the cut cell. The root of an
interior cell is the cell itself. Thus, there is a one-to-one mapping between aggregates (including interior
cells) AK ∈ Thagg and the root cut cell K ∈ Thin . As a result, we can use the same index for the aggregate
and the root cell. We build the aggregates in Thagg with Algorithm 2.1. In any case, other aggregation
algorithms could be considered, e.g., touching in the first step of the algorithm not only the interior cells,
but also cut cells without the small cut cell problem. It can be implemented by defining the quantity
.
ηK = |K∩Ω|
|K| and touch in the first step not only the interior cells but also any cut cell with ηK > η0 > 0
for a fixed value η0 .
Algorithm 2.1 (Cell aggregation scheme).
(1) Mark all interior cells as touched and all cut cells as untouched.
(2) For each untouched cell, if there is at least one touched cell connected to it through a facet F such
that F ∩ Ω 6= ∅, we aggregate the cell to the touched cell belonging to the aggregate containing
the closest interior cell. If more than one touched cell fulfills this requirement, we choose one
arbitrarily, e.g., the one with smaller global id.
(3) Mark as touched all the cells aggregated in 2.
(4) Repeat 2. and 3. until all cells are aggregated.
touched

(a) Step 1

untouched

(b) Step 2

Aggregates’ boundary

(c) Step 3

∂Ω

(d) Step 4

Figure 2. Illustration of the cell aggregation scheme defined in Algorithm 2.1.
Fig. 2 shows an illustration of each step in Algorithm 2. The black thin lines represent the boundaries
of the aggregates. Note that from step 1 to step 2, some of the lines between adjacent cells are removed,
meaning that the two adjacent cells have been merged in the same aggregate. The aggregation schemes
can be easily applied to arbitrary spatial dimensions. As an illustrative example, Fig. 3 shows some of
the aggregates obtained for a complex 3D domain.
In the forthcoming sections, we need an upper bound of the size of the aggregates generated with
Algorithm 2.1. To this end, let us consider the next lemma.

AGGREGATED UNFITTED FINITE ELEMENT METHOD

5

Figure 3. 3D aggregates.
Lemma 2.2. Assume that from any cut cell K0 ∈ Thact there is a cell path {K0 , K1 , . . . , Kn } that satisfies:
1) two consecutive cells share a facet F such that F ∩ Ω 6= ∅; 2) Kn is an interior cell; 3) n ≤ γmax , where
γmax is a fixed integer. Then, the maximum aggregate size is at most (2γmax + 1)hmax .
Proof. By construction, an aggregate can grow at most at a rate of one layer of elements per each iteration.
Thus, after n iterations the aggregate size will be at most (2n + 1)hmax considering that the aggregate can
potentially grow in all spatial directions. It is obvious to see that the aggregation scheme finishes at most
after γmax iterations. Thus, the aggregate size will be less or equal than (2γmax + 1)hmax .

From Lemma 2.2, it follows that the aggregate size will be bounded if so is the value of γmax . In what
follows, we assume that γmax is fixed, e.g., eliminating any cut cell that would violate property 3) in
Lemma 2.2. One shall assume that each cut cell shares at least one corner with an interior cell (this is
usually true if the grid is fine enough to capture the geometry). In this situation, we can easily see that
γmax = 2 for 2D and γmax = 3 for 3D. Then, by Lemma 2.2, the aggregate size is at most 5hmax in 2D and
7hmax in 3D. Even though it is not used in the proof of Lemma 2.2, the fact that we aggregate cut cells
to the touched cells belonging to the aggregate containing the closest interior cell (see step 2 in Algorithm
2.1) contributes to further reduce the aggregate size. Indeed, the actual size of the aggregates generated
in the numerical examples (cf. Section 6) is much lower than the predicted by these theoretical bounds.
In 2D, the aggregate size tends to 2hmax as the mesh is refined, whereas it tends to 3hmax in the 3D case.
This shows that the aggregation scheme produces relative small aggregates in the numerical experiments.
3. Aggregated unfitted Lagrangian finite element spaces
Our goal is to define a FE space using the cell aggregates introduced above. To this end, we need to
.
introduce some notation. In the case of n-simplex meshes, we define the local FE space V (K) = Pq (K),
i.e., the space of polynomials of order less or equal to q in the variables x1 , . . . , xd . For n-cube meshes,
.
V (K) = Qq (K), i.e., the space of polynomials that are of degree less or equal to k with respect to each
variable x1 , . . . , xd . In this work, we consider that the polynomial order q is the same for all the cells in
the mesh. We restrict ourselves to Lagrangian FE methods. Thus, the basis for V (K) is the Lagrangian
basis (of order q) on K. We denote by N (K) the set of Lagrangian nodes of order q of cell K. There is a
one-to-one mapping between nodes a ∈ N (K) and shape functions φa (x); it holds φa (xb ) = δab , where xb
are the space coordinates of node b. We assume that there is a local-to-global DOF map such that the
resulting global system is C 0 continuous. This process can be elaborated for hp-adaptivity as well, but it
is not the purpose of this work.
With this notation, we can introduce the active FE space associated with the active portion of the
background mesh
.
Vhact = {v ∈ C 0 (Ωact ) : v|K ∈ V (K), for any K ∈ Thact }.
We could analogously define the interior FE space
.
Vhin = {v ∈ C 0 (Ωin ) : v|K ∈ V (K), for any K ∈ Thin }.

AGGREGATED UNFITTED FINITE ELEMENT METHOD

6

The active FE space Vhact (see Fig. 4c) is the functional space typically used in unfitted FE methods (see,
e.g., [16, 15, 6]). It is well known that Vhact leads to arbitrary ill conditioned systems when integrating the
FE weak form on the physical domain Ω only (if no extra technique is used to remedy it). It is obvious
that the interior FE space Vhin (see Fig. 4a) is not affected by this problem, but it is not usable since it
is not defined on the complete physical domain Ω. Instead, we propose an alternative space Vhagg that is
defined on Ω but does not present the problems related to Vhact . We can define the set of nodes of Vhin and
.
Vhact as Nhin and Nhact , respectively (see Fig. 4). We define the set of outer nodes as Nhout = Nhact \ Nhin
(marked with red crosses in Fig. 4b). The outer nodes are the ones that can lead to conditioning problems
due to the small cut cell problem (see (9)). The space Vhagg is defined taking as starting point Vhact , and
adding judiciously defined constraints for the nodes in Nhout .

◦ nodes in Nhin
• nodes in Nhact
× nodes in Nhout

(a) Vhin

(b) Vhagg

(c) Vhact

Figure 4. Finite Element spaces.
In order to define Vhagg we observe that, in nodal Lagrangian FE spaces, there is a one-to-one map
between DOFs and nodes (points) of the FE mesh (for vector spaces, the same is true for every component
of the vector field). On the other hand, we can define the owner vertex, edge, or face (VEF) of a node
as the lowest-dimensional VEF that contains the node. Furthermore, we can construct a map that for
every VEF F such that F 6⊂ Ω, gives a cell owner among all the cells that contain it. This map can be
arbitrarily built. E.g., we can consider as cell owner the one in the smallest aggregate. As a result, we
have a map between DOFs and (active) cells. Every active cell belongs to an aggregate, which has its own
root (interior) cell. So, we also have a map between DOFs and interior cells. This map between b ∈ Nhout
and the corresponding interior cell is represented with K(b) (see Fig. 5).

×

aggregate
node in Nhout
node to cell map

Figure 5. Map from outer nodes to interior cells.
The space of global shape functions of Vhin and Vhact can be represented as {φb : b ∈ Nhin } and
: b ∈ Nhact }, respectively. Functions in these FE spaces are uniquely represented by their nodal
in
values. We represent the nodal values of uh ∈ Vhin as uin ∈ R|Nh | , whereas the nodal values of uh ∈ Vhact
act
as uact ∈ R|Nh | . Considering, without loss of generality, that the interior nodal values are labeled the
out
same way for both FE spaces, we have that uact = [uin , uout ]T , where uout ∈ R|Nh | .
{φb

AGGREGATED UNFITTED FINITE ELEMENT METHOD

7

Now, we consider the following extension operator. Given uh ∈ Vhin and the corresponding nodal values
we compute the outer nodal values as follows:
X
φa (xb )uin a , for b ∈ Nhout .
uout b =
(1)

uin ,

a∈N (K(b))

That is, the value at an outer node b ∈ Nhout is computed by extrapolating the nodal values of the interior
cell K(b) associated with it. In compact form, we can write it as uout = Cuin , where C is the global
in
act
matrix of constraints. We define the global extension matrix E : R|Vh | → R|Vh | as Euin = [uin , Cuin ]T .
Let us also define the extension operator E : Vhin → Vhact , such that, given uh ∈ Vhin represented by its
nodal values uin , provides the FE function E(uh ) ∈ Vhact with nodal values Euin . We define the range of
.
this operator as Vhagg = range(E(Vhin )) ⊂ Vhact . This FE space is called the aggregated FE space since the
map K(·) between outer nodes and interior cells is defined using the aggregates in Thagg . The motivation
behind the construction of such space is to have a FE space covering Ωact (and thus Ω) with optimal
approximability properties and without the ill-conditioning problems of Vhact .
As one can observe, the new space is defined only by interior nodal values, whereas the conflictive outer
nodes are eliminated via the constraints in (1). These constraints are cell-wise local. Thus, they can be
readily applied at the assembly level in the cell loop, making its implementation very simple, even for nonadaptive codes that cannot deal with non-conforming meshes. We consider as basis for Vhagg the extension
of the shape functions of Vhin , i.e., {E(φa )}a∈N in . The fact that it is a basis for Vhagg is straightforward,
h
due to the fact that the extension operator is linear. The extension of a shape function is easily computed
as follows:
X
E(φa ) = φa +
Cba φb , for a ∈ Nhin ,
b∈C(a)

where C(a) represents the set of outer nodes in Nhout that are constrained by a.
Remark 3.1. We note that one could consider an alternative aggregated space,
Vhagg,∗ = {v ∈ C 0 (Ω) : v|A ∈ V (A), for any A ∈ Thagg },
where V (A) denotes the space of q order Lagrangian polynomials on n-simplices or n-cubes. It is obvious
to check that in fact Vhagg,∗ ⊂ Vhagg , but it is not possible to implement the inter-element continuity for this
space using standard FE techniques. On the other hand, the FE space Vhagg has the same size as the interior
problem and the implementation in existing FE codes requires minimal modifications. Furthermore, it is
also easy to check that the two approaches coincide for DG formulations, where all DOFs belong to the
cells. In fact, a DG method with Vhagg,∗ has been proposed in [20].
4. Approximation of elliptic problems
For the sake of simplicity, we consider the Poisson equation with constant physical diffusion as a model
problem, even though the proposed ideas apply to any elliptic problem with H 1 -stability, e.g., the linear elasticity problem and heterogeneous problems. The Poisson equation with Dirichlet and Neumann
boundary conditions reads as (after scaling with the diffusion term): find u ∈ H 1 (Ω) such that
− ∆u = f

in Ω,

u = gD

on ΓD ,

∇u · n = g N

on ΓN ,

(2)

where (ΓD , ΓN ) is a partition of the domain boundary (the Dirichlet and Neumann boundaries, respectively), f ∈ H −1 (Ω), g D ∈ H 1/2 (ΓD ), and g N ∈ H −1/2 (ΓN ).
For the space discretization, we consider H 1 -conforming FE spaces on the conforming mesh Thact that
are not necessary aligned with the the physical boundary ∂Ω. For simplicity, we assume that, for any cut
cell K ∈ Thact , either K ∩ Γ ⊂ ΓD or K ∩ Γ ⊂ ΓN . We consider both the usual FE space Vhact as well as
the new aggregated space Vhagg in order to compare their properties. We will simply use Vh when it is not
necessary to distinguish between Vhact and Vhagg .
For unfitted grids, it is not clear to include Dirichlet conditions in the approximation space in a strong
manner. Thus, we consider Nitsche’s method [23, 24] to impose Dirichlet boundary conditions weakly on

AGGREGATED UNFITTED FINITE ELEMENT METHOD

8

ΓD . It provides a consistent numerical scheme with optimal converge rates (also for high-order elements)
that is commonly used in the embedded boundary community [6]. We define the FE-wise operators:
Z
Z
.
(τK uv − v (n · ∇u) − u (n · ∇v)) dS,
∇u · ∇v dV +
AK (u, v) =
ΓD ∩K
K∩Ω
Z

.
`K (v) =
τK vg D − (n · ∇v) g D dS,
ΓD ∩K

defined for a generic cell K ∈ Thact . Vector n denotes the outwards normal to ∂Ω. The bilinear form
AK (·, ·) includes the usual form resulting from the integration by parts of (2) and the additional term
associated with the weak imposition of Dirichlet boundary conditions with Nitsche’s method. The righthand side operator `K (·) includes additional terms related to Nitsche’s method. The coefficient τK > 0 is
a mesh-dependent parameter that has to be large enough to ensure the coercivity of AK (·, ·).
The global FE operator A : Vh → Vh 0 and right-hand side term ` ∈ Vh 0 are stated as the sum of the
element contributions, i.e.,
. X
. X
A(u, v) =
AK (u, v), `(v) =
`K (v),
for u, v ∈ Vh .
(3)
K∈Thact

K∈Thact

We will make abuse of notation, using the same symbol for a bilinear form, e.g., A : Vh → Vh 0 ,
.
and its corresponding linear operator, i.e., hAu, vi = A(u, v). Furthermore, we define b : Vh 0 → Vh as
.
b(v) = f (v) + g N (v) + `(v), for v ∈ Vh . With this, the global problem can be stated as: find uh ∈ Vhagg
such that
(4)
A(uh , vh ) = b(vh ),
for any vh ∈ Vhagg .
By definition, this problem can analogously be stated as: find uh ∈ Vhin such that A(E(uh ), E(vh )) =
b(E(vh )) for any vh ∈ Vhagg . After the definition of the FE basis (of shape functions) that spans Vhagg , or
alternatively the extension operator E(·), the previous problem leads to a linear system to be solved.
A sufficient (even though not necessary) condition for A to be coercive is to enforce the element-wise
constant coefficient τK to satisfy
BK (v, v)
.
τK ≥ CK = sup
,
(5)
v∈V (K) DK (v, v)
for all the mesh elements K ∈ Thact intersecting the boundary ΓD . In the previous formula, DK (·, ·) and
BK (·, ·) are the forms defined as
Z
Z
.
.
DK (u, v) =
∇u · ∇v dV, and BK (u, v) =
(n · ∇u) (n · ∇v) dS.
K∩Ω

ΓD ∩K

Since Vh is finite dimensional, and DK (·, ·) and BK (·, ·) are symmetric and bilinear forms, the value CK
.
(i.e., the minimum admissible coefficient τK ) can be computed numerically as CK = λ̃max , being λ̃max the
largest eigenvalue of the generalized eigenvalue problem (see [15] for details): find uK ∈ Vh |K and λ̃ ∈ R
such that
BK (uK , vK ) = λ̃DK (uK , vK ) for all vK ∈ Vh |K .
(6)
For standard FEs for body-fitted meshes, it is enough to compute coefficient τK as τK = β/h to
satisfy condition (5), where β is a sufficiently large (mesh independent) positive constant (see, e.g., [25]).
However, for standard unfitted FE methods using the usual space Vhact without any additional stabilization,
coefficient τK cannot be computed a priori; in fact, the minimum cell-wise value that assures coercivity is
not bounded above. In this case, a value for τK ensuring coercivity has to be computed for each particular
setup using the cell-wise eigenvalue problem (6). The introduction of the new space Vhagg solves this
problem and τK is bounded again in terms of the element size as expected in the body-fitted case (see
Section 5.4 for more details). In this case, we have taken τK = 100/hK in the numerical experiments
below.
The linear system matrix that arises from (4) can be defined as
.
Aab = A(E(φa ), E(φb )),
for a, b ∈ Nhin .
(7)

AGGREGATED UNFITTED FINITE ELEMENT METHOD

The mass matrix related to the aggregated FE space Vhagg is analogously defined as
Z
.
Mab =
E(φa )E(φb ),
for a, b ∈ Nhin .

9

(8)

Ω

It is well known that the usual FE space Vhact is associated with conditioning problems due to cut cells.
The condition number of the discrete system without the aggregation, i.e., considering Vhact instead of
Vhagg in (4), scales as
−(2q+1−2/d)
κ(A) ∼ min ηK
,
(9)
K∈Thact

where κ(A) is the 2-norm condition number of A (see [15] for details). Thus, arbitrarily high condition
numbers are expected in practice since the position of the interface cannot be controlled and the value ηK
can be arbitrarily close to zero. This problem is solved if the new aggregated space Vhagg is used instead
of Vhact (cf. Corollary 5.9).
5. Numerical analysis
In this section, we analyze the well-posedness of the agregated unfitted FE method (4), the condition
number of the arising linear system, and a priori error estimates. As commented above, we assume
that the background mesh is quasi-uniform. Therefore, the number of neighboring cells of a given cell is
bounded above by a constant ncell independently of h. In a mesh refinement analysis, we also assume that
the coarser mesh level-set function already represents the domain boundary.
In the following analysis, all constants being used are independent of h and the location of the cuts in
cells, i.e., ηK . They may also depend on the threshold η0 in the aggregation algorithm if considered; we
have considered η0 = 1 for simplicity. The constants can depend on the shape/size of Ω and ΓD , the order
of the FE space, and the maximum aggregation distance γmax , which are assumed to be fixed in this work.
In turn, due to Lemma 2.2, the maximum size of an aggregate is bounded by a constant times h. As a
result, the following results are robust with respect to the so-called small cut cell problem. When we have
that A ≤ cB for a positive constant c, we may use the notation A . B; analogously for &.
For the analysis below, we need to introduce some extra notation. Given a function uh ∈ Vhin (or Vhact ),
the nodal vector u will be used without any superscript, as soon as it is clear from the context. For
a given cell K, the cell-wise coordinate vector is represented with uK . On the other hand, given a FE
function uh ∈ Vhin , for every interior cell K ∈ Thin , let us define define the cell-wise extension operator
EAK u = [uK , CÂK uK ], where CÂK is the cell-wise constraint matrix, whose entries can be computed in
P
the reference space (see (1)), such that Cu · Cu = K∈T in CÂK uK · CÂK uK . We denote with k · k2 the
h
Euclidean norm of a vector and the induced matrix norm. Standard notation is used to define Sobolev
spaces (see, e.g., [26]). Given a Sobolev space X, its corresponding norm is represented with k · kX .
5.1. Stability of the coordinate vector extension matrix. We start the analysis of the scheme by
proving bounds for the norm of the global and cell-wise coordinate vector extension matrix. Therefore,
their norms can be bounded independently of the cut location and the size of the aggregate.
Lemma 5.1. The cell-wise and global coordinate vector extension matrices hold the following bounds:
1 ≤ kEAK k22 ≤ 1 + kCÂK k22 ,

for every AK ∈ Thagg ,

and

1 ≤ kEk22 ≤ 1 + kCk22 ≤ Ce ,

for a positive constant Ce .
Proof. Using the definition of the extension operator in Section 3, we have that kEuk22 = kuk22 + kCuk22 .
We proceed analogously for the cell-wise result, to get kEAK uk22 = kuk22 + kCÂK uk22 . It proves the first
result. On the other hand, we have,
X
X
kCuk22 =
kCÂK uK k22 ≤
kCÂK k22 kuK k22 ≤ ncell sup kCÂK k22 kuk22 ,
AK ∈Thagg

AK ∈Thagg

AK ∈Thagg

where we have used the fact that the constraint matrix is aggregate-wise and that the maximum number of
cell neighbors of a vertex/edge/face is bounded above by a constant ncell . The value supAK ∈T agg kCÂK k22
h
(or an upper bound) can explicitly be computed prior to the numerical integration and its entries are
independent of the aggregate cut and the geometrical mapping, i.e., h. In fact, given a polynomial order

AGGREGATED UNFITTED FINITE ELEMENT METHOD

10

and γmax , one can precompute the maximum value of kCÂK k22 among all possible aggregate configurations
and explicitly obtain an upper bound Ce of the global extension matrix norm. It proves the lemma. 
5.2. Mass matrix condition number. In order to provide a bound for the condition number of the
mass matrix, we rely on the maximum and minimum eigenvalues of the local mass matrix in the reference
cell K̂:
λ− kuK k22 ≤ kuh k2L2 (K̂) ≤ λ+ kuK k22 ,

for uh ∈ V (K).

The values of λ− and λ+ only depend on the order of the FE space and can be computed for different
orders on n-cubes or n-simplices (see [27]). Using typical scaling arguments, one has the following bound
for the local mass matrix of the physical cell:
λ− hdK kuK k22 ≤ kuh k2L2 (K) ≤ λ+ hdK kuK k22 .

(10)

In the next lemma, we prove the equivalence between the L2 (Ωact ) norm and the interior DOF Euclidean
norm, for functions in Vhagg .
Lemma 5.2. The following bounds hold:
hd kuk22 . kE(uh )k2L2 (Ωact ) . hd kuk22 ,

for any uh ∈ Vhin .

(11)

Proof. By definition, every function in Vhagg can be expressed as E(uh ) for some uh ∈ Vhin . Using (10), the
fact that Ωin ⊂ Ω, and the quasi-uniformity of the background mesh, we obtain the lower bound in (11)
as follows:
X
X
kE(uh )k2L2 (Ωact ) ≥ kuh k2L2 (Ωin ) =
kuh k2L2 (K) ≥
hdK λ− kuK k22 & hd kuk22 .
K∈Thin

K∈Thin

On the other hand, using Ω ⊂ Ωact , Lemma 5.1, (10), and the fact that the number of surrounding cells
of a node is bounded above by a positive constant, we get:
X
X
kE(uh )k2L2 (Ωact ) =
kE(uh )k2L2 (K) ≤
hdK λ+ kuK k22 . hd kEuk22 . hd kuk22 .
K∈Thact

K∈Thact



It proves the lemma.

The upper and lower bounds in (11) lead to the continuity of the extension operator and a bound for
the condition number of the mass matrix of the aggregated FE space.
Corollary 5.3 (Continuity of the extension operator). The extension operator satisfies the following
bound:
kE(uh )k2L2 (Ωact ) . kuh k2L2 (Ωin ) ,

for any uh ∈ Vhin .

Corollary 5.4 (Mass matrix condition number). The mass matrix M in (8), related to the aggregated FE
space Vhagg , is bounded by κ(M) ≤ C, for a positive constant C.
5.3. Inverse inequality. In order to prove the condition number bound for the system matrix arising
from (4), we need to prove first an extended inverse inequality. We rely on the fact that an inverse
inequality holds for the FE space Vhact , i.e.,
k∇uh kL2 (Ωact ) . h−1 kuh kL2 (Ωact ) ,

for any uh ∈ Vhact .

(12)

This standard result for conforming meshes can be found, e.g., in [22, p. 111].
Lemma 5.5 (Inverse inequality). The following inverse inequality holds:
k∇E(uh )kL2 (Ωact ) . h−1 kuh kL2 (Ωin ) ,

for any uh ∈ Vhin .

Proof. Using the fact that Ωin ⊆ Ω ⊆ Ωact , E(uh ) ∈ Vhact , the standard inverse inequality (12), and the
stability of the extension operator in Lemma 5.3, we get:
k∇E(uh )kL2 (Ωact ) . h−1 kE(uh )kL2 (Ωact ) . h−1 kuh kL2 (Ωin ) .
It proves the lemma.



AGGREGATED UNFITTED FINITE ELEMENT METHOD

11

5.4. Coercivity and Nitsche’s coefficient. In this section, we consider a trace inequality that is needed
to prove the coercivity of the bilinear form in (3). Given a cell K ∈ Thact , let us consider the set of
constraining interior cells K1 , . . . , KmK , mK ≥ 1, i.e., the interior S
cells that constraint at least one DOF
.
.
K
of the cut cell. Let us also define Kcut = K ∩ Ω and ΩK = Kcut ∪ m
i=1 Ki ⊂ Ω.
Lemma 5.6. For any uh ∈ Vhagg and K ∈ Thact , the following bound holds
−1

kn · ∇uh kL2 (ΓD ∩K) ≤ C∂ hK 2 k∇uh kL2 (ΩK ) ,
for a positive constant C∂ .
Proof. For interior cells, the left-hand side is zero and the bound trivially holds. Let us consider a cut cell
.
K. Let us also consider a FE function uh ∈ Vhagg and its gradient ξ h = ∇uh . Assuming that all the cells
have the same order, we have that ξ h belongs to the discontinuous Lagrangian FE space of order q − 1,
and we represent the corresponding coordinate vector with ξ K .
First, we use the equivalence of norms in finite dimension and a scaling argument to get:
kξ h k2L2 (ΓD ∩K) . |ΓD ∩ K|kξ h k2L∞ (K) ,
where the constant can only depend on the FE space order. Analogously, we have kξ h k2L∞ (K̂) . kξ K k22 .
Following the same ideas as above, ξ can be expressed as an extension of the corresponding nodal values of
the q − 1 order FE spaces on top of the interior cells Ki , represented with ξ K ; we represent this extension
i

with the matrix DK , i.e., ξ K = DK [ξ 1 , . . . , ξ m ]T . Using an analogous reasoning as above for matrix C,
P K
2
the norm of this matrix cannot depend on the cut or h. Thus, we have that kξ K k22 . m
i=1 kξ Ki k2 . On the
other hand, using again the equivalence of norms in finite dimension, we get kξ K k2 . kξ h kL2 (K̂i ) . As a
i
result, using typical scaling arguments, and using the fact that |K| . |Ki | . |K| for constants independent
of mesh size or order, we get:
mK
X
2
−1
kξ h kL∞ (K) . |K|
kξ h k2L2 (Ki ) .
i=1

Combining these results, we get:
kn ·

ξ h k2L2 (Γ∩AK )

≤ |Γ ∩

K||K|kξ h k2L∞ (K)

.

h−1
K

mK
X

kξ h k2L2 (Ki ) ,

i=1

where we have used the fact that |Γ ∩
lemma.

K||K|−1

.

h−1
K

holds for a quasi-uniform mesh. It proves the


5.5. Well-posedness of the unfitted FE problem. In this section, we prove coercivity and continuity
of the bilinear form (3). First, we prove coercivity with respect to the following mesh dependent norm in
Vhagg :
X
.
2
|||uh |||2h = k∇uh k2L2 (Ω) +
βK h−1
for uh ∈ Vhagg ,
K kuh kL2 (ΓD ∩K) ,
K∈Thact

which is next proved to bound the L2 (Ω) norm.
Theorem 5.7. The aggregated unfitted FE problem in (4) satisfies the following bounds:
i) Coercivity:
A(uh , uh ) & |||uh |||2h , for any uh ∈ Vhagg ,
(13)
ii) Continuity:
A(uh , vh ) . |||uh |||h |||vh |||h , for uh , vh ∈ Vhagg ,
(14)
if βK > C, for some positive constant C. In this case, there exists one and only one solution of (4).
Proof. For cut cells, we use
Z
−1
2
−1
2
uh (n · ∇uh ) ≤ αK C∂ h−1
hK kn · ∇uh k2L2 (ΓD ∩K)
K kuh kL2 (ΓD ∩K) + αK C∂
ΓD ∩K

−1
2
2
≤ αK C∂ h−1
K kuh kL2 (ΓD ∩K) + αK k∇uh kL2 (ΩK ) .

(15)

AGGREGATED UNFITTED FINITE ELEMENT METHOD

12

Using the fact that the mesh is quasi-uniform and that the number of neighboring cells and γmax is
bounded, one can take a value for αK large enough (but uniform with respect to h and the cut location)
such that:
Z
X
1
2
2
2
uh (n · ∇uh ) ≤
αK C∂ h−1
K kuh kL2 (ΓD ∩K) + k∇uh kL2 (Ω) .
2
ΓD
act
K∈Th

As a result, we get:
X
1
2
(βK − αK C∂ ) h−1
A(uh , uh ) ≥ k∇uh k2L2 (Ω) +
K kuh kL2 (ΓD ∩K) .
2
act
K∈Th

For, e.g., βK > 2αK C∂ , A(uh , uh ) is a norm. By construction, this lower bound for βK is independent of
the mesh size h and the intersection of Γ and Thact . It proves the coercivity property in (13). Thus, the
bilinear form is non-singular. The continuity in (14) can readily be proved by repeated use of the CauchySchwarz inequality and inequality (15). Since the problem is finite-dimensional and the corresponding
linear system matrix is non-singular, there exists one and only one solution of this problem.

Lemma 5.8. If Ω has smoothing properties, the following bound holds:
for any uh ∈ Vhagg .

kuh kL2 (Ω) . |||uh |||h ,

Proof. Let us consider uh ∈ Vhagg and let ψ ∈ H01 (Ω) solve the problem −∆ψ = uh with the boundary
conditions ψ = 0 on ΓD and n·∇ψ = 0 on ΓN . Using the fact that the domain Ω has smoothing properties,
it holds kψkH 2 (Ω) . kuh kL2 (Ω) . We have, after integration by parts:
Z
Z
Z
kuh kL2 (Ω) = − uh ∆ψ =
∇uh · ∇ψ −
uh n · ∇ψ.
(16)
Ω

Ω

ΓD

The first term in the right-hand side of (16) is easily bounded using the Cauchy-Schwarz inequality:
Z
∇uh · ∇ψ ≤ k∇uh kL2 (Ω) k∇ψkL2 (Ω) . k∇uh kL2 (Ω) kuh kL2 (Ω) .
Ω

On the other hand, the following trace inequality holds kn · ∇ψk2L2 (ΓD ) . |ψ|2H 2 (Ω) for a constant that
depends on the size of ΓD (see [26]). Using the Cauchy-Schwarz inequality and the previous trace inequality,
we readily get:

1 
1
2
2
Z
X
X
−1
2
2
 

uh n · ∇ψ ≤ 
h kuh k 2
−
hK kn · ∇ψk 2
K

ΓD

L (ΓD ∩K)

L (ΓD ∩K)

K∈Thact

K∈Thact

1


≤

2

X

2

h−1
K kuh kL2 (ΓD ∩K)

kn · ∇ψkL2 (ΓD )

K∈Thact

1


.

2

X

2

h−1
K kuh kL2 (ΓD ∩K)

kuh kL2 (Ω) .

K∈Thact

Combining these bounds, we prove the lemma.



Corollary 5.9 (Stiffness matrix condition number). The condition number of the linear system matrix A
in (7) is bounded by κ(A) . h−2 .
Proof. To prove the corollary, we have to bound u · Au = A(E(uh ), E(uh )) above and below by kuk22
times some constant. The lower bound follows from the coercivity property in Th. 5.7, Lemma 5.8, the
lower bounds in Lemmas 5.2 and 5.1, which lead to A(E(uh ), E(uh )) & kE(uh )k2L2 (Ω) & hd kuk22 . The upper
bound is readily obtained from the continuity property in Lemma 5.7 and the upper bound in Lemma 5.8,
i.e., u · Au = A(E(uh ), E(uh )) . |||E(uh )|||2h . Using scaling arguments and the equivalence of norms for
2
finite-dimensional spaces, we get kE(uh )k2L2 (ΓD ∩K) . hd−1
K kEuK k2 . Adding up for all cells, invoking the

AGGREGATED UNFITTED FINITE ELEMENT METHOD

13

fact that the number of neighbour cells is bounded, and using the upper bound of the coordinate vector
extension operator in 5.8, we obtain:
X
d−2
2
βK h−1
kuk22 .
(17)
K kE(uh )kL2 (ΓD ∩K) . h
K∈Thact

Using the inverse inequality in Lemma 5.5 and the upper bound in Lemma 5.2, we obtain:
k∇E(uh )k2L2 (Ω) . h−2 kE(uh )k2L2 (Ω) . hd−2 kuk22 .
Combining (17)-(18), we get u · Au ≤ chd−2 kuk22 . It proves the corollary.

(18)


5.6. Error estimates. In this section, we get a priori error estimates for the aggregated FE scheme (4).
In order to do that, we prove first approximability properties of the corresponding spaces.
Lemma 5.10. Let us consider an aggregated FE space of order q, m ≤ q, 1 ≤ s ≤ m ≤ q + 1, 1 ≤ p ≤ ∞,
and m > dp . Given a function u ∈ H m (Ω), it holds:
inf

uh ∈Vhagg

ku − uh kWps (Ω) . hm−s |u|Wpm (Ω) .

Proof. Under the assumptions of the lemma, we have that the following embedding Wpm (Ω) ⊂ C 0 (Ω̄) is
continuous (see, e.g., [28, p. 486]). Thus, given a function u ∈ Wpm (Ω), let represent with σ(u) the
vector ofP
nodal values in Ωin , i.e., σ(u)a = uh (xa ) for a ∈ N (Ωin ). We define the interpolation operator
.
Ih (u) = a∈N (Ω) φa [Eσ(u)]a .
Given a cut cell K ∈ Thact , the fact that its DOFs values only depend on interior DOFs in Ω̄K , and since
m (K) ⊆ W m (K), it follows from the upper bound of the norm of the
each shape function φa belongs to W∞
p
nodal extension operator in Lemma 5.1 that kIh (u)kWpm (K) ≤ CkukC 0 (Ω̄K ) (see also [22, Lemma 4.4.1]).
On the other hand, we consider an arbitrary function π(u) ∈ Wpm (Ω) such that π(u)|K ∈ Pq (ΩK ). We
note that, by construction, π(u)|K = Ih (π(u))|K . Thus, we have:
ku − Ih (u)kWpm (K) ≤ ku − π(u)kWpm (K) + kIh (π(u) − u)kWpm (K)
. ku − π(u)kWpm (K) + kπ(u) − ukC 0 (Ω̄K ) . ku − π(u)kWpm (ΩK ) ,
where we have used in the last inequality the previous continuous embedding. Since ΩK is an open bounded
domain with Lipschitz boundary by definition, one can use the Deny-Lions lemma (see, e.g., [28]). As a
result, the π(u) that minimizes the right-hand side holds:
ku − Ih (u)kWpm (ΩK ) . |u|Wpm (ΩK ) .


Using standard scaling arguments, we prove the lemma.

Theorem 5.11. If Ω has smoothing properties and the solution u of the continuous problem (2) belongs
to Wpm (Ω) for dp < m ≤ q, the solution uh ∈ Vhagg of (4) satisfies the following a priori error estimate:
ku − uh kH 1 (Ω) ≤ hm−1 |u|H m (Ω) .
Proof. Combining the consistency of the numerical method, i.e., A(u, vh ) = `(vh ), and the continuity and
coercivity of the bilinear form in Th. 5.7, we readily get, using standard FE analysis arguments:
|||wh − uh |||2h . A(wh − uh , wh − uh ) = A(wh − u, wh − uh ) . |||wh − u|||h |||wh − uh |||h ,
for any wh ∈ Vhagg . On the other hand, we use the trace inequality (see [29])
kψk2L2 (∂T ) . |∂T |−1 kψk2L2 (T ) + |∂T |kψk2H 1 (T ) ,

for any ψ ∈ H 1 (T ).

Using this trace inequality, we get:
−1
−2
2
2
2
2
h−1
K ku − wh kL2 (ΓD ∩K) ≤ hK ku − wh kL2 (∂ΩK ) . hK ku − wh kL2 (ΩK ) + ku − wh kH 1 (ΩK ) .

Combining the previous bound with the approximability property in Lemma 5.10, we readily get
|||wh − u|||h . hm−1 |u − wh |H m (Ω) .
It proves the theorem.



AGGREGATED UNFITTED FINITE ELEMENT METHOD

14

6. Numerical experiments
6.1. Setup. The numerical examples below consider as a model problem the Poisson equation with nonhomogeneous Dirichlet boundary conditions. The value of the source term and the Dirichlet function are
defined such that the PDE has the following manufactured exact solution:
 
1/2 
2
2
2
,
u(x, y, z) = sin 4π (x − 2.3) + y + z
(x, y) ∈ Ω ⊂ R2 , z = 0 in 2D,

(x, y, z) ∈ Ω ⊂ R3 in 3D.

We consider two different geometries, a 2D circle and a 3D complex domain with the shape of a popcorn
flake (see Fig. 6). These geometries are often used in the literature to study the performance of unfitted
FE methods (see, e.g., [4], where the definition of the popcorn flake is found). In all cases, we use the
cuboid [0, 1]d , d = 2, 3 as the bounding box on top of which the background Cartesian grid is created. For
the sake of illustration, Fig. 6 displays both the considered geometries, numerical solution and bounding
box.

(a) 2D case

(b) 3D case

Figure 6. View of the problem geometries, numerical solution and bounding box.
The main goal of the following tests is to evaluate the (positive) effect of using the aggregation-based FE
space Vhagg instead of the usual one Vhact . In the next plots, the results for the usual (un-aggregated) FE
space are labeled as standard, whereas the cases with the aggregation are labeled as aggregated (or aggr.
in its short form). In all the examples, we use Lagrangian reference FEs with bi-linear and bi-quadratic
shape functions in 2D, and tri-linear and tri-quadratic ones in 3D.
Both the standard and the aggregated formulations have been implemented in the object-oriented HPC
code FEMPAR [17]. The system of linear equations resulting from the problem discretization are solved
within FEMPAR with a sparse direct solver from the MKL PARDISO package [30]. Condition number
estimates are computed outside FEMPAR using the MATLAB function condest.2 For the standard
unfitted FE space Vhact , we expect very high condition numbers that can hinder the solution of the discrete
system using standard double precision arithmetic. To address this effect and avoid the breakdown of sparse
direct solvers, we bound from below the minimum distance between the mesh nodes and the intersection
of edges with the boundary Γ to a small numerical threshold Dmin proportional to the cell size, namely
Dmin = εh, where ε is a (mesh independent) user defined tolerance. If the edge cut-node distance is
below this threshold , the edge cut is collapsed with the node, perturbing the geometry. In the numerical
experiments, we take ε = 10−6 and ε = 10−3 in 2D and 3D respectively. Using the fact that ηK ∼ εd , we can
rewrite the condition number estimate (9) in terms of the user-defined tolerance ε as κ(A) ∼ ε−d(2q+1−2/d) .
2MATLAB is a trademark of THE MATHWORKS INC.

AGGREGATED UNFITTED FINITE ELEMENT METHOD

15

For instance, we have κ(A) ∼ ε−7 and κ(A) ∼ ε−13 for first and second order interpolations, respectively,
in 3D. This illustrates that the condition numbers expected for second order interpolation are extremely
high as it is confirmed below unless very large values of ε are considered. However, the value of ε cannot be
increased without affecting the numerical error, since it perturbs the geometry, and destroys at some point
the order of convergence of the numerical method. Similar perturbation-based techniques with analogous
problems have been used in the frame of the finite cell method in [8]. We note that the tolerance ε is not
needed at all when using the aggregated FE space.
6.2. Moving domain experiment. In the first numerical experiment, we study the robustness of the
unfitted FE formulations with respect to the relative position between the unfitted boundary and the
background mesh. To this end, we consider two moving domains that can travel along one of the diagonals
of the bounding box (see Fig. 7). The considered geometries are obtained by scaling down the circle
and the popcorn flake depicted in Fig. 6 by a factor of 0.25. In both cases, the position of the bodies is
controlled by the value of the parameter ` (i.e., the distance between the center of the body and a selected
vertex of the box). As the value of ` varies, the objects move and their relative position with respect to
the background mesh changes. In this process, arbitrary small cut cells can show up, leading to potential
conditioning problems. In this experiment, we consider a background mesh with element size h = 2−5 .

(a) 2D case

(b) 3D case

Figure 7. Setup of the moving domain experiment.
Fig. 8 shows the condition number estimate of the underlying linear systems varying the position of
the physical domain Ω. The plot is generated using a sample of 200 different values of `. It is observed
that the condition numbers are very sensitive to the position of the domain for the standard unfitted
FE formulation, whereas the condition numbers are nearly independent of the position when using the
aggregation-based FE spaces. Note that the standard formulation leads to very high condition numbers
specially for second order interpolations and the 3D case. Moving from 1st order to 2nd order leads to a
rise in the condition number between 10 and 15 orders of magnitude. The same disastrous effect is observed
when moving from 2D to 3D. In contrast, the condition number is nearly insensitive to the number of
space dimensions, and mildly depends on the interpolation order (as for body-fitted methods) when using
aggregation-based FE spaces.
From the results shown in Fig. 8, it is clear that the aggregation-based FE spaces are able to dramatically
improve the condition numbers associated with the standard unfitted FE formulation. The next question
is how cell aggregation impacts on the accuracy of the numerical solution. In order to quantify this effect,
Fig. 9 shows the computed energy norm of the discretization error. It is observed that the error is slightly
increased when using the aggregation-based FE spaces. This is because the considered meshes in this
moving domain experiment are rather coarse. The error increments become negligible for finer meshes

AGGREGATED UNFITTED FINITE ELEMENT METHOD

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.

log10(condest(A))

25
20
15
10

30
p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.

25

log10(condest(A))

30

16

5

20
15
10
5

0

0
0.3

0.4

0.5

0.6

0.7

0.3

0.4

0.5

ℓ

0.6

0.7

ℓ

(a) 2D case

(b) 3D case

Figure 8. Condition number vs. domain position.
(see Section 6.3 below). In this example, we cannot compute a solution for all the values of ` for 3D and
2nd order interpolation without using cell aggregation (see the discontinuous fine red curve in Fig. 9b).
The condition numbers are so high (order 1030 ) that the system is intractable, even with a sparse direct
solver, using standard double precision floating point arithmetic.
1

1

log10(Error energy norm)

0.5
0
-0.5
-1
-1.5
-2
-2.5
-3
0.3

0.4

0.5

0.6

0.7

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.

0.5

log10(Error energy norm)

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.

0
-0.5
-1
-1.5
-2
-2.5
-3
0.3

0.4

0.5

ℓ

(a) 2D case

0.6

0.7

ℓ
(b) 3D case

Figure 9. Error energy norm vs. domain position.

6.3. Convergence test. The second experiment is devoted to study the asymptotic behavior of the
methods as the mesh is refined. To this end, we consider the geometries and bounding boxes displayed in
Fig. 6, which are discretized with uniform Cartesian meshes with element sizes h = 2−m , m = 3, 4, . . . , 9
in 2D, and m = 3, 4, 5, 6 in 3D.
First, we study how the size of the aggregates scales when the mesh is refined. Fig. 10 shows that the
aggregate size is 2h in 2D, whereas it tends to 3h in the 3D case. These results agrees with the theoretical
bounds for the aggregate size discussed in Section 2.
Then, we study the scaling of the condition numbers with respect to the mesh size (see Fig. 11). For the
aggregation-based FE spaces, the condition numbers of the stiffness matrix scales as h−2 , like in standard
FE methods for body fitted meshes. This confirms the theoretical result of Corollary 5.9. Conversely,
the condition number has an erratic behavior if cell aggregation is not considered. The reason is that, as
shown in the previous experiment (cf. Section 6.2), the standard unfitted FE formulation leads to condition
numbers very sensitive to the position of the unfitted boundary. Several configurations of cut cells can show
up when the mesh is refined, leading to very different condition numbers. As in the previous experiment,
the condition number is very sensitive to the interpolation order and number of space dimensions for the
standard unfitted FE formulations. This effect is reverted when using cell aggregates.

AGGREGATED UNFITTED FINITE ELEMENT METHOD

17

8

H/h

6
4
2D
3D
Bound in 2D
Bound in 3D

2
0
-3.5

-3

-2.5

-2

-1.5

-1

log10(h)

Figure 10. Scaling of the ratio maximum aggregate size (H) / element size (h).
40

40
p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.
slope -2

35

30

30

25

25

log10(condest(A))

log10(condest(A))

35

20

15

20

15

10

10

5

5

0

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.
slope -2

0
-2.5

-2

-1.5

-1

log10(h)

-2.5

-2

-1.5

-1

log10(h)

(a) 2D case

(b) 3D case

Figure 11. Scaling of the condition number upon mesh refinement.
Finally, we study the convergence of the discretization error. To this end, Figs. 12 and 13 report the
discretization errors measured both in the energy norm and in the L2 norm. Like in the previous experiment
(cf. Section 6.2), the discrete system could not be solved when using the finest meshes in 3D for 2nd order
interpolation without using cell-aggregation due to extremely large condition numbers (see the incomplete
curve in Fig. 12b). The results show that the error increment associated with the aggregation-based FE
space becomes negligible when the mesh is refined. Moreover, the theoretical results of Section 5.6 are
confirmed: optimal order of convergence is (asymptotically) achieved in all cases when using aggregationbased FE spaces both for the energy and the L2 norm, for 1st and 2nd order interpolations, and for 2 and
3 spatial dimensions.
7. Conclusions
We have proposed a novel technique to construct FE spaces designed to improve the conditioning
problems associated with unfitted FE methods. The spaces are defined using cell aggregates obtained
by merging the cut cells to interior cells. In contrast to related methods in the literature, the proposed
technique is easy to implement in existing FE codes (it only involves cell-wise constraints) and it is general
enough to deal with both continuous and DG formulations. Another novelty with respect to previous works
is that we include the mathematical analysis of the method. For elliptic problems, we have proved that 1)
the novel FE space leads to condition numbers that are independent from small cut cells, 2) the condition

1

1

0

0

-1

-1

-2

-2

log10(Error energy norm)

log10(Error energy norm)

AGGREGATED UNFITTED FINITE ELEMENT METHOD

-3

-4

-5

-6

-3

-4

-5

-6

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.
slope 1
slope 2

-7

18

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.
slope 1
slope 2

-7

-8

-8
-2.5

-2

-1.5

-1

-2.5

log10(h)

-2

-1.5

-1

log10(h)

(a) 2D case

(b) 3D case

1

1

0

0

-1

-1

-2

-2

log10(Error L2 norm)

log10(Error L2 norm)

Figure 12. Convergence of the discretization error in energy norm.

-3

-4

-5

-6

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.
slope 2
slope 3

-7

-8

-3

-4

-5

-6

p=1, standard
p=2, standard
p=1, aggr.
p=2, aggr.
slope 2
slope 3

-7

-8
-2.5

-2

-1.5

log10(h)

(a) 2D case

-1

-2.5

-2

-1.5

-1

log10(h)

(b) 3D case

Figure 13. Convergence of the discretization error in L2 norm.

number of the resulting system matrix scales with the inverse of the square of the size of the background
mesh as in standard FE methods, 3) the penalty parameter of Nitsche’s method is bounded from above,
and 4) the optimal FE convergence order is recovered. These theoretical results are confirmed with 2D
and 3D numerical experiments using both first and second order interpolations.

AGGREGATED UNFITTED FINITE ELEMENT METHOD

19

References
[1] S. Badia, F. Nobile, and C. Vergara.
Fluid–structure partitioned procedures based on
Robin transmission conditions. Journal of Computational Physics, 227(14):7027–7051, 2008.
doi:10.1016/j.jcp.2008.04.006.
[2] D. Kamensky, M.-C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, and
T. J. R. Hughes. An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves. Computer Methods in Applied Mechanics and Engineering, 284:
1005–1053, 2015. doi:10.1016/j.cma.2014.10.040.
[3] T. Belytschko, N. Moës, S. Usui, and C. Parimi. Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering, 50(4):993–1013, 2001. doi:10.1002/10970207(20010210)50:4<993::AID-NME164>3.0.CO;2-M.
[4] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing. CutFEM: Discretizing geometry
and partial differential equations. International Journal for Numerical Methods in Engineering, 104
(7):472–501, 2015. doi:10.1002/nme.4823.
[5] R. Mittal and G. Iaccarino. Immersed Boundary Methods. Annual Review of Fluid Mechanics, 37
(1):239–261, 2005. doi:10.1146/annurev.fluid.37.061903.175743.
[6] D. Schillinger and M. Ruess. The Finite Cell Method: A Review in the Context of Higher-Order
Structural Analysis of CAD and Image-Based Geometric Models. Archives of Computational Methods
in Engineering, 22(3):391–455, 2014. doi:10.1007/s11831-014-9115-y.
[7] Y. Sudhakar, J. P. Moitinho de Almeida, and W. A. Wall. An accurate, robust, and easy-to-implement
method for integration over arbitrary polyhedra: Application to embedded interface methods. Journal
of Computational Physics, 273:393–415, 2014. doi:10.1016/j.jcp.2014.05.019.
[8] J. Parvizian, A. Düster, and E. Rank. Finite cell method: H- and p-extension for embedded domain
problems in solid mechanics. Computational Mechanics, 41(1):121–133, 2007. doi:10.1007/s00466007-0173-y.
[9] Y. Saad. Iterative Methods for Sparse Linear Systems. Other Titles in Applied Mathematics. Society
for Industrial and Applied Mathematics, 2003.
[10] W. Briggs, V. Henson, and S. McCormick. A Multigrid Tutorial, Second Edition. Other
Titles in Applied Mathematics. Society for Industrial and Applied Mathematics, 2000.
doi:10.1137/1.9780898719505.
[11] S. Badia, A. F. Martín, and J. Principe. Multilevel Balancing Domain Decomposition at Extreme
Scales. SIAM Journal on Scientific Computing, pages C22–C52, 2016. doi:10.1137/15M1013511.
[12] A. Menk and S. P. A. Bordas. A robust preconditioning technique for the extended finite element method. International Journal for Numerical Methods in Engineering, 85(13):1609–1632, 2011.
doi:10.1002/nme.3032.
[13] L. Berger-Vergiat, H. Waisman, B. Hiriyur, R. Tuminaro, and D. Keyes. Inexact Schwarz-algebraic
multigrid preconditioners for crack problems modeled by extended finite element methods. International Journal for Numerical Methods in Engineering, 90(3):311–328, 2012. doi:10.1002/nme.3318.
[14] B. Hiriyur, R. Tuminaro, H. Waisman, E. Boman, and D. Keyes. A Quasi-algebraic Multigrid Approach to Fracture Problems Based on Extended Finite Elements. SIAM Journal on Scientific Computing, 34(2):A603–A626, 2012. doi:10.1137/110819913.
[15] F. de Prenter, C. V. Verhoosel, G. J. van Zwieten, and E. H. van Brummelen. Condition number
analysis and preconditioning of the finite cell method. Computer Methods in Applied Mechanics and
Engineering, 316:297–327, 2017. doi:10.1016/j.cma.2016.07.006.
[16] S. Badia and F. Verdugo. Robust and scalable domain decomposition solvers for unfitted finite element
methods. Computational and Applied Mathematics, In press, 2017.
[17] S. Badia, A. F. Martín, and J. Principe. FEMPAR: An object-oriented parallel finite element framework. Archives of Computational Methods in Engineering, In press, 2017.
[18] E. Burman and P. Hansbo.
Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method. Applied Numerical Mathematics, 62(4):328–341, 2012.
doi:10.1016/j.apnum.2011.01.008.
[19] C. Helzel, M. Berger, and R. Leveque. A High-Resolution Rotated Grid Method for Conservation
Laws with Embedded Geometries. SIAM Journal on Scientific Computing, 26(3):785–809, 2005.

AGGREGATED UNFITTED FINITE ELEMENT METHOD

20

doi:10.1137/S106482750343028X.
[20] F. Kummer. Extended Discontinuous Galerkin methods for multiphase flows: The spatial discretization. Center for Turbulence Research Annual Research Briefs, pages 319–333, 2013.
[21] O. Marco, R. Sevilla, Y. Zhang, J. J. Ródenas, and M. Tur. Exact 3D boundary representation in
finite element analysis based on Cartesian grids independent of the geometry. International Journal
for Numerical Methods in Engineering, 103(6):445–468, 2015. doi:10.1002/nme.4914.
[22] S. C. Brenner and R. Scott. The Mathematical Theory of Finite Element Methods. Springer, softcover
reprint of hardcover 3rd ed. 2008 edition, 2010.
[23] R. Becker. Mesh adaptation for Dirichlet flow control via Nitsche’s method. Communications in
Numerical Methods in Engineering, 18(9):669–680, 2002. doi:10.1002/cnm.529.
[24] J. Nitsche. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von
Teilräumen, die keinen Randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen
Seminar der Universität Hamburg, 36(1):9–15, 1971. doi:10.1007/BF02995904.
[25] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of discontinuous Galerkin
methods for elliptic problems. SIAM Journal on Numerical Analysis, 39(5):1749–1779, 2002.
[26] S. Agmon and N. S. F. (U.S.). Lectures on elliptic boundary value problems. Van Nostrand mathematical studies. Van Nostrand, 1965.
[27] H. C. Elman, D. J. Silvester, and A. J. Wathen. Finite elements and fast iterative solvers: with
applications in incompressible fluid dynamics. Oxford University Press, 2005.
[28] A. Ern and J.-L. Guermond. Theory and Practice of Finite Elements. Springer, 2004.
[29] D. Arnold. An Interior Penalty Finite Element Method with Discontinuous Elements. SIAM Journal
on Numerical Analysis, 19(4):742–760, 1982. doi:10.1137/0719052.
[30] Intel MKL PARDISO - Parallel Direct Sparse Solver Interface. https://software.intel.com/
en-us/articles/intel-mkl-pardiso.
Department of Civil and Environmental Engineering.
Girona 1-3, Edifici C1, 08034 Barcelona, Spain.

Universitat Politècnica de Catalunya, Jordi

CIMNE – Centre Internacional de Mètodes Numèrics en Enginyeria, Parc Mediterrani de la Tecnologia,
UPC, Esteve Terradas 5, 08860 Castelldefels, Spain.

