1

The Unfitted Discontinuous Galerkin Method for
Solving the EEG Forward Problem

arXiv:1601.07810v2 [] 2 May 2016

Andreas Nüßing, Carsten H. Wolters, Heinrich Brinck, Christian Engwer

Abstract— Objective: The purpose of this study is to introduce
and evaluate the unfitted discontinuous Galerkin finite element
method (UDG-FEM) for solving the electroencephalography
(EEG) forward problem. Methods: This new approach for source
analysis does not use a geometry conforming volume triangulation, but instead uses a structured mesh that does not resolve the
geometry. The geometry is described using level set functions and
is incorporated implicitly in its mathematical formulation. As no
triangulation is necessary, the complexity of a simulation pipeline
and the need for manual interaction for patient specific simulations can be reduced and is comparable with that of the FEM for
hexahedral meshes. In addition, it maintains conservation laws
on a discrete level. Here, we present the theory for UDG-FEM
forward modeling, its verification using quasi-analytical solutions
in multi-layer sphere models and an evaluation in a comparison
with a discontinuous Galerkin (DG-FEM) method on hexahedral
and on conforming tetrahedral meshes. We furthermore apply
the UDG-FEM forward approach in a realistic head model
simulation study. Results: The given results show convergence
and indicate a good overall accuracy of the UDG-FEM approach.
UDG-FEM performs comparable or even better than DG-FEM
on a conforming tetrahedral mesh while providing a less complex
simulation pipeline. When compared to DG-FEM on hexahedral
meshes, an overall better accuracy is achieved. Conclusion: The
UDG-FEM approach is an accurate, flexible and promising
method to solve the EEG forward problem. Significance: This
study shows the first application of the UDG-FEM approach to
the EEG forward problem.

I. I NTRODUCTION
The localization of current sources in the human brain from
surface electroencephalography (EEG) measurements (the inverse problem) requires a model for the forward problem, i.e.,
the determination of surface potentials from current sources
in the cortical sheet of the human brain [1], [2]. Several
different approaches have been proposed to solve the EEG
forward problem. When approximating the head by a multicompartment sphere model, quasi-analytical series expansion
formulas are available [3]. More realistic methods are based on
head volume conductor representations obtained from quasinoninvasive magnetic resonance imaging (MRI). Under these
Andreas Nüßing is with the Institute for Bioinformatics and Chemoinformatics, Westfälische Hochschule, Recklinghausen, Germany, with the Institute
for Computational and Applied Mathematics, University of Münster, Germany
and with the Institute for Biomagnetism and Biosignalanalysis, University of
Münster, Germany.
C. H. Wolters is with the Institute for Biomagnetism and Biosignalanalysis,
University of Münster, Germany.
H. Brinck is with the Institute for Bioinformatics and Chemoinformatics,
Westfälische Hochschule, Recklinghausen, Germany.
Christian Engwer is with the Institute for Computational and Applied
Mathematics and the Cluster of Excellence EXC 1003, Cells in Motion, CiM,
University of Münster, Germany.

approaches are boundary element methods (BEM) [4]–[7],
finite volume methods (FVM) [8], finite difference methods
(FDM) [9]–[11] and -finite element methods (FEM) [12]–
[20]. The method presented in this paper is closely related
to FEM approaches. The latter were shown to produce highly
accurate solutions [19]. In combination with transfer matrices
and efficient linear solvers [2], [12], [14], the computational
effort could be reduced. By using volumetric meshes, they are
able to handle anisotropic conductivities and complex model
geometries [21].
A common FEM approach is to use constrained Delaunay
tetrahedralization (CDT) to generate a conforming tetrahedral
mesh which is constructed from non-intersecting surface triangulations [21]. Although the method allows smooth tissue
surface representations, the process of constructing proper
non-intersecting surfaces might become an involved task. In
order to allow CDT meshing, unrealistic model features might
be introduced, such as the artificial closing of holes (e.g.
optical canals or foramen magnum). Furthermore, the surfaces
are required to be nested, while they are touching in a realistic
scenario (e.g. the inner skull surface and the brain surface).
Such limitations can be circumvented by using hexahedral
models. They can be directly generated from voxel-based MRI
images and are used in several source analysis applications
[13], [22], [23]. The hexahedral approach keeps the advantages
of the FEM mentioned above, while reducing the effort for
creating individual head models.
However, one problem of FEM modeling using regular
hexahedral meshes is the possibility of the occurrence of skull
leakages: In areas where the skull is very thin, for example
the temporal bone, where skull thickness is 2 mm or even less
[24, Table 2], physically unsound leakages of current through
the skull might occur in regular hexahedral approaches with
insufficient resolution [25], [26]. It was shown in [27] that
an appropriate skull modeling is of special importance for
accurate EEG forward modeling. The effects of skull leakages
might thus lead to especially significant errors when hexahedral FEM modeling is used in combination with insufficient
resolution. A way to alleviate this problem has been proposed
in [26], where instead of the node based approach of the
Lagrange or Continuous Galerkin FEM (CG-FEM), a cell
based Discontinuous Galerkin FEM (DG-FEM) method was
introduced. In this paper, we will therefore focus on DG-FEM
approaches.
A second problem of regular hexahedral FEM approaches
is the stair-case like representation of smooth head tissue
surfaces, producing geometrical inaccuracies that lead to modeling errors, especially in combination with lower resolutions.

2

Therefore, methods have been developed to reduce such geometrical errors. In [28], the structured mesh was modified
to obtain a smoother interface between different tissue compartments. However, this method might still suffer from skull
leakages and geometry adaptation is limited by the need for
positive Jacobian determinants. In [17], a so-called immersed
FEM has been proposed for the EEG forward problem. It
combines the simplicity of a hexahedral CG-FEM approach
with the accuracy of a conforming tetrahedral approach with
regard to the modeling of smooth tissue surfaces. Using a level
set function to represent the surface, the local basis functions
of the hexahedral approach were modified to conform to the
tissue compartments. It was shown that the accuracy of a
conforming tetrahedral approach could be met, while offering
a simpler simulation pipeline.
In this paper, we present the unfitted discontinuous Galerkin
finite element method (UDG-FEM) to solve the EEG forward
problem. The method was first introduced in the context of
micro-scale simulations in porous media [29] and has since
then found different applications [30]–[32]. It combines the
advantages of the DG-FEM approach on a structured hexahedral mesh with implicit surface representations given by level
set functions and thereby the avoidance of skull leakages and
inappropriate modeling of smooth tissue surfaces. As we will
show, it outperforms a state-of-the-art DG-FEM approach and
leads to simpler forward modeling pipelines because it is not
restricted to nested compartments and can handle intersecting
surfaces such as touching surfaces of inner skull and brain.
Our paper is structured as follows: After a thorough description of the theory in section II, we provide methodological
aspects of our verification and evaluation in section III. Results
are then presented in section IV and discussed in section V.
The paper ends with concluding remarks in chapter VI.
II. T HEORY
A. A discontinuous Galerkin (DG-FEM) method for solving
the EEG forward problem
The EEG forward problem can be solved by providing a
solution to Poisson’s equation which results from the quasistatic Maxwell’s equations [1], [2]. Let Ω ⊂ R3 denote the
head domain and ∂Ω its surface. The task is then to find the
electric scalar potential u for which
∇ · σ∇u = f

in Ω

(1)

σ∇u · n = 0

in ∂Ω

(2)

holds. σ : Ω → Sym3 (R) denotes the tissue conductivity
tensor of second rank and f the source current density. For
each x ∈ Ω, σ(x) is assumed to be symmetric and positive
definite, but it is allowed to be discontinuous over tissue
boundaries. A common source model, which we will also use
here, is the mathematical dipole. It is represented by a position
x0 ∈ Ω and a moment vector M ∈ R3 and given as:
f (x) = M · ∇δ(x − x0 )

(3)

where δ denotes Dirac’s delta distribution. Note that the
mathematical dipole model is not a function and that its
divergence has to be considered in distributional sense. Note

also that for (1) to hold in a strong sense, more rigorous
regularity assumptions on u, σ and f would be required, so
that (1) should be seen in a symbolic sense.
In [26], a discontinuous Galerkin (DG-FEM) formulation
for (1) has been derived, which we will briefly recall here.
For a detailed description and an analysis of the properties
of DG-FEM, we refer to [26]. Note that this formulation also
serves as the foundation of the unfitted discontinuous Galerkin
method (UDG-FEM) presented below. First, we introduce a
volume triangulation of Ω:
Th = {Ei |i ∈ I = {0, . . . , N − 1}}
[
Ei = Ω
Ei ∩ Ej = ∅ ∀i 6= j,

(4)
(5)

i∈I

of open sets Ei ⊂ R3 . For such an open set E, E denotes its
closure, i.e. E = E ∪ ∂E. The mesh width h ∈ R is defined
as h := max{diam(E) : E ∈ Th }. In the following, we will
restrict the description to the case that all elements are either
tetrahedrons or hexahedrons. The skeleton Γh of Th is defined
as Γh := {γi,j = E i ∩ E j : Ei , Ej ∈ Th , i 6= j, |γi,j | > 0}.
Let Vhk = {u ∈ L2 (Ω) : u|E ∈ P k (E) ∀E ∈ Th } denote the
broken polynomial space on Th . P k (E) denotes a space of
polynomials on E of degree k ∈ N. In the following, we will
assume that σ is constant on each Ei and denote its value by
σi . The discontinuous Galerkin method for solving (1), which
we will use, then reads: Find uh ∈ Vhk such that for all test
functions vh ∈ Vhk
a(uh , vh ) + J(uh , vh ) = l(vh )
holds. The bilinear forms a and J are given as
Z
Z
a(uh , vh ) =
σ∇uh · ∇vh dx −
Juh K · hσ∇vh ids
Ω
Γh
Z
−
Jvh K · hσ∇uh ids
Z Γh
τγ
J(uh , vh ) =η
Juh K · Jvh Kds
h
Z Γh γ
l(vh ) =
f vh dx

(6)

(7)
(8)
(9)
(10)

Ω

The jump Juh K∈ R3 on the intersection between two elements
Ei and Ej with unit outer normals ni ∈ R3 and nj ∈ R3 ,
respectively, is defined as Juh K := uh |Ei ni + uh |Ej nj . The
weighted average of the flux of uh on the interface is denoted
by hσ∇uh i ∈ R3 . With δi := nti σi ni and δj := ntj σj nj , this
can be defined as
δj
δi
σi ∇uh |Ei +
σj ∇uh |Ej
(11)
hσ∇uh i :=
δi + δj
δi + δj
Note that both the jump of the potential and the weighted
average of the flux are vector valued quantities and that the
normals ni and nj are opposing vectors, i.e. ni = −nj . The
factor τγ scales the penalty term J at conductivity jumps on
an edge γ = γi,j . It is defined as the harmonic average of δi
and δj , i.e. τγ := 2δi δj /(δi + δj ). Note that the bilinear form
a+J is symmetric. This variant of a DG-FEM approach is also
called symmetric weighted interior penalty Galerkin (SWIPG)
method [33]. For a sufficiently large η > 0, the problem has

3

a unique solution. Note that each integral can be replaced by
a sum over element local contributions.
In the following, we will use P1 (E) on tetrahedral and
Q1 (E) on hexahedral elements for the local polynomial spaces
in Vhk . The polynomial space of each element Ej ∈ Th is
spanned by Nb ∈ N local basis functions ϕj0 , . . . , ϕjNb −1
(Nb = 4 for tetrahedrons, Nb = 8 for hexahedrons). On the
reference element, these local basis functions can be given
as 1, x, y, z for tetrahedrons and 1, x, y, z, xy, xz, yz, xyz for
hexahedrons, respectively. By mapping each local basis function on each element to a unique global index, we obtain the
global basis functions ϕ0 , . . . , ϕn−1 , where n = Nb · N . Note
that the support of each global basis function spans only a
single element of the mesh.
B. An unfitted discontinuous Galerkin (UDG-FEM) method
for solving the EEG forward problem
The former version of the DG-FEM approach uses a triangulation, which resolves the geometry. The unfitted discontinuous
Galerkin (UDG-FEM) method takes a different approach and
describes the geometry using level set functions. We assume
that the head domain is embedded in a larger domain Ω̂ ⊂ R3 .
A level set function for a subdomain Ω ⊂ Ω̂ is a scalar
continuous function Φ : Ω̂ → R with the property


< 0 if x ∈ Ω,
(12)
Φ(x) = 0 if x ∈ ∂Ω,


> 0 if x ∈ Ω̂ \ Ω.
As an example, a level set function for the unit sphere
can be defined as Φ(x) = kxk − 1. On a structured grid
Th , the level set function is approximated as a piecewise
multilinear Q1 function Φh by evaluating Φ at each grid
node. In the following, we will call this structured mesh
fundamental mesh. By employing multiple level set functions
Φ0h , . . . , ΦhL−1 , L ∈ N, we can differentiate between multiple
domains. A level set can represent the boundary between two
tissue compartments or delimit the domain of a single compartment (see e.g. [34]). Each such level set function separates
Ω̂ into two parts with respect to its sign. We will denote these
parts by Ωi,− and Ωi,+ for the negative and positive side
of the level set function i. From these parts, we can create
domains D0 , . . . , DD−1 , D ∈ N, consisting of intersections
of negative
T and positive sides of the level set functions, i.e.
Dj ⊂ { (i,p)∈I Ωi,p : I ⊂ {0, . . . , L − 1} × {−, +}} for
j ∈ {0, . . . , D − 1}. Each domain can represent a tissue
compartment. For a domain Dj , we define its support
S Ω(Dj )
as the union of these intersections, i.e. Ω(Dj ) := Ω̃∈Dj Ω̃.
We require that the supports of all domains are pairwise
disjoint, i.e., Ω(Di ) ∩ Ω(Dj ) = ∅ holds for i 6= j. For
each element Ei , we set its set of intersecting domains to
D(Ei ) := {Dj : Ei ∩ Ω(Dj ) 6= ∅}. From this domain
information, we generate a cut cell triangulation:
[
Th=
{Ei ∩ Ω(Dj ) : Dj ∈ D(Ei )}
(13)
|
{z
}
E ∈T
i

h

=:Eij

The elements Eij of T h are called cut cells. Eij is thus defined
as the part of element Ei belonging to domain Dj . In the

Fig. 1: Construction of the subtriangulation for a single
bilinear level set function. In the left image, two domains
are delimited by a bilinear level set function and in the right
image, the resulting discrete domains are shown. The dashed
lines delimit the elements of the subtriangulation which are
only used for integration.

following, we will require that the conductivity tensor σ is
constant on each cut cell. On T h , we can again define a broken
polynomial space and pose the DG problem (6). The local
polynomial spaces are now defined on each cut cell. They
can be obtained by simply restricting the local basis function
of the element to the cut cell, i.e. setting their value to zero
outside of the cut cells domain. The support of the global
basis functions thus spans only a single cut cell. The integrals
in (6) thus reduce to integrals of the local basis functions over
cut cells. Note that an element of the fundamental mesh can
contain multiple sets of local basis functions, depending on the
number of different domains in that element. An example can
be seen in the left part of Fig. 1, where, restricted to a single
element of the fundamental mesh, two domains are separated
by a bilinear level set function, thus forming two cut cells.
For each of these two cut cells, a set of local basis functions
is introduced.
For the evaluation of integrals over a cut cell and its
boundary, several methods have been proposed [35]. We use
an extended marching cubes algorithm [36]. The domain of
a cut cell is approximated by a first order subtriangulation
into simple elements. As the discrete level set functions Φih
are Q1 functions, they are completely defined by their values
at the mesh nodes. The subtriangulation can be computed
automatically solely based on these values. The computation
uses precomputed look up tables and produces topologically
correct triangulations. An example for such a subtriangulation
in 2D can be seen in the right part of Fig. 1 (the subtriangulation of the green subdomain is marked by the dashed lines).
For a detailed description of the extended marching cubes
algorithm, we refer to [36]. When an element is cut by multiple
level set functions, we apply the marching cubes algorithm
recursively on the elements of the subtriangulation. Note that
we introduce a slight error in the interface reconstruction
by this recursive application. However, this error decreases
with decreasing h and should be negligible when using linear
basis functions. The integration over a cut cell E can be
replaced by integrations over the simple elements E j of the

4

subtriangulation:
Z
XZ
g(x)dx =
E

j

=

XZ
j

g(x)dx

(14)



g µj (x̂) det Jµj (x̂) dx̂

(15)

Ej

Ê j

with the affine map µj : Ê j 7→ E j mapping local coordinates
of the reference element Ê j to global coordinates of the
element of the subtriangulation E j . The local integrals over
the reference elements of the subelements can be carried out
using common quadrature rules. Besides the subtriangulation
of the volume, the marching cubes algorithm also produces a
subtriangulation of the interface. This can be used to compute
the corresponding skeleton integrals.
When inserting the global basis functions into the bilinear
form and the right hand side (6), we get a sparse linear
equation system M y = f , with M ∈ Rn×n and y, f ∈ Rn
where n = Nb · N denotes the global number of degrees of
freedom. M = (mij ) and f = (fij ) are defined as
η
mij = a(ϕi , ϕj ) + J(ϕi , ϕj ), fi = f (ϕi )
(16)
h
Since all basis functions have support in at most one element
or cut cell, the entries of M and f can be blocked element
or cut cell wise. After solvingP
this linear system for y, the
potential u is given as u(x) = i yi ϕi (x).
III. M ETHODS
A. Source model
For standard CG-FEM, several different source models have
been proposed, such as the partial integration [12], [17], [19],
the Saint-Venant [2], [20], the Whitney or Raviart Thomas
[18], [37] or the subtraction approach [13], [38]. In principle,
the presented method is applicable to any given source model.
Due to its simplicity and accuracy when compared to the subtraction approach in DG-FEM simulations [26], [39], we use
here the partial integration approach to model a mathematical
dipole. When tested with a basis function ϕi , the following
right-hand side results:
(
−M · ∇ϕi (x0 ) if x0 ∈ support(ϕi )
fi = f (ϕi ) =
(17)
0
else
For the standard DG-FEM discretization, the support of the
discrete source model is completely contained in a single
element. For the UDG-FEM approach, the same holds with
respect to a single cut cell. This is especially advantageous
when solving the EEG forward problem for many dipoles
using a fast transfer matrix approach (see below). Note that
for tetrahedral elements, where we use the P1 basis, the
gradient of the basis function is constant. In this case, the
source model does not depend on the local position within its
supporting element. In addition, if the dipole lies exactly on a
cell boundary, we shift it slightly towards the interior of one
of the neighboring elements.

B. Transfer matrix
In order to reduce the computational load when solving
the EEG forward problem for many sources, we use a fast
transfer matrix approach [2], [12], [14] , which is also related
to the adjoint method [40]. In most cases, we are not interested
in the potential in the whole volume conductor, but only in
the potential difference at a set of Ne electrode positions
p1 , . . . , pNe ∈ R3 with respect to a reference electrode p0 . We
will denote the potential differences by U ∈ RNe . These values
can be obtained, by multiplying the solution vector y ∈ Rn
with a restriction matrix R ∈ RNe ×n : U = Ry. The entries of
the restriction matrix R = (rk,i ) are given as rk,i = ϕi (p0 ) −
ϕi (pk ). Replacing y by M −1 f we get U = RM −1 f = T f ,
with the transfer matrix T = RM −1 ∈ RNe ×n . T can be
computed by solving M T t = Rt (considering the symmetry
of M ). The latter can be carried out for each column of T t and
Rt separately. For a given discrete source model, the potential
differences at the electrode positions can now be computed by
a simple matrix vector multiplication.
C. Implementation
The method is implemented using the Distributed and
Unified Numeric Environment1 (DUNE) [41]–[44]. DUNE
is a general purpose, open source C++ library for solving
partial differential equations using mesh-based methods. It
is extensible by providing a modular structure and offering
generic interfaces and separation between data structures and
algorithms. For representing tetrahedral and hexahedral conforming meshes, we use the DUNE-ALUGrid module [45].
The discretization of the partial differential equation uses the
DUNE-PDELab module [46].
D. Solver methods
In order to solve the different linear equation systems, we
use different solvers for DG-FEM and UDG-FEM. For the
DG-FEM approach we use a conjugate gradient method with
an algebraic multigrid preconditioner (AMG) as presented in
[47]. Within this multigrid method, an SSOR preconditioner
is used as a smoother and a direct sparse solver is employed
as a coarse solver [48]. We used a V-cycle for the recursive
scheme and applied two pre- and two post-smoothing steps
on each level. For the UDG-FEM approach, a similar AMG
approach should in principle be possible but has not yet been
implemented so far. For UDG-FEM, we use here a conjugate
gradient method with a block incomplete LU decomposition as
a preconditioner on the cut cell blocks [43]. The iteration for
all methods was stopped at a relative reduction of the residual
kf − M yk2 of 10−8 .
E. Multi-layer sphere model
To verify and evaluate the new approach, we compared
the numerical solution to a quasi-analytical solution in a four
compartment sphere model [3]. The radii of the four layers and
the corresponding conductivity values are shown in Table I.
For the conductivity values, we followed the recommendations
of [27].
1 http://www.dune-project.org

5

TABLE I: Sphere radii, tissue labels and conductivity values
from outer to inner compartment.
Radius

92 mm

86 mm

80 mm

78 mm

Tissue

skin

skull

CSF

brain

Conductivity

0.43 S/m

0.01 S/m

1.79 S/m

0.33 S/m

TABLE II: Number of cut cells and degrees of freedom for
the different UDG-FEM discretizations. Note that the number
of degrees of freedom is given as DOFs = 8 · total.
Nd

skin

skull

CSF

brain

total

DOFs

24

1.328

1.120

872

1.520

4.840

38.720

25

6.872

5.776

4.064

10.552

27.264

218.112

26

37.864

32.800

18.184

78.000

166.848

1.334.784

F. Meshing aspects
For studying the convergence behavior of the UDG-FEM
approach as the number of elements is increased, we generated
a sequence of meshes with different resolutions. We started
with a fundamental mesh with Nd = 24 elements in each
dimension and refined the mesh by doubling the number
of elements in each dimension. Tab. II shows the resulting
number of cut cells and degrees of freedom (DOF; Nb = 8
DOFs on each cut cell). We denote these meshes by the
number of DOFs as: UDG 39k, UDG 218k and UDG 1335k
(see Figs. 4, 5 and 6). The mesh widths h of the three
fundamental meshes are approximately 12.13 mm, 6.06 mm
and 3.03 mm respectively.
For a comparison study of the UDG-FEM approach with
the DG-FEM approach presented by [26], we constructed for
the DG-FEM approach a conforming mesh of the multi-layer
sphere model with tetrahedral elements such that both methods
used a similar number of DOFs. The conforming mesh for
DG-FEM was generated using a constrained Delaunay tetrahedralization, as implemented in the tetgen software [21], [49].
The resulting tetrahedral DG-FEM mesh (DG tet 1447k had
1.446.804 DOFs, which is close to the highest mesh resolution
used for the UDG-FEM (UDG 1335k).
Furthermore, we created a regular hexahedral mesh of the
multi-layer sphere model with 2 mm resolution, resulting in
3.056.904 DOFs (DG hex 3057k). This regular hexahedral
mesh approximates the smooth surfaces of the spherical compartments only in a staircase like manner. The goal here is to
compare DG-FEM on a higher resolution regular hexahedral
mesh (DG hex 3057k) with UDG-FEM using much lower
resolution (UDG 39k) to show the significant contribution
of UDG-FEM to more accurately represent smooth tissue
surfaces and thus better approximate realistic head geometries.
Fig. 2 shows sections of the conforming tetrahedral and
hexahedral meshes as well as the finest cut cell mesh. For both
methods and all meshes, we computed the transfer matrices
and calculated the potentials at the electrodes for all sources
as explained in the following.

G. Sources
We generated two sets of dipoles with unit strength, one
with tangential and the other with radial orientation. For each
set, we generated 1000 dipoles at 10 different eccentricities in
the inner compartment. The location of the dipoles were computed randomly within the given eccentricity. The orientation
of the tangential dipoles were also chosen randomly in the
tangential plane. An eccentricity value of 0 denotes the center
of the sphere and a value of 1 a location at the boundary of the
inner compartment. Because it is well known that numerical
errors increase with increasing eccentricity (a reasoning for
this effect has been given in [38]), the chosen eccentricities
were scaled logarithmically with increasing eccentricity and
range from 0.1666 to 0.9939. The latter corresponds to a
distance of only 0.48 mm to the inner sphere surface and
thus very high eccentricity (where thus also higher numerical
errors have to be expected [38]). With regard to the application,
distances between 1 mm (eccentricity of 0.9872) and 2.5 mm
(0.9678) or 3 mm (0.9615) seem to be the most important,
because, depending on the location, the cortex has a thickness
of about 2 to 6 mm [50] and source locations should be chosen
to be in the middle of the grey matter compartment. Among the
10 chosen eccentricities, the values between 0.82 mm (0.9895)
and 2.45 mm (0.9686) come closest to this range, so that we
will later focus on these values and on especially the middle
value of this range, i.e., 1.42 mm (0.9818).
H. Error measures
The analytical and numerical potential solutions for each
dipole were evaluated at 200 electrodes on the outermost
surface of the sphere model, denoted by Uana ∈ R200 and
Unum ∈ R200 , respectively [3]. The error between both
solutions was computed using two common measures, the
relative difference measure (RDM%) and the magnitude error
(MAG%) [19]:
Unum
Uana
−
(18)
RDM%(Uana , Unum ) = 50
kUana k2
kUnum k2 2


kUnum k2
MAG%(Uana , Unum ) = 100
−1
(19)
kUana k2
The measures are scaled such that 0 ≤ RDM% ≤ 100
and −100 ≤ MAG% < ∞ hold. Note that the accuracy
of a method is better, the closer both measures are to zero.
Since the quality of the simulation is known to depend on the
local mesh geometry as well as on the intra-element position
of a dipole, RDM% and MAG% are analyzed statistically
for the above described sample of sources as suggested in
[19], including results for the different source eccentricities in
separate box-plots. This statistical analysis includes maximum
and minimum, indicated by upper and lower error bars, and
thereby the total range (TR). Furthermore, it includes the
interval between upper and lower quartile, i.e., the interquartile
range (IQR), also known as the spread, which is marked by a
box with a black dash showing the median.
I. Realistic head model
To construct a realistic four compartment head volume
conductor model, we used the data of [21] with an adapted

6

(a) Conforming tetrahedral mesh

(b) Conforming hexahedral mesh

(c) Cut cell mesh

Fig. 2: Sections of the different meshes used in the multi-layer sphere verification. From left to right, the images show the DG
tet 1447k, DG hex 3057k and UDG 1335k models. The different colors represent the different conductivity values.

preprocessing pipeline: T1-weighted (T1w-) and T2-weighted
(T2w-) MRI scans of a healthy 25-year-old male subject were
acquired in a 3 T MR scanner (MagnetomTrio, Siemens,
Munich, Germany) with a 32-channel head coil. For the T1wMRI, an MP-RAGE pulse sequence (TR/TE/TI/FA = 2300
ms/3.03 ms/ 1100 ms/8◦ , FOV = 256 × 256 × 192 mm,
voxel size = 1 × 1 × 1mm) with fat suppression and GRAPPA
parallel imaging (acceleration factor = 2) was used. For the
T2w image, an SPC pulse sequence (TR/TE = 2000 ms/307
ms, FOV = 255 × 255 × 176 mm, voxel size = 0.99 ×
1.0 × 1.0 mm interpolated to 0.498 × 0.498 × 1.00 mm)
was used. MR images were resampled to 1 mm isotropic
resolution. The T2w-MRI was registered onto the T1w-MRI
using a rigid registration approach and mutual information
as cost-function as implemented in the FSL-toolbox2 . The
skin and skull compartments were segmented by applying a
gray-value based active contour approach [51]. Subsequently,
the segmentation was manually corrected and, because of the
importance of modeling skull holes for source analysis [52],
[53], the foramen magnum and the two optic canals were
correctly modeled as skull openings. The model was not cut
off directly below the skull but realistically extended at the
neck. CURRY73 was used to extract high-resolution surfaces
of skin and skull compartments. A Taubin smoothing was
applied to remove staircase-like effects [54]. Cortex segmentation and surface reconstruction were performed using the
FreeSurfer-toolbox4 . From these surfaces, we computed signed
distance functions on a 2 mm structured mesh, which were
subsequently used as level set functions for defining the four
different tissue compartments skin, skull, csf and brain. The
same conductivity values were used for these compartments
as in the multi-layer sphere model. The basis functions were
defined on a mesh with a coarser resolution of 4 mm. A
section of the resulting mesh can be seen in Fig. 3. The
resulting discretization had 106.031 cut cells and thus 848.248
DOFs and skin, skull, csf and brain compartments consisted
of 38.541, 25.587, 18.147 and 23.756 cut cells, respectively
(head model REA 848k, see Fig. 7).
2 http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL
3 http://www.neuroscan.com
4 https://surfer.nmr.mgh.harvard.edu

IV. R ESULTS
A. Verification and evaluation in multi-layer sphere model
a) Convergence of the UDG-FEM approach: In a first
study, we verified the new UDG-FEM approach and examined
its convergence behavior for increasing mesh resolution.
The RDM% and MAG% statistical errors can be seen
in Fig. 4 (note: x-axes are logarithmically scaled). With
increasing source eccentricity, we observe an overall increase
in the RDM% error and an increasing total range (TR) and
interquartile range (IQR) of the MAG% error for all meshes.
With regard to the RDM%, the finest mesh achieves a
maximal error of 2.5 % for radial sources and 2.0 % for
tangential sources. Both of these values are observed at the
maximal eccentricity of 0.9939 (i.e., only 0.48 mm from the
inner sphere surface). The median value over all eccentricities
is 0.6 % for radial and for tangential sources. At the distance
of 1.42 mm (eccentricity value 0.9818), from coarse to fine
mesh, the IQR decreases from 7.7 % over 1.1 % down to 0.5
% for radial sources and from 5.5 % over 0.9 % down to 0.5
% for tangential sources. The TR behaves similarly.
For the MAG% error on the finest mesh, we observe
a maximal absolute value over all eccentricities of 8.8 %
for radial and 3.4 % for tangential sources. For radial and
tangential sources, the median value over all eccentricities is
0.8 %. When examining the MAG% at the distance of 1.42
mm (0.9818), from coarse to fine mesh, the IQR decreases
from 16.5 % over 1.0 % down to 0.8 % for radial sources and
from 4.9 % over 1.3 % down to 0.8 % for tangential sources.
The TR behaves again similarly.

Fig. 3: A section of the cut cell mesh for the realistic
head model: the thick black lines show the cut cells on the
fundamental mesh on which the basis functions are defined,
while the dashed black lines depict the subtriangulation used
for integration in the UDG-FEM approach.

7

20

(c) MAG% error for radial dipoles

0.9939

0.9818
0.9818

eccentricity

0.9686

0.9459

0.9065

0.9939

0.9895

0.9818

0.9686

0.9459

0.9065

15

0.8384

15

0.7209

10
0.8384

5

10
0.5177

0.9686

0

0.7209

5

eccentricity

0.9939

5

0

0.1666

UDG 39k
UDG 218k
UDG 1335k

10

MAG%

5

0.9895

15

UDG 39k
UDG 218k
UDG 1335k

10

MAG%

(b) RDM% error for tangential dipoles

0.1666

15

eccentricity

0.9895

(a) RDM% error for radial dipoles

0.9459

0.9939

0.9895

0.9818

eccentricity

0.9686

0.9459

0.9065

0.8384

0

0.7209

0

0.5177

5
0.1666

5

0.9065

10

0.8384

10

15

0.7209

15

0.5177

RDM%

20

UDG 39k
UDG 218k
UDG 1335k

25

0.5177

25

RDM%

30

UDG 39k
UDG 218k
UDG 1335k

0.1666

30

(d) MAG% error for tangential dipoles

Fig. 4: Verification and convergence study for the UDG-FEM approach: RDM% (upper row) and MAG% (lower row) errors
for radial (left column) and tangential (right column) sources for the three different resolutions 39k (green), 218k (red) and
1335k (blue). Note that the x-axis is logarithmically scaled.

b) Comparison of UDG-FEM to DG-FEM on conforming mesh with tetrahedral elements: In Fig. 5, we compare the
UDG-FEM approach using model UDG 1335k (blue) with the
DG-FEM approach on the conforming mesh with tetrahedral
elements using model DG tet 1447k (green). We show the
RDM% (upper row) and MAG% (lower row) errors for radial
(left column) and tangential (right column) sources (note: xaxes are again logarithmically scaled). The source model for
the DG-FEM approach with tetrahedral elements does not
depend on the local position of the dipole within its supporting
element, as the right hand side only involves the gradient of
the basis functions, which is constant for linear functions. We
therefore indicate by three vertical yellow lines the maximal
(6.70 mm at 0.9141), mean (4.77 mm at 0.9389) and minimal
(3.27 mm at 0.9581) element diameter of those elements in
the inner compartment which touch the brain/csf compartment
boundary. Note that this information is only relevant for the
DG-FEM approach with tetrahedral elements. As no dipoles
were located at exactly those values, we indicated instead their
closest eccentricities.
Since the UDG-FEM results were already discussed above,
we focus now on the DG tet 1447k results of Fig. 5. We

observe an overall increasing RDM% error up to the aforementioned left vertical line. A similar observation can be made
for the TR and IQR of the MAG% error. The maximal RDM%
error over all eccentricities is 6.3 % for radial dipoles and
7.8 % for tangential dipoles. For both source orientations, the
mean value over all eccentricities is 1.5 %. For radial sources,
the maximal absolute value of the MAG% error is 4.0 % and
for tangential sources it is 5.3 %.
Most importantly, with regard to RDM%, over all eccentricities and for both source orientations, with an average median
value of 0.7 %, a TR of 1.5 % and an IQR of 0.4 %, UDG
1335k performs better than DG tet 1447k, which only achieves
an average median value of 1.5 %, a TR of 5.8 % and an IQR
of 0.8 %. For the MAG% the situation is less clear with, for
UDG 1335k, an average median value of -0.7 %, a TR of 4.4
% and an IQR of 0.8 % and, for DG tet 1447k, an average
median value of -0.2 %, a TR of 6.0 % and an IQR of 1.4 %.
At the distance of 1.42 mm (0.9818), with regard to RDM%,
we find for radial sources for UDG-FEM a median of 0.7 %,
an IQR of 0.5 % and a TR of 1.9 % and a median of 1.1 %,
an IQR of 0.5 % and a TR of 4.3 % for the DG-FEM. For
tangential sources, a median of 0.9 %, an IQR of 0.5 % and

0.9895

0.9939

0.9895

0.9939

0.9818

0.9686

0.9459

0.9065

0.8384

eccentricity

(b) RDM% error for tangential dipoles

4

4

2

2

MAG%

0

2

0.9818

eccentricity

0.9686

0.9459

0.9065

0.8384

0.9939

0.9895

0.9818

0.9686

0.9459

eccentricity

(c) MAG% error for radial dipoles

0.7209

DG tet 1447k
UDG 1335k

4
0.9065

0.8384

0.7209

0.1666

4

0.5177

DG tet 1447k
UDG 1335k

0.5177

2

0

0.1666

MAG%

(a) RDM% error for radial dipoles

0.7209

DG tet 1447k
UDG 1335k

0.5177

8
7
6
5
4
3
2
1
0

0.1666

0.9939

0.9895

0.9818

eccentricity

0.9686

0.9459

0.9065

0.8384

0.7209

RDM%

DG tet 1447k
UDG 1335k

0.5177

8
7
6
5
4
3
2
1
0

0.1666

RDM%

8

(d) MAG% error for tangential dipoles

Fig. 5: Comparison between the UDG-FEM approach using model 1,335k (blue) and the DG-FEM approach on a conforming
mesh with tetrahedral elements using model 1,447k (green): RDM% (upper row) and MAG% (lower row) errors for radial
(left column) and tangential (right column) sources. Note that the x-axis is logarithmically scaled. With vertical yellow lines,
the maximal, mean and minimal (from left to right) element diameter is indicated of those elements in the inner compartment
touching the compartment boundary. This diameter is given relative to the radius of the inner sphere. Note that this information
is only relevant for the DG-FEM approach.

a TR of 1.7 % can be noted for UDG-FEM and a median of
1.6 %, an IQR of 0.7 % and a TR of 6.8 % for DG-FEM.
With regard to MAG%, we find for radial sources a median
of 0.2 %, an IQR of 1.4 % and a TR of 6.5 % for UDG-FEM
and a median of -0.5 %, an IQR of 0.7 % and a TR of 3.5 %
for the DG-FEM. For tangential sources we find a median of
-0.9 %, an IQR of 0.9 % and a TR of 3.6 % for UDG-FEM
and a median of 0.5 %, an IQR of 2.4 % and a TR of 6.0 %
for DG-FEM.
In summary, in the sphere model verification study the
UDG-FEM overall outperforms the DG-FEM. Note in addition, that in the case of a realistic head, the model generation
for the tetrahedral DG-FEM is far more sophisticated and,
because touching brain and skull surfaces are not allowed, the
volume conductor representation might be even less appropriate than it is for UDG-FEM.
c) Comparison of low-resolution UDG-FEM and highresolution DG-FEM on a structured mesh with hexahedral
elements: In Fig. 6, we compare the UDG-FEM approach

(green) with the DG-FEM approach on a structured mesh with
2 mm hexahedral elements (yellow) and show the RDM%
(upper row) and MAG% (lower row) errors for radial (left
column) and tangential (right column) sources. Note that the
x-axes are again logarithmically scaled. For the UDG-FEM
approach, we use the results from the coarsest mesh (39k
DOFs) in the first test, while the DG-FEM approach uses about
79 times the number of DOFs (3057k).
With regard to RDM%, even if for source eccentricities
below 2.45 mm (0.9686) DG hex 3057k has slightly lower
median, IQR and TR errors, for the more important higher
eccentricities, both methods perform similarly. With regard
to MAG%, the UDG-FEM performs even significantly better
than the much higher resolution DG-FEM: For radial dipoles
and source eccentricities below 2.45 mm (0.9686), the mean
of all median values is 1.7 % for UDG 39k, while it is 14.2
% for DG hex 3057k and for tangential sources, we find 0.5
% and 13.7 %, respectively. For radial sources at 1.42 mm
distance to the csf (0.9818), the MAG% median error of DG

9

50

(c) MAG% error for radial dipoles

0.9939
0.9939

0.9686

0.9818
0.9818

eccentricity

0.9686

0.9459

0.9065

0.9939

0.9895

0.9818

0.9686

0.9459

0.9065

100

0.8384

100

0.7209

50

0.5177

50

0.8384

0

0.7209

MAG%

0

0.1666

DG hex 3057k
UDG 39k

50

eccentricity

0.9895

100

DG hex 3057k
UDG 39k

50

MAG%

(b) RDM% error for tangential dipoles

0.1666

100

eccentricity

0.9895

(a) RDM% error for radial dipoles

0.9459

0.9939

0.9895

0.9818

eccentricity

0.9686

0.9459

0

0.9065

0

0.8384

10
0.7209

10
0.5177

20

0.1666

20

0.9065

30

0.8384

30

40

0.7209

RDM%

40

0.5177

50

DG hex 3057k
UDG 39k

60

0.5177

60

RDM%

70

DG hex 3057k
UDG 39k

0.1666

70

(d) MAG% error for tangential dipoles

Fig. 6: Comparison between the UDG-FEM approach with 38.720 DOF (green) and the DG-FEM approach on a structured
mesh with hexahedral elements and 3.056.904 DOF (yellow): RDM% (upper row) and MAG% (lower row) errors for radial
(left column) and tangential (right column) sources. Note that the x-axis is logarithmically scaled.

hex 3057k is at -20.6 %, while it is only at 9.1 % in case of
UDG 39k. Similarly, the IQR of DG hex 3057k of 75.3 % and
TR of 304.6 % is reduced to 18.6 % and 78.6 % in case of
UDG 39k, respectively. For tangential sources, the median of
12.5 % of DG hex 3057k is reduced to only 0.5 % for UDG
39k, the IQR of 20.2 % down to 8.1 % and the TR of 142.8
% down to 67.8 % in case of UDG 39k.
In summary, the UDG-FEM with only 39k DOFs achieves
already better or at least comparable results than the 2 mm
regular hexahedral DG-FEM approach with about 79 times
the number of DOFs (3057k). This test-case shows the high
importance of a good geometry-approximation as performed
by the UDG-FEM. With regard to application, higher resolution is needed to reduce the numerical errors especially for
highly eccentric sources. The UDG approach thus has to be
combined with a sufficient resolution, as also shown in the
first study (Fig. 4).
B. Forward simulation in the realistic head model
In the last study, we computed an example solution of the
EEG forward problem with the UDG-FEM approach for a
current dipole in the realistic four compartment head model

REA 848k. The result is shown in Fig. 7. The current dipole
is located in the auditory cortex and oriented normally to
the grey matter surface. The two images on the left show
the potential distribution on the scalp surface and the two
images on the right the different tissue compartments and
the potential distribution on a sagittal slice through the head
model. A smooth representation of the different tissue surfaces,
obtained from the level set functions, can be observed. The
forward simulation leads to a negativity in frontal areas and
a positivity in the lower left neck area which is in line
with practical findings for auditory evoked potentials [1]. The
sequential runtime for the solution of this example simulation
was approximately 2 minutes with a preceding 4 minutes of
setup time on a single core of a conventional laptop.
V. D ISCUSSION
In this paper, we presented the unfitted discontinuous
Galerkin finite element method (UDG-FEM) for solving the
EEG forward problem. In Section II, we provided a mathematical formulation of the underlying model and in Section
III a description of the overall method along with information
regarding its implementation.

10

Fig. 7: Exemplary EEG forward solution for an auditory source computed by UDG-FEM in the realistic four compartment
head model REA 848k. The two images on the left show the potential at the scalp surface and the two images on the right
the four tissue compartments and the potential on a sagittal slice through the head model.

In our verification and evaluation studies using multi-layer
sphere models in Section IV-A, the results first of all show
convergence of the UDG-FEM numerical solutions towards
the analytical solutions when the mesh resolution is increased
and, thereby, a better geometric approximation of the sphere
model is achieved (Fig. 4). Furthermore, we showed that
the accuracy of a discontinuous Galerkin FEM (DG-FEM)
approach on a conforming mesh with tetrahedral elements
could be met by our new UDG-FEM approach in models with
a comparable number of degrees of freedom (DOFs) (Fig. 5).
Most importantly, an even better accuracy was achieved by the
UDG-FEM when compared to DG-FEM on a structured mesh
with hexahedral elements, even though the UDG-FEM mesh
had a factor of about 79 times less DOFs than the DG-FEM
mesh (Fig. 6).
The comparison with the DG-FEM approach in Fig. 6
indicated the importance of the geometric representation of
the domain boundary. This is in accordance to the findings in
[26], where a first DG-FEM approach was presented for the
EEG forward problem based on the subtraction approach [38]
and where the new approach was verified in multi-layer sphere
models on structured meshes with hexahedral elements. In
contrast to Lagrange or Continuous Galerkin FEM (CG-FEM)
approaches [12]–[20], conservation properties could be proven
for the new DG-FEM [26]. As the DG-FEM formulation is
part of the UDG-FEM approach, the latter directly inherits the
conservation properties of the former. Furthermore, in [26],
the geometric representation of the domain boundary could
be identified as a main source of error for the numerical
scheme. Here, we reproduced this DG-FEM result in Fig. 6
and showed, that even on a very coarse mesh with about 79
times less DOFs, the UDG-FEM approach achieved a similar
accuracy when compared to DG-FEM. Therefore, when taking
the convergence results of Fig. 4 into account, it gets clear
that on finer meshes UDG-FEM outperforms DG-FEM on
structured meshes with hexahedral elements.
In [17], a trilinear immersed finite element method (TIFEM) has been proposed to solve the EEG forward problem.
Similar to the UDG-FEM approach, it uses a level set function
to represent the geometry and the basis functions are defined
on a structured mesh. Since it is based on a continuous
Galerkin approach, for the context of this paper, it is now
abbreviated as TI-CG-FEM. Like UDG-FEM, also TI-CG-

FEM modifies the local basis functions to incorporate the
geometric information, but in contrast to UDG-FEM and like
all other CG-FEM approaches, it does not fulfill the conservation property of electric charge. An additional important
difference is that the TI-CG-FEM approach does not change
the DOFs when compared with its corresponding CG-FEM
approach, while UDG-FEM leads to an increase in DOFs when
compared to its corresponding DG-FEM approach. In future
studies, besides investigating the implications of fulfilling or
not the conservation property and increasing or not the number
of DOFs, it might furthermore be interesting to compare UDGFEM and TI-CG-FEM with regard to thin skull compartments
and limited FEM resolutions. This was done for DG-FEM and
CG-FEM in [26], where skull leakage problems, appearing in
CG-FEM, could be strongly alleviated with DG-FEM. Such
scenarios might get relevant, e.g., in infant studies [53] or
for temporal bone areas, where skull thickness is 2 mm or
even less [24, Table 2]. When considering the comparison
CG-FEM and TI-CG-FEM on the one hand and DG-FEM
and UDG-FEM on the other, we find similar results with
regard to accuracies. Both, TI-CG-FEM and UDG-FEM show
better accuracies than their counterparts on structured hexahedral meshes. Furthermore, when compared to the results
on conforming tetrahedral meshes, the accuracies are met
by both approaches. A more direct comparison is, however,
difficult, since the parametrization of the multi-layer sphere
model for verifying TI-CG-FEM in [17] is different from ours
with regard to number of compartments (4 here versus 3 in
[17]), radii of the compartments, chosen conductivities and
source eccentricities (in [17], dipoles where considered along
a single axis up to 2 mm to the next conductivity jump, i.e.,
in the terms used here, up to an eccentricity of 0.977).
Another method to improve the geometric representation of
the standard CG-FEM approach on a structured hexahedral
mesh has been presented in [28]. In order to reduce the effect
of a staircase approximation of the different compartment
surfaces, the position of the mesh nodes is geometrically
adapted based on the different conductivities in the surrounding elements. The nodes are shifted towards the centroids of
their minority elements by a fixed global shift length. This
length has to be chosen such that the determinant of the
Jacobian of the element transformation remains positive. As
nodes are moved solely based on the voxel segmentation, the

11

topology of the mesh is restricted to the one given by the
segmentation and no further geometric information is included.
In [28], the geometry-adapted hexahedral CG-FEM method is
verified in a 4-layer sphere model, parameterized like here with
the only exception that skull conductivity was chosen 0.0042
S/m (note, the lower the skull conductivity, the higher the
expected numerical errors). A 2 mm mesh resolution was used
and sources were modeled with subtraction and Saint-Venant
approaches. Source positions where chosen along a single axis
up to 2mm to the next conductivity jump (eccentricity of
0.974), so we compare the results with the median values
of the 2 mm UDG-FEM model. As we used the partial
integration approach in the verification of UDG-FEM, we
compare the results with the Saint-Venant approach, which is
methodologically closer than the subtraction approach. With
respect to the RDM%, both methods perform reasonably well,
with slightly better accuracy for the UDG-FEM approach.
The difference of the MAG% error is more pronounced,
especially for radial sources. For the geometry adapted CGFEM, the MAG% error for radial sources is approximately
3.5 - 4 %, while it is at approximately -1 % (at eccentricity
0.9686) for the UDG-FEM approach. As presented in [26],
the CG-FEM method with hexahedral elements might suffer
from skull leakage currents on coarser meshes. The most
important difference might thus be that this can be alleviated
by employing the DG-FEM approach [26] and especially the
here presented UDG-FEM, but not by the CG-FEM geometryadaptation of [28].
A different approach to solving the EEG forward problem
is the boundary element method (BEM) [4]–[7]. This method
uses surface meshes and only takes the interfaces between
different homogenized tissue compartments into account. In
principle, almost arbitrary partitioning of the domain into different tissue compartments is possible [55]. In practice, however, mainly nested compartments with approximated holes are
implemented, as the head model generation process is a difficult problem (see, e.g., [53]). For the UDG-FEM method, no
such restrictions apply to the level set functions representing
the different tissue boundaries. UDG-FEM is especially suited
for such complex shaped domains as presented here for source
analysis and, for other applications, in [56]. Although it was
not directly part of this study, UDG-FEM can, like all FEM
approaches, handle realistic brain tissue anisotropy, which is
not the case for BEM.
We did not evaluate the computational load of the UDGFEM approach in this paper. When compared to the DGFEM approach on a mesh with the same resolution, the
UDG-FEM approach introduces more local basis functions
on cells which are cut by a level set. This increases the
computational load as more DOFs have to be considered.
However, the additional DOFs are well invested in areas of
special geometrical significance and as the UDG-FEM approach can achieve higher accuracies even on coarser meshes
(see Fig.6), the overall mesh resolution can thus even be
strongly reduced. As a consequence, it then leads to a drastic
reduction of the computational effort. Another source of
increased computational load can be found in the assembly of
the linear system. As the integration over cut cells introduces a

subtriangulation, more quadrature points are generated and the
bilinear form has to be evaluated more often. However, since
for the FEM transfer matrix approach, the linear system has
to be assembled only once and then solved for each electrode
[2], the resulting slightly increased overall setup time is rather
negligible. In a subsequent study, an evaluation of different
solution techniques for UDG-FEM, including multigrid methods and a thorough investigation of the performance compared
to competitive methods should be carried out. In this study,
we only performed an example simulation in a realistically
shaped head model. Due to the promising results in the multilayer sphere models, a thorough evaluation of UDG-FEM
in different realistic scenarios, especially with respect to the
above mentioned limitations of conforming tetrahedral and
hexahedral meshes, will be carried out in a subsequent study.
In addition, an investigation of the volumetric error distribution
in selected dipole and volume conductor configurations in
sphere models for the different methods, including CG-FEM
and TI-CG-FEM, using the analytic solution [3] would provide
further insight into the different model properties and will be
subject of our future work.
VI. C ONCLUSION
We presented theory, verification and evaluation of a new
unfitted discontinuous Galerkin finite element method (UDGFEM) for solving the EEG forward problem. In a multilayer sphere model, the accuracy shows proper convergence
behavior with increasing mesh resolution. When compared to
the discontinuous Galerkin FEM (DG-FEM) discretization on
a conforming tetrahedral mesh, the method performs similarly
or more accurately. On a structured mesh with hexahedral
elements, the UDG-FEM approach outperforms the DG-FEM
approach with respect to accuracy. The EEG forward simulation in a realistic head model for an auditory source resulted in
a smooth potential distribution over the head surface which is
in line with practical findings for auditory evoked potentials.
ACKNOWLEGMENT
This work was partially supported by the Cluster of Excellence 1003 of the Deutsche Forschungsgemeinschaft (DFG
EXC 1003 Cells in Motion), by EU project ChildBrain (Marie
Curie Innovative Training Networks, grant agreement no.
641652) and by the Priority Program 1665 of the Deutsche
Forschungsgemeinschaft (DFG) (WO1425/5-1).
R EFERENCES
[1] M. Hämäläinen et al., “Magnetoencephalographytheory, instrumentation, and applications to noninvasive studies of the working human
brain,” Reviews of Modern Physics, vol. 65, no. 2, p. 413, 1993.
[2] J. de Munck et al., EEG and MEG – forward modelling.
In R.Brette & A.Destexhe (eds.): Handbook of Neural Activity
Measurement. Cambridge University Press, ISBN 978-0-521-516228, 2012, chapter doi: 10.1017/CBO9780511979958.006. [Online].
Available: http://www.di.ens.fr/∼brette/HandbookMeasurement/
[3] J. de Munck and M. J. Peters, “A fast method to compute the potential
in the multisphere model,” IEEE Trans.Biomed.Eng., vol. 40, no. 11, pp.
1166–1174, 1993.
[4] J. C. Mosher et al., “EEG and MEG: forward solutions for inverse
methods,” IEEE Trans.Biomed.Eng., vol. 46, no. 3, pp. 245–259, 1999.

12

[5] Z. A. Acar and S. Makeig, “Neuroelectromagnetic forward head modeling toolbox,” J. Neurosci. Meth., vol. 190, no. 2, pp. 258–270, 2010.
[6] A. Gramfort et al., “Forward field computation with OpenMEEG,”
Comp.Intell.Neurosci., vol. 2011, 2011.
[7] M. Stenroos and J. Sarvas, “Bioelectromagnetic forward problem: isolated source approach revis (it) ed,” Phys. Med. Biol., vol. 57, no. 11,
p. 3517, 2012.
[8] M. J. Cook and Z. J. Koles, “A high-resolution anisotropic finite-volume
head model for EEG source analysis,” in Engineering in Medicine and
Biology Society, 2006. EMBS’06. 28th Annual International Conference
of the IEEE. IEEE, 2006, pp. 4536–4539.
[9] K. Wendel et al., “The influence of CSF on EEG sensitivity distributions
of multilayered head models.” IEEE Trans. Biomed. Eng., vol. 55, no. 4,
pp. 1454–1456, 2008.
[10] F. Vatta et al., “Solving the forward problem in EEG source analysis by
spherical and fdm head modeling: a comparative analysis-biomed 2009.”
Biomed.Sci.Instrument., vol. 45, pp. 382–388, 2008.
[11] V. Montes-Restrepo et al., “Influence of skull modeling approaches on
EEG source localization,” Brain Topogr., vol. 27, pp. 95–111, 2014.
[12] D. Weinstein et al., “Lead-field bases for electroencephalography source
imaging,” Ann. Biomed.Eng., vol. 28, no. 9, pp. 1059–1066, 2000.
[13] P. H. Schimpf et al., “Dipole models for the EEG and MEG,” IEEE
Trans. Biomed. Eng., vol. 49, no. 5, pp. 409–418, 2002.
[14] N. Gencer and C. Acar, “Sensitivity of EEG and MEG measurements
to tissue conductivity.” Phys.Med.Biol., vol. 49, pp. 701–717, 2004.
[15] C. Wolters et al., “Numerical approaches for dipole modeling in finite
element method based source analysis.” Int.Cong.Ser., vol. 1300, pp.
189–192, June 2007.
[16] S. Lew et al., “Improved EEG source analysis using low resolution
conductivity estimation in a four-compartment finite element head
model.” Human Brain Mapping, vol. 30, no. 9, pp. 2862–2878, 2009.
[17] S. Vallaghé and T. Papadopoulo, “A trilinear immersed finite element
method for solving the electroencephalography forward problem,” SIAM
J.Sci.Comp., vol. 32, no. 4, pp. 2379–2394, 2010.
[18] S. Pursiainen et al., “Forward simulation and inverse dipole localization
with the lowest order Raviart-Thomas elements for electroencephalography,” Inverse Problems, vol. 27, no. 4, p. 045003, 2011.
[19] J. Vorwerk et al., “Comparison of boundary element and finite element
approaches to the EEG forward problem,” Biomed. Eng., vol. 57, no.
SI-1, pp. 795–798, 2012.
[20] T. Medani et al., “FEM method for the EEG forward problem and
improvement based on modification of the Saint Venants method,”
Progress Electromag. Res., vol. 153, pp. 11–22, 2015.
[21] J. Vorwerk et al., “A guideline for head volume conductor modeling in
EEG and MEG,” NeuroImage, vol. 100, pp. 590–607, 2014.
[22] E. Paraskevopoulos et al., “Musical expertise is related to altered
functional connectivity during audiovisual integration,” PNAS, vol. 112,
no. 40, pp. 12 522–27, 2015.
[23] Ü. Aydin et al., “Combined EEG/MEG can outperform single modality
EEG or MEG source reconstruction in presurgical epilepsy diagnosis,”
PLoS ONE, vol. 10, no. 3, p. e0118753, 2015.
[24] J.-H. Kwon et al., “The thickness and texture of temporal bone in
brain CT predict acoustic window failure of transcranial doppler.” J.
Neuroimag., vol. 16, pp. 347–52, 2006.
[25] H. Sonntag et al., “Leakage effect in hexagonal FEM meshes of the EEG
forward problem,” Clinical EEG and Neuroscience, p. P022, 2013.
[26] C. Engwer et al., “A discontinuous Galerkin method for the EEG
forward problem,” ArXiv e-prints, Nov. 2015, submitted. [Online].
Available: arXiv:1511.04892
[27] M. Dannhauer et al., “Modeling of the human skull in EEG source
analysis,” Human Brain Mapping, vol. 32, no. 9, pp. 1383–1399, 2011.
[28] C. H. Wolters et al., “Geometry-adapted hexahedral meshes improve
accuracy of finite-element-method-based EEG source analysis.” IEEE
Trans. Biomed. Eng., vol. 54, no. 8, pp. 1446–1453, 2007.
[29] P. Bastian and C. Engwer, “An unfitted finite element method using
discontinuous Galerkin,” Int. J. Numer. Methods Eng., vol. 79, no. 12,
pp. 1557–1576, 2009.
[30] F. Heimann et al., “An unfitted interior penalty discontinuous Galerkin
method for incompressible Navier–Stokes two-phase flow,” Int. J. Numer.
Meth. in Fluids, vol. 71, no. 3, pp. 269–293, 2013.
[31] C. Engwer and S. Westerheide, “Heterogeneous coupling for implicitly
described domains,” in Domain Decomposition Methods in Science and
Engineering XXI. Springer, 2014, pp. 809–817.
[32] E. Burman et al., “A Cut Discontinuous Galerkin Method for
the Laplace-Beltrami Operator,” ArXiv e-prints, Jul. 2015. [Online].
Available: arXiv:1507.05835

[33] A. Ern et al., “A discontinuous Galerkin method with weighted averages
for advection–diffusion equations with locally small and anisotropic
diffusivity,” IMA J. Num. Anal., 2008.
[34] X. Han et al., “CRUISE: cortical reconstruction using implicit surface
evolution,” NeuroImage, vol. 23, no. 3, pp. 997–1012, 2004.
[35] T. Papadopoulo and S. Vallaghé, “Implicit meshing for finite element
methods using levelsets,” in Computer Vision, 2007. ICCV 2007. IEEE
11th International Conference on. IEEE, 2007, pp. 1–8.
[36] C. Engwer and A. Nüßing, “Geometric integration over irregular
domains with topologic guarantees,” ArXiv e-prints, Jan. 2016,
submitted. [Online]. Available: arXiv:1601.03597
[37] I. O. Tanzer et al., “Representation of bioelectric current sources using
Whitney elements in the finite element method,” Phys.Med.Biol., vol. 50,
no. 13, p. 3023, 2005.
[38] C. Wolters et al., “Numerical mathematics of the subtraction method
for the modeling of a current dipole in EEG source reconstruction using
finite element head models.” SIAM J.Sci.Comp., vol. 30, no. 1, pp. 24–
45, 2007.
[39] J. Vorwerk, “New finite element methods to solve the EEG/MEG forward problem,” Dissertation, Fachbereich Mathematik und Informatik,
Westfälische Wilhelms-Universität Münster, Germany, February 2016.
[40] S. Vallaghé et al., “The adjoint method for general EEG and MEG
sensor-based lead field equations,” Physics in Medicine and Biology,
vol. 54, no. 1, p. 135, 2008.
[41] P. Bastian et al., “A generic grid interface for parallel and adaptive
scientific computing. Part I: Abstract framework,” Computing, vol. 82,
no. 2–3, pp. 103–119, 2008.
[42] P. Bastian et al, “A generic grid interface for parallel and adaptive
scientific computing. Part II: Implementation and tests in DUNE,”
Computing, vol. 82, no. 2–3, pp. 121–138, 2008.
[43] M. Blatt and P. Bastian, “The iterative solver template library,” in
Applied Parallel Computing. State of the Art in Scientific Computing,
ser. Lecture Notes in Computer Science, B. Kågström, E. Elmroth,
J. Dongarra, and J. Waśniewski, Eds., vol. 4699. Springer, 2007, pp.
666–675.
[44] C. Engwer and F. Heimann, “DUNE-UDG: a cut-cell framework
for unfitted discontinuous Galerkin methods,” in Advances in DUNE.
Springer, 2012, pp. 89–100.
[45] A. Dedner et al., “The DUNE-ALUGrid module,” arXiv preprint
arXiv:1407.6954, 2014.
[46] P. Bastian et al., “Generic implementation of finite element methods in
the distributed and unified numerics environment (DUNE),” Kybernetika,
vol. 46, no. 2, pp. 294–315, 2010.
[47] P. Bastian et al., “Algebraic multigrid for discontinuous Galerkin
discretizations of heterogeneous elliptic problems,” Numerical Linear
Algebra with Applications, vol. 19, no. 2, pp. 367–388, 2012.
[48] T. A. Davis, “Algorithm 832: Umfpack v4.3—an unsymmetric-pattern
multifrontal method,” ACM Trans. Math. Softw., vol. 30, no. 2, pp. 196–
199, Jun. 2004.
[49] H. Si, “Tetgen, a delaunay-based quality tetrahedral mesh generator,”
ACM Trans. Math. Softw., vol. 41, no. 2, pp. 11:1–11:36, Feb. 2015.
[50] G. Li et al., “Measuring the dynamic longitudinal cortex development
in infants by reconstruction of temporally consistent cortical surfaces,”
NeuroImage, vol. 90, pp. 266–279, 2014.
[51] L. A. Vese and T. F. Chan, “A multiphase level set framework for image
segmentation using the Mumford and Shah model,” Int.J.Comp.Vis.,
vol. 50, no. 3, pp. 271–293, 2002.
[52] R. Oostenveld and T. F. Oostendorp, “Validating the boundary element
method for forward and inverse EEG computations in the presence of a
hole in the skull,” Human Brain Mapping, vol. 17, no. 3, pp. 179–192,
2002.
[53] Roche-Labarbe et al., “High-resolution electroencephalography and
source localization in neonates,” Human Brain Mapping, vol. 29, no. 2,
pp. 167–176, 2008.
[54] G. Taubin, “A signal processing approach to fair surface design,” in
Proceedings of the 22nd annual conference on Computer graphics and
interactive techniques. ACM, 1995, pp. 351–358.
[55] J. Kybic et al., “Generalized head models for MEG/EEG: Boundary
element method beyond nested volumes,” Phys.Med.Biol., vol. 51, no. 5,
p. 1333, 2006.
[56] C. Engwer, “An unfitted discontinuous Galerkin scheme for micro-scale
simulations and numerical upscaling,” Ph.D. dissertation, Heidelberg
University, 2009.

