Technical Article

Recent Advances of Isogeometric Analysis in Computational
Electromagnetics

arXiv:1709.06004v1 [] 18 Sep 2017

Z. Bontinck, J. Corno, H. De Gersem, S. Kurz, A. Pels, S. Schöps, F. Wolf, C. de Falco, J. Dölz, R. Vázquez and U. Römer

Abstract — In this communication the advantages and drawbacks
of the isogeometric analysis (IGA) are reviewed in the context of
electromagnetic simulations. IGA extends the set of polynomial
basis functions, commonly employed by the classical Finite Element Method (FEM). While identical to FEM with Nédélec’s basis
functions in the lowest order case, it is based on B-spline and NonUniform Rational B-spline basis functions. The main benefit of
this is the exact representation of the geometry in the language of
computer aided design (CAD) tools. This simplifies the meshing
as the computational mesh is implicitly created by the engineer
using the CAD tool. The curl- and div-conforming spline function
spaces are recapitulated and the available software is discussed. Finally, several non-academic benchmark examples in two and three
dimensions are shown which are used in optimization and uncertainty quantification workflows.

I. Introduction
In electrical engineering numerical simulations have become an
inevitable tool for designing new components, such as electrical
machines or antennas, as well as for understanding their behavior
and their interaction with the environment, as e.g. commonly
analyzed in electromagnetic compatibility simulations. However, due to exponentially increasing computational resources
and increasing accuracy expectation of the simulation engineers,
the numerical models become correspondingly more and more
complex.
There are many trends leading to this increased complexity and
the following list is obviously not complete. For example, multiphysical phenomena are considered to be increasingly important and they lead to many challenges in the modeling and the
numerical analysis, see e.g. [1]. On the other hand, the geometrical design of components to be simulated has become
very complex. It is often directly imported from a Computer
Aided Design (CAD) software into the engineering tools that
are the main interest of the Compumag community. Geometrical simplifications of the design, e.g. removing screws, require a
lot of manual labor and sometimes even insights into the actual
field distribution within the component, which should have been
gained by the simulation in the first place. According to [2],
approximately 80 % of the overall analysis time nowadays can
be accounted to mesh generation in the automotive, aerospace,
and ship building industries. At the same time, the higher accuracy demands require even finer meshes or higher order ansatz
functions, which in turn necessitate the creation of meshes with
curved elements. Now, if the initial mesh generation for realistic
computational models is already troublesome, how can one robustly carry out shape optimization or uncertainty quantification
of geometrical manufacturing imperfections?
In this communication, we advertise Isogeometric Analysis
(IGA) as proposed by Tom Hughes in [2], see e.g. [3]. Its
distinct feature is to mitigate the meshing step and immediately
work on the CAD geometry description. The idea of more
sophisticated parametric, i.e. non-polynomial, mappings is not
new but gained further interest also by CAD software vendors
due to the recent works on IGA. Eventually, the mesh shall be

C (ξ)
w=1

O
Fig. 1. A NURBS curve in R2 is the projection of a B-spline curve in R3 , this
allows for the definition of circles and other conic sections.

automatically created by the engineer using the CAD software.
Currently, this is not entirely true as boundary representations
(BRep) of volumetric objects are used in most CAD tools. The
trivariate mappings of the interior have still to be created manually. However, this is an active and ongoing field of research
[4].
Eventually, the Finite Element Method (FEM) is the work
horse in low-frequency and many high-frequency electromagnetic simulations [5, 6]. In the last decades, many variants and
improvements of the Finite Element Method have been proposed.
IGA proposes to generalize the set of polynomial basis functions, employed by FEM, with the introduction of more general
B-spline basis functions and Non-Uniform Rational B-splines
(NURBS). These functions are well known in the design community and are the basic ingredient of nowadays CAD software.
They allow for the exact parametrization of common curves and
surfaces such as conic sections, and are extremely flexible and
intuitive when dealing with more complex shape creation and
deformation. By using them for the discretization process of
general Partial Differential Equations (PDEs), it is possible to
inherit these properties and directly use the CAD design as the
computational domain, without the need of a meshing step.
It is also worth pointing out that, since IGA operates in the
same Galerkin framework as FEM, in most instances pre-existent
finite element codes can be easily modified in order to work in an
isogeometric setting by changing the basis function construction
routines only.
The possibility to straightforwardly deal with geometrical
changes makes of IGA a powerful tool when used in the context of shape optimization problems or shape sensitivity analysis [3, 7, 8]. Moreover, given the higher global regularity of the
basis function space, the isogeometric approach has been shown
to present several advantages over FEM in addition to the better
handling of geometries, like a faster convergence with respect to
the number of degrees of freedom [9] and the possibility to treat
high order differential operators [10, 11].
These properties have made IGA appealing for a wide variety of
applications. For the application to computational electromagnetics, in particular, IGA shows the ability of consistently discretizing complexes of differential forms [12] which is a property

1
Bi1 (ξ)

Bi2 (ξ)

1

0

0
0

1/3

2/3

1

0.2

0.4

0.6

0.8

1

ξ

1

2

y

Bi2 (ξ)
0

0

0

1/3

2/3

1

1

ξ
Fig. 2. B-spline basis functions of degree 1 and 2 on open, uniform knot
vectors (Ξ = [0, 0, 1/3, 2/3, 1, 1] on top and Ξ = [0, 0, 0, 1/3, 2/3, 1, 1, 1] at
the bottom).

of great importance for achieving spectrally correct discretization of the Maxwell PDEs [13, 14].
In this communication we present a general introduction to the
isogeometric framework highlighting some differences, advantages and drawbacks with respect to FEM. We also introduce
the curl- and div-conforming spaces necessary for the correct
discretization of Maxwell’s equations as given in [13]. We then
proceed to present several applications of IGA to different kind
of problems both in 2D and in 3D. Some final remarks on an
isogeometric boundary element method are also given.
II. Isogeometric Discretization of Maxwell’s Equations
The differential form of Maxwell’s equations is given by
∂B
=0
∂t
∇·D=ρ
∂D
∇×H−
=J
∂t
∇ · B = 0,
∇×E+

(1a)
(1b)
(1c)
(1d)

with E and D the electric field strength and electric flux density,
H and B the magnetic field strength and magnetic flux density, ρ
the electric charge density and J the current density. The system
is completed by the material relations
D = E + P
B = µ (H + M)

(2a)
(2b)

where the electric permittivity  and the magnetic permeability
µ are, in general, non-linear tensor valued functions of position
expressing the material properties, and the quantities P and M
are the electric polarization and the magnetization respectively.
Depending on the problem at hand, different assumptions and
simplifications can be made. We focus here on high-frequency
problems in frequency domain and low-frequency problems that
can be considered quasi-static or even static. We postpone the detailed discussion of the specific formulations and models to section IV. In general, however, the numerical solution of such problems is nowadays typically performed with the FEM. The spaces
arising from the weak formulation are non-standard spaces of
square integrable functions (i.e. in L2 ) with weakly defined curl
in L2 ; this space is commonly denoted by H (curl; ·). In order to
properly approximate the solution field, the discrete spaces need
to mimic several important properties that hold on the continuous level (see sub-section B). For classical polynomial FEM it is

0

0

1

2

3

x
Fig. 3. Continuity can be controlled in single points (see the black circle) by
knot repetition (here: Ξ = [0, 0, 0, .2, .4, .4, .6, .8, 1, 1, 1]).

well known that a sequence of discrete spaces with conforming
discretization can be obtained using (high-order) Nédélec-type
elements, which allocate the degrees of freedom on the edges or
the faces of the mesh in order to ensure consistency [5, 6].
A. IGA Basis Functions
In order to represent complex shapes, the use of polynomials or
rational segments may often be inadequate or imprecise. On the
other hand, B-spline and NURBS functions enjoy some major
advantages that make them extremely convenient for surface
representation and are therefore the most common choice for
solid geometry modeling on which nowadays CAD tools are
based. Of the main advantages we care to mention, e.g., that they
can exactly represent all conic sections (i.e. circles, ellipses, etc.),
see e.g. Fig. 1, that they can be generated by many efficient and
numerically stable algorithms, and that they can easily handle
specified continuity in single points [15].
The main idea behind the isogeometric approach [2] is to discretize the problem unknowns with the same set of basis functions that CAD employs for the construction of geometries.
Let p be the prescribed degree and


Ξ = ξ1 . . . ξn+p+1

(3)

be a vector that partitions [0, 1] into elements (ξi ∈ Ω̂ = [0, 1]).
Then, the Cox-de Boor’s formula [16] defines n one-dimensional
B-spline basis functions {Bip (ξ)} with ξ ∈ Ω̂ as
(
1 if ξi ≤ ξ < ξi+1
Bi0 (ξ) =
(4)
0 otherwise
Bip (ξ) =

ξi+p+1 − ξ
ξ − ξi
B p−1 (ξ) +
B p−1 (ξ) ,
ξi+p − ξi i
ξi+p+1 − ξi+1 i+1

with i = 1, . . . , n. An example of B-spline basis of degree
p = 1, 2 are shown in Fig. 2. One can notice that the support
of a B-spline of degree p is always p + 1 knot spans and, as a
consequence, each p-th degree function has p − 1 continuous
derivatives across the element boundaries (i.e. across the knots)
if they are not repeated. Repetition of knots can be exploited to
prescribe the regularity.
NURBS of degree p are defined as rational B-splines
wi Bip (ξ)
Nip (ξ) = P
,
p
j wj Bj (ξ)

(5)

be due, for example, to topological reasons or to the presence
of different materials. In these cases it is common practice to
resort to the so called multipatch approach [3]. Here the physical
domain is split into simpler subdomains Ωk such that each one
of them is the image of the reference space through a map Fk of
the type given in (9). An example of such a situation is depicted
in Fig. 4 on the right.

Ω3

Ω4

Ω2
Ω1

F

b
Ω

b
Ωk = Fk (Ω)

Fig. 4. The mesh in the physical space is given by the transformation through the
NURBS mapping F of the knot lines in the reference domain (on the left). On the
rigth an example of a multipatch mapping to represent complicated structures.

with wi a weighting parameter associated with the i-th basis
function. It is clear that B-splines are NURBS with all weights
equal to one, due to the partition of unity property.
Multi-dimensional B-splines and NURBS are constructed using
a tensor product approach. For example, let Ξd be the knot
vectors, pd the degrees and nd the number of basis functions
(with d = 1, 2, 3), a trivariate B-spline is given by
p3
p2
p1
Bp
i (ξ) = Bi1 (ξ1 ) Bi2 (ξ2 ) Bi3 (ξ3 ) ,

(6)

where p = (p1 , p2 , p3 ) and i = (i1 , i2 , i3 ) is a multi-index in
the set
I = {(i1 , i2 , i3 ) : 1 ≤ id ≤ nd } .
(7)
We will denote by Sp (Ξ) the multivariate B-spline space
spanned by the basis functions Bp
i .
To build a B-spline or NURBS curve, we start by defining a set of
control points. These act as weights for the linear combination
of the basis functions, giving the mapping to the physical space.
In particular, given n one-dimensional basis functions Nip and n
control points P ∈ Rd , i = 1, . . . , n, a curve parameterization
is given by:
n
X
Pi Nip (ξ) .
(8)
C (ξ) =
i=1

The control points define the so called control mesh but this does
not, in general, conform to the actual geometry. On the contrary,
the physical mesh is a decomposition of it: the mesh element
edges are the image of the knot lines through the mapping (8)
(see Fig. 4 on the left).
It is straightforward to notice that the continuity of a CAD curve
is controlled by the basis functions (by knot repetitions in particular), while the control points define the shape without altering
the curve continuity. Moreover, as a consequence of the locality
of the basis functions, moving a single control point can affect
the geometry of no more than p + 1 elements of the curve.

The main advantage of using NURBS over B-spline curves is
the possibility to exploit both the control points and the weights
in (5) to control the local shape: as wi increases, the curve is
pulled closer to the control point Pi , and vice versa. While
non-rational splines (or Bézier curves) can approximate a circle,
they are unable to represent it exactly. Rational splines, however,
overcome this issue.
In many real-world applications, the computational domain may
be too complicated to be represented by a single NURBS mapping from the reference domain to the physical space. This could

Finally, we care to mention that the possibility to use the same
space for the parameterization of geometry and solution (the so
called isoparametric approach), is particularly interesting when
dealing with shape deformations, either coming from the solution of mechanical problems or from optimization or sensitivity
analysis procedures. The deformation degrees of freedom can
simply be added to the geometry control points to obtain the deformed domain, and the internal parameterization automatically
follows with no need of cumbersome remeshing procedures.
However, we will see that in the case of Maxwell the isoparametric concept must be relaxed as the derivative of a NURBS
is not a NURBS function. Therefore, one uses NURBS for the
mapping and B-splines as weighting and ansatz functions.

B. IGA Conforming Spaces
Analogously to classical FEM, IGA is (typically) based on a
Galerkin approach: the equations are written in their variational
formulation, and the solution is sought in a finite dimensional
space with the correct approximation properties. In IGA, however, the basis function space is inherited from the one used to
parametrize the geometry.
Let us now consider a domain Ω ∈ Rd that can be exactly
parametrized with a mapping F of the type in (8), i.e.
F : Ω̂ → Ω,

(9)

with Ω̂ the reference domain [0, 1]d . We will denote by DJ the
Jacobian of the transformation.
We define the following pull-back functions
ι0 (v) := v ◦ F

v ∈ H 1 (Ω)

>

ι1 (v) := (DF) (v ◦ F)

−1

ι2 (v) := det (DF) (DF)

ι3 (v) := det (DF) (v ◦ F)

(10a)

v ∈ H (curl; Ω) (10b)

(v ◦ F) v ∈ H (div; Ω)
2

v ∈ L (Ω) .

(10c)
(10d)

It can be shown that (10b) and (10c) preserve the curl and the
divergence, respectively, from the reference domain to the physical one [5, 6]. Due to this property, the following commuting
de Rham diagram holds:
 
H 1 Ω̂

ˆ
∇

ι0

H 1 (Ω)



H curl; Ω̂

ˆ
∇×

ι1

∇

H (curl; Ω)



H div; Ω̂

ˆ
∇·

ι3

ι2

∇×

H (div; Ω)

 
L2 Ω̂

∇·

(11)

L2 (Ω)

Here we have also introduced the usual Sobolev space H 1 of
L2 functions with square integrable gradient and the space
H (div; ·) of functions with weak divergence in L2 .

To obtain a conforming discretization of H(curl, Ω), we search
an analogous diagram in the discrete setting. By exploiting the
fact that the derivative of a B-spline function is still a B-spline
function [16], we start defining the sequence of B-spline spaces

on the reference domain Ω̂
S 0 (Ω̂) = Sp (Ξ)

(12)

1

S (Ω̂) = Sp1 −1,p2 ,p3 (Ξ) × Sp1 ,p2 −1,p3 (Ξ)

(13)

× Sp1 ,p2 ,p3 −1 (Ξ)

2

S (Ω̂) = Sp1 ,p2 −1,p3 −1 (Ξ) × Sp1 −1,p2 ,p3 −1 (Ξ)

(14)

× Sp1 −1,p2 −1,p3 (Ξ)

3

S (Ω̂) = Sp−1 (Ξ).

(15)

It has been proven [13] that, using these spaces, a discrete counterpart to (11) can be constructed
S 0 (Ω̂)

ˆ
∇

S 1 (Ω̂)

ˆ
∇×

S 2 (Ω̂)

ˆ
∇·

S 3 (Ω̂).

(16)

To define the spaces in the physical domain, we use the pullbacks (10), i.e. a conforming discretization of H (curl; Ω) is
given by
n
o
1
S 1 (Ω) = v = ι−1
(v̂),
v̂
∈
S
(
Ω̂)
.
(17)
1
In case of a multipatch domain, e.g. Fig. 4, a global discretization space is constructed. Neighbouring patches are required to
share either a full edge or a full face, i.e. no T-junctions are
allowed. On each patch the matrix assembly is then performed
independently and the common degrees of freedom are matched
one-to-one through static condensation. It is clear that the multipatch approach reduces the global regularity to C 0 , although
inside each patch the discretization remains highly smooth.
III. Available Software
For those researchers interested in IGA and wanting to test how
it works in practice, a perfect way to start is GeoPDEs [19,
20], an open source and freely distributed package written in
MATLAB/Octave language [21, 22]. The GeoPDEs package
was developed with a double aim. First, to serve as a didactic
tool to introduce IGA to other researchers and students. For
this reason the package contains a long list of examples, and
all the functions include a detailed documentation accessible
from MATLAB with the help command. Secondly, to serve
as a research tool for fast prototyping and for testing new ideas
and methods in IGA, and that is the main reason why the code
is developed in MATLAB, which is a de facto standard for
prototyping of numerical algorithms.
GeoPDEs contains all what is required for the implementation
of IGA: evaluation of B-splines and NURBS functions, matrix
assembly, imposition of boundary conditions, etc. It also contains functions to export the numerical results to ParaView [23]
for post-processing. Solving a new problem can often be easily
accomplished by calling existing functions, which only require
minor modifications with respect to the already existing examples. For example, most applications we present in section IV
have been implemented with GeoPDEs.
One of the most interesting features is that, as far as we know,
GeoPDEs is the only MATLAB software that contains the spline
complex of Section II, which is crucial for the application of IGA
in computational electromagnetism. Moreover, GeoPDEs is
under constant development, and new features appear from time
to time. For instance, a very recent addition is the introduction of
IGA adaptive methods based on hierarchical B-splines (see [24]
for the details).
In Listing 1, we report a simple example of how to solve
Maxwell’s eigenvalue problem on a geometry given in the file

geo_file.mat. First we extract the knot vectors defining the
geometry, we create a refinement and modify them like in (13).
Then we construct the mesh and the curl-conforming space (17)
(line 23). The functions op_curlu_curlv_tp and op_u_v_tp
are responsible for the assembly of the curl-curl matrix K and
the mass matrix M respectively, exploiting the tensor product
structure. In lines 33-38 we exclude from the computation the
PEC degrees of freedom at the boundary and finally we call eig
to compute the eigenmodes.
GeoPDEs provides a compromise between clarity and efficiency. It can be applied to the solution of relatively large
three-dimensional problems such as those in Section IV. For
researchers interested in even larger problems, or more efficient
implementations, there are several libraries written in C or C++.
Here, we care to mention, igatools [25], G+Smo [26] and
PetIGA [27] (which recently added curl- and div-conforming
spline discretizations [28]). A more detailed list of available
IGA software can be found in [29].
Listing 1: A simple example of how to solve Maxwell’s eigenvalue problem in GeoPDEs
1
2
3
4
5

% Load Geometry
geo = geo_load (’geo_file .mat ’);

% Define Hcurl conforming knot v e c t o r s
[knots , zeta] = kntrefine (geo.nurbs .knots ,
nsub , degree , regularity );
6 [ knots_hcurl , degree_hcurl ] = knt_derham (
knots , degree , ’Hcurl ’);
7
8 % C o n s t r u c t t h e Mesh
9 rule = msh_gauss_nodes (nquad);
10 [qn , qw] = msh_set_quad_nodes (zeta , rule);
11 msh = msh_cartesian (zeta , qn , qw , geo);
12
13 % C o n s t r u c t t h e S p a c e
14 scalar_spaces = cell (msh.ndim , 1);
15 f o r idim = 1: msh.ndim
16
scalar_spaces {idim} = sp_bspline (
knots_hcurl {idim}, degree_hcurl {idim},
msh);
17 end
18 space = sp_vector ( scalar_spaces , msh , ’curl
- preserving ’);
19
20 % A s s e m b l e t h e m a t r i c e s
21 K = op_curlu_curlv_tp (space , space , msh);
22 M = op_u_v_tp (space , space , msh , @(x,y) mu0
*eps0* o n e s ( s i z e (x)));
23
24 % A p p l y PEC Boundary C o n d i t i o n s
25 drchlt_dofs = [];
26 f o r iside = 1: numel (space. boundary )
27
drchlt_dofs = union ( drchlt_dofs , space .
boundary (iside).dofs);
28 end
29 int_dofs = setdiff (1: space.ndof ,
drchlt_dofs );
30
31 % S o l v e t h e E i g e n v a l u e Problem
32 eigv = e i g ( f u l l (K(int_dofs , int_dofs )),
f u l l (M(int_dofs , int_dofs )));

IV. Applications
While IGA is already wide spread in the mechanical engineering
community, real-world applications in the context of electrical
engineering are still rare. This is particularly true for three-

Fig. 6. Spyplots of mass matrices for different degrees. On the left the FEM
lowest order case (in this case IGA and FEM coincide). The ratio between
non-zero elements and total number of elements in the matrix is 0.0045. On
the right C 2 B-spline basis functions of degree p = 3. In this case the ratio is
0.0199.
Fig. 5. Lorentz detuning in a TESLA cavity cell: In green the design geometry,
in color the displaced walls. The magnitude of the displacement u is enhanced
by a factor 5 · 105 for visibility.

dimensional problems that require the more complicated curland div-conforming spline spaces from equations (13) and (14).
In [30] the two-dimensional shape optimization of a magnetic
density separator was proposed, and in [31] the optimization
of ferromagnetic materials in a magnetic actuator was shown.
Recently, researchers have also started to investigate Isogeometric Analysis of integral equations in the context of electrical
engineering, e.g. [32, 33, 34, 35].
In the following subsections the application of IGA is demonstrated for several real-world examples in two and three dimensions and an outlook on isogeometric Boundary Element Methods (BEMs) is given. The aim is not to give a precise description
but rather to illustrate that IGA is indeed a useful tool for the
Compumag community.

A. Radio Frequency Cavities
To achieve acceleration of the particle bunches in particle accelerator Radio Frequency (RF) cavities, the electromagnetic field
has to oscillate at a very specific frequency, synchronously to
the movement of the charges. The eigenfrequency is determined
by the shape of the cavity walls, which is therefore critical for
the design of any cavity. However, the high-energy field exerts
a radiation pressure on the walls, which impresses a mechanical
deformation of the domain. Albeit small, this deformation may
lead to a significant shift of the resonance frequency. This effect is known as Lorentz detuning [36, 37, 38] and needs to be
predicted with high precision in order to achieve a robust cavity
design.
Applying standard FEMs present two main problems: the domain boundary is approximated by polynomials and the deformation of the cavity walls may require an interpolation and
remeshing step or an ad-hoc mesh movement procedure. In [39]
the MpCCI Multiphysics Interface [40] was used for exchanging
geometry and solution between the CST Microwave predecessor
MAFIA based on the Finite Integration Technique and the software package ParaFep based on Finite Elements [41, 42]. IGA
is able to overcome these issues allowing an exact representation of the geometry, leading to higher accuracy, and a direct
application of the computed deformation to the design shape,
without any further approximation. Finally it offers the possibility to obtain highly smooth solutions, which can prove extremely
valuable for particle tracking applications [43].

Maxwell’s Eigenvalue Problem. As a first example we consider the academic case of the computation of the eigenmodes
in a cylindrical resonating cavity (pill-box), where the eigenmodes can be computed analytically. The fields are assumed to
be time-harmonic and oscillating in vacuum ( = 0 , µ = µ0 ),
with no charges or currents. The first order system (1) can then
be rewritten as a second order equation for the electric field E
only
∇ × ∇ × E = µ0 0 ω 2 E.

(18)

The solution of problem (18) is a sequence of eigenmodes
2
, Em ) which represents the excitable modes in the cavity.
(ωm
The quantities ωm are the resonance frequencies fm = ωm /2π.
Equation (18) can be discretized using B-splines belonging to the
space (17) in order to obtain a generalized eigenvalue problem
Ke = ω 2 Me,

(19)

with K and M the curl-curl and mass matrix respectively, and
e the vector containing the electric field degrees of freedom.
For the software implementation one can follow the lines of
Listing 1. The matrices obtained with the IGA discretization
typically present a larger bandwidth compared to their FEM
counterparts (see Fig. 6), but, as previously mentioned, they
are also smaller for a given accuracy. We present a comparison between the two approaches in terms of efficiency of the
solution of the eigenvalue problem. A set of matrices, obtained from meshes with increasing refinement, was generated
for order 2 and 3 polynomial FEM basis functions using the
proprietary software CST [41], and exported to MATLAB. The
same Arnoldi/Lanczos solver is used for the IGA matrices and
for the FEM ones. In Fig. 7 we report the computational time
necessary for attaining a prescribed level of accuracy. Given the
higher accuracy-per-degree-of-freedom of IGA, it is possible to
achieve a considerable speed-up [44].

Lorentz Detuning. The procedure we propose is straightforward. First, we compute the fields in the cavity by solving (18).
As a second step we solve the mechanical problem for the cavity
walls. Given the small deformations involved we use the linear
elasticity model


∇ · 2η∇(S) u + λI∇ · u = pn

(20)

for the deformation u, with η, λ the Lamé constants of niobium
and pn the pressure applied. The symbol ∇(S) denotes the
symmetric gradient operator. The electromagnetic problem (18)

|Ez |[V m−1 ]

1,000

t[s]

·108

FEM - p = 2
IGA - p = 2
FEM - p = 3
IGA - p = 3

500

4

2

0

−0.4

0
10−11 10−10 10−9 10−8 10−7 10−6 10−5
Relative Error

where Epk and Hpk are the field peak values, nc is the outside normal to the cavity and ∗ denotes the complex conjugate
operator.
As mentioned above the computed deformation u can be directly
applied to the control polygon of the cavity geometry to obtain
the deformed cavity. The solution of Maxwell’s eigenvalue
problem on this domain gives the frequency shift.
For the computation of the frequency shift in the case of a pillbox cavity and a comparison with a sensitivity analysis procedure
we refer the interested reader to [44].
In Fig. 5 we depict the deformation for the more realistic case
of a 1.3 GHz TESLA cavity [45]. This deformation is typically
of the order of tenths of nm, nevertheless, this leads to a measurable frequency shift in the order of hundredths of Hz that
needs to be addressed during operation. The computed shift is
approximately 1 kHz [46] which is in good agreement with the
literature.
Field Flatness Tuning. When dealing with multi-cell cavities
such as the TESLA design, small variations between the cells’
shape are sufficient to substantially alter the field profile. Of
particular interest for the correct operation of RF structures is
to achieve a uniform energy distribution in each cell. This
presents two advantages: it maximizes the accelerating voltage
of the cavity (i.e. the net energy gained by the particles) and
it minimizes the peak surface fields, which are responsible for
electric field emission or quench in the superconducting walls
[47].
We set z as the longitudinal axis of the cavity and we denote by
Epk,j the peak value of Ez (r = 0, z) in the j-th cell. The field
flatness is typically measured by two quantities:

4
|Ez |[V m−1 ]

couples into the mechanical problem by the radiation pressure


1 
prad = − 0 Epk · nc E∗pk · nc
(21)
4


1 
+ µ0 Hpk × nc · H∗pk × nc
4

1 − (maxj |Epk,j | − minj |Epk,j |)
E (|Epk,j |)
std (Epk,j )
η2 = 1 −
,
E (|Epk,j |)

0

0.2

0.4

0.2

0.4

z[m]

8

Fig. 7. Computational time required to solve the eigenvalue problem with
ARPACK for finding the first accelerating mode in the pill-box cavity, with
a prescribed accuracy. The IGA implementation was performed in GeoPDEs
while for the FEM simulation CST EM STUDIO [41] was used.

η1 =

−0.2

·10

3
2
1
0

−0.4

−0.2

0
z[m]

Fig. 8. Absolute value of the longitudinal electric field Ez in the untuned
(top) and tuned (bottom) TESLA cavity. The computed value for the field
flatness improve from η1 = 76.81 % and η2 = 92.39 % to η1 = 97.76 % and
η2 = 99.15 %.

In practice, the tuning is performed through mechanical deformation. The field flatness is measured and, using a circuit model
for the cells as capacitively coupled LC oscillators, the required
frequency shift for each cell is computed [47]. To lower (resp.
increase) the frequency, each cell is shortened (elongated) along
the z axis by a tuning machine which clamps the cavity between
cells and applies forces to obtain a permanent deformation [48].
Numerically, the tuning is performed through a multi-objective
shape optimization procedure that aims at maximizing the field
flatness parameters (22) and achieving a resonance frequency
of 1.3 GHz. The geometry parameters that are typically chosen
for the optimization are the length of the two end-cups half-cells
and the equatorial radius of the cavity. The first two parameters strongly influence the field flatness, while the equatorial
radius is mainly responsible for fixing the frequency of the cavity. The three parameters are used for a non-linear constrained
optimization procedure using the SQP method. The bounds for
the parameters are set to ±1.5 mm.
In Fig. 8 the longitudinal electric field along the cavity axis
before and after the tuning is depicted.

The advantage of IGA when dealing with shape optimization
is the possibility to affect both the geometry and the mesh at
the same time by simply moving the control points. For each
configuration there is no need of a remeshing step that might
introduce undesired noise in the computation. Furthermore,
NURBS geometries are well suited for the computation of shape
derivatives in order to obtain gradient information.

(22a)
B. Electric Machine Simulation
(22b)

where E and std denote the expected value and the standard
deviation operators. Both η1 and η2 are typically required to be
≥ 0.95 for a well tuned cavity.

For the modeling of electric machines, a common approach is
the magnetostatic formulation where the eddy currents and the
displacement currents are neglected, see e.g. [18]. Under these
assumptions, by introducing the magnetic vector potential A,

FEM
IGA

30

EMF [V]

20

10

0
0

5

10
15
20
Harmonic Order

25

30

Fig. 10. Spectrum of the first 32 modes of the EMF of the PMSM, cf. [49, 51].
Fig. 9. Computational domain Ω of a PMSM model as discussed in [49, 51].

0.2
Γ

with B = ∇ × A, one obtains
∇ × (ν∇ × A) = J

(23)

on the domain Ω with appropriate boundary and gauging conditions. Furthermore, given the motor structure, it is often
sufficient to solve the problem only on a 2D cross section.
One obtains a Poisson problem for the longitudinal component
A = [0, 0, u]>
−∇ · (ν∇u) = Jz ,

-0.2

Ωcon
Ωair

0.2

-0.2
Fig. 11. Quadrupole as given in Figure 7.1 of [54] on the left. The red pole tips
are modeled as interfaces subject to uncertainty. On the right, discretization in
terms of multipatch NURBS. All units in meter. Image based on [57].

(24)

where ν is the space-dependent reluctivity, and Jz represents the
current excitations due to the presence of coils and/or permanent
magnets.
In Fig. 9 the geometry of a Permanent Magnet Synchronous
Machine (PMSM) is depicted. Although only two dimensional,
the topology and the presence of different materials requires the
construction of a high number of patches (12 for the rotor, 78 for
the stator), which is still low compared to the number of elements
in a FEM discretization. In [49] we propose a harmonic statorrotor coupling method to overcome this issue and to allow for
easy handling of the machine rotation.
The application of IGA to machine simulation is particularly
interesting for the possibility to exactly represent its circular
shape and, even more so, for the smoothness of the computed
fields. The evaluation of torques and electromotive force (EMF)
are often calculated from the fields in the air gap using the
Maxwell’s stress tensor. Since the results obtained through
this approach are very sensitive to the representation and to the
discretization of the air gap [50], the higher continuity of the
isogeometric solutions has been proved beneficial.
In Fig. 10 the spectrum of the first 32 harmonics of the EMF for
the PMSM is depicted. The results show good agreement with
the ones obtained with a classical FEM simulation, but the IGA
system is considerably smaller [49, 51].
One further interesting possibility for electric machine simulation taking into account the machine rotation would be a coupling strategy using classical IGA on the stator and the rotor
and an IGA-BEM discretization in the air gap region (see e.g.
[52]). In sub-section D we briefly introduce the BEM setting
and show that it can be readily applied in conjunction with the
IGA concept.

C. Accelerator Magnets
In particle accelerators normal and superconducting magnets
are used for focusing and bending the particle beams. The
simulation of these devices is an important task during the computer aided design process as these devices are operated at their
physical limits but also as to deliver highly accurate field maps
for beam dynamics simulations and, finally, the simulation of
quench protection [58, 59].
Uncertainty Quantification. The first magnet example is devoted to illustrate the treatment of shape uncertainties, see e.g.
[53], in the context of magnet design. On the left side of Fig. 11
the two-dimensional geometry of a model quadrupole magnet
adopted from Fig. 7.1 of [54] is depicted, where the pole tips are
assumed to be affected by uncertainty, e.g. due to manufacturing
imperfections. The initial pole shape is described by hyperbolas,
i.e., through the relation x/a2 − y/b2 = 1 in local coordinates,
where for the initial shape we have set a = 0.05 and b = 0.056,
respectively. In a stochastic setting, uncertain interfaces can be
modeled by choosing a or b as random variables, or equivalently
by using any other CAD standard shape representation with random parameters. Then, a truncated Karhunen-Loève expansion
can be used to obtain a reduced number of uncorrelated random inputs [55]. A multi-patch configuration, consisting of 36
patches is chosen to represent the geometry as depicted in Fig.
11 on the right side. Each patch is described by a mapping
of the type (9) using second degree NURBS basis functions
with C 1 -continuity. In particular, this allows for an exact representation of hyperbolic shapes. Focusing on shape variations,
we model the material to be linear with a constant permeability µ = 4π10−2 H m−1 . A total piecewise constant current of

Table I. Quadrupole convergence at different mesh levels. Error
estimate of the shape gradient and finite difference approximation for
different number of perturbed control points. The reference for g is
computed with 367 641 DOF.

Parameter
∆h g /g0
∆FD δg

level 1

level 2

level 3

level 4

0.56%

0.05%

< 0.01%

< 0.01%

5.85 × 10−5

3.96 × 10−5

7.29 × 10−5

1.38 × 10−4

15 MA is supplied for each of the four conductor parts of Fig.
11. Homogeneous Dirichlet boundary conditions are applied on
the whole boundary Γ. Following the iso-parametric concept,
a basis for the discrete subspace of H 1 is constructed in the
same space as the mapping F. The multipole coefficients are
evaluated at a reference radius of r0 = 20 mm.
We focus on the uncertainty in the quadrupole gradient, defined
as g := 2B2 /r02 where Bn is the n-th normal multipole coefficient. As we are concerned with several parameters and only
one cost function, adjoint techniques are well suited to this end
[56, 57]. Here, no assumptions, despite the C 1 -smoothness, are
made for the shape uncertainty and, therefore, we resort to a
worst-case analysis. The maximum deviation of g from its design value g0 = 19.73 T m−2 is investigated for different levels
of shape parametrization, characterized by the number of free
control points. On the coarsest level the hyperbola is described
by three control points, however, the end points are kept fixed
to avoid variability in the singular points, as this would require
a more general shape calculus as presented here. Hence, only
one control point per pole is subject to uncertainty. Through
mesh refinement this number is increased by one on each level,
up to level four. After refining the geometry, in each direction
every mesh cell is divided by twenty. For an evaluation of the
accuracy of the numerical approximation of the quadrupole gradient g, see Table I. There, the error on each level is estimated
with respect to a fine discretization consisting of 367 641 total
DOF, denoted as ∆h g. Additionally, the numerical computation of the shape gradient is verified. To this end in Table I the
maximum deviation with respect to a finite difference gradient
computation, denoted ∆FD , is given. Both estimated errors are
found to be sufficiently small. We emphasize that no correlation
is imposed here. If knowledge of the shape perturbations were
available, e.g., in terms of measurements, the correlation structure could be incorporated by means of convex constraints in a
worst-case scenario context as outlined in [60, p.10]. In Table
II numerical results for the different parametrization levels are
presented. Not surprisingly, a significantly smaller worst-case
estimate is obtained for the coarse parametrization. In this case,
large perturbations in the control point are necessary to obtain a
comparable shape perturbation to the finer parametrization levels. We compute the worst-case scenario by a Taylor expansion
wcsL (g) and by directly solving the optimization problem, denoted wcs∗ (g). Here, for the latter case MATLAB’s fmincon
routine is used to carry out sequential quadratic programming.
We observe a difference of about 15%, and infer that the problem
is rather sensitive to shape perturbations and first order approximations should be used to obtain rough estimates of the output
uncertainties, solely.
Shape Optimization. Another application in which IGA has
proven advantageous is the optimization of a Stern-Gerlach magnet [61]. A Stern-Gerlach magnet (see Fig. 12) is used to magnetically separate a beam of atoms or atom clusters in order to
experimentally determine their angular momentum and its spatial quantization. For this purpose, the atoms are accelerated and
shot through the aperture of the magnet. To deflect the particles,

Table II. Quadrupole error by linearization and error estimated by
means of MATLAB’s fmincon. wcsL (g) indicates the worst-case
scenario computed by Taylor expansion, wcs∗ (g) the one obtained by
solving the optimization problem. Perturbation magnitude of s = 0.1.

num. free CP
4
8
12

wcsL (g) / g0
7.82
16.77
18.57

wcs∗ (g) / g0
8.48
19.77
21.77

rel. diff.
7.78%
15.17%
14.70%

Fig. 12. 3D model of one half of the Stern-Gerlach magnet with Rabi-type pole
tips (modeled with CST EM STUDIO [41]), see [61].

the magnet should provide both, a high magnetic field strength
to cause a precession of the magnetic dipoles, and a high magnetic field gradient to deflect the atoms. To obtain an acceptable
resolution, an additional requirement is a high homogeneity of
the magnetic field gradient in the beam area. These quantities
are strongly influenced by the geometry of the pole tips. In this
example, Rabi-type pole tips are considered. The optimization
process consists of finding an optimal geometry for these pole
tips to satisfy the before-mentioned requirements.
The Stern-Gerlach magnet is operated with a DC current. Thus
a 3D non-linear magnetostatic formulation of the Maxwell’s
equations is suited to calculate the magnetic field
∇ × (ν(B)∇ × A) = J,

(25)

where ν(B) is the non-linear reluctivity and J the exciting current density. Due to symmetry only one half of the magnet is
considered in the optimization process, as depicted in Fig. 12.
To obtain an efficient simulation, the magnet is modeled by using a combination of a magnetic equivalent circuit and a field
model discretized by IGA: the outer yoke and coils are modeled
by the magnetic equivalent circuit which is extracted from a full
3D simulation.
Finally, the area of the pole tips is discretized by a 2D IGA
model, i.e. equation (24) on the domain Ωp as shown in Fig. 13.
This enables a smooth representation of the pole shapes by using
NURBS and simplifies the optimization process, as the deformation of the poles is easily achieved by shifting control points
and adjusting weights. The control points used for the NURBS
representation of the poles are depicted in Fig. 14. The partial
model is divided into 3 patches, namely the left pole region, gap
and right pole region. Both the magnetic equivalent circuit and
the partial model are connected by a field-circuit coupling. To
avoid a non-linear evaluation of this coupled model and thus save
computational effort, a further simplification is employed: the
permeability in the IGA partial model is frozen. This is possible
as the variations in the permeability distribution due to changes
in the geometry are small.
For the optimization, the quantities of interest in the beam area

8
Original geometry
Original control points
Optimized geometry

7

Coil

y [mm]

6
Iron

5
4
3
2

Coil

1
0
Partial model

Air gap

-6

-4

-2

x
Symmetry plane (Prescribed A=0)

Ωbeam of the magnet are the average magnetic field gradient and
the inhomogeneity which are given by
Z
1
τav =
τ (x, y) dΩ
(26)
|Ωbeam |
Ωbeam

1
|Ωbeam |

Z

Ωbeam



2
τ (x, y)
−1
dΩ ,
τav

4

(BEM) gained a lot of attraction in research. The idea of BEM
is the representation of the solution u of a PDE through a representation formula, which expresses the solution via functions
on the boundary of the problem domain only.
There exist different such representations, each with their specific strengths and drawbacks. As a simple example, the so
called indirect representation formula for the Laplace Dirichlet
problem
in Ω or ΩC

−∆u = 0
u|Γ = g

(27)

on Γ,

(30)

is given by

d|B|
dx

where τ (x, y) =
is the magnetic field gradient. The goal
function used for the optimization process is
τw
τw
+−
,
f (x, y, w) =
|τav |
|τav |

2

Fig. 14. Optimized and original geometry of the pole tips including beam area
(gray) and labeled control points of original geometry, see also [61].

Beam area Ωbeam

Fig. 13. 2D schematic of one half of the Stern-Gerlach magnet with Rabi-type
pole tips, see [61].

v
u
u
=t

0

x [mm]

Ωp

y

and

-8

(28)

where τw is a free parameter, and τav ,  are the quantities of
interest defined above. As the requirements of high average
magnetic field gradient and homogeneity are, to a certain degree,
antithetic, the free parameter τw controls the preference of one
over the other quantity. The overall optimization problem can
be written as
min f (x, y, w),
(29)
x,y,w

where x, y, w are the coordinates and weights of the control
points defining the poles. These quantities are subject to certain
geometrical limits to ensure the validity of the geometry.
For the simulation, the model is implemented using GeoPDEs.
B-splines of order 5 are used to represent the magnetic vector
potential A. This way, the magnetic flux density B and even
the average magnetic field gradient τav are smooth functions
opposed to conventional FE simulations.
The original and optimized geometry of the pole shoes is depicted in Fig. 14. The improvements of the average magnetic
field gradient and homogeneity after the optimization are shown
in Table III. For validation, the optimized geometry is imported
into a 3D model of the magnet in the commercial software
CST EM STUDIO [41]. Note that the presented optimization
approach using IGA and field-circuit coupling is considerably
faster than the commercial software, while offering results which
are in very good agreement [61].
D. IGA and the Boundary Element Method
The framework of Isogeometric Analysis is not limited to
FEM. In recent years, Isogeometric Boundary Element Methods

e
:=
u(x) = V(µ)(x)

Z

Γ

1
µ(y) dΓy ,
4π|x − y|

(31)

on R3 \Γ, where µ denotes some unknown density on the boundary. This density can be found by solving the arising variational
problem: find µ in X(Γ) such that
hV(µ), νiΓ = hg, νiΓ ,

∀ν ∈ X(Γ),

(32)

where X(Γ) denotes an appropriate function space on the
e to the boundary. Choosing
boundary and V is the restriction of V
a suitable 2D spline space S(Γ), which can be defined patchwise
as the 2D analogue of (15) mapped in the physical domain as
in (17). This leads to the discrete formulation of finding a
µh ∈ S(Γ) such that
hV(µh ), νh iΓ = hg, νh iΓ ,

∀νh ∈ S(Γ).

(33)

After inserting the canonical basis of spline space S(Γ) we receive the linear system
Vw = g

(34)

where g encodes the discretization of the Dirichlet data given in
(30), and w encodes the desired density function. Plugging the
discrete density back into (31), one can evaluate the solution u
of (30) in any given point x ∈
/ Γ, c.f. Fig. 15.

This representation of the solution counteracts one of the greatest
weaknesses of IGA FEM: where in many cases trivariate spline
mappings need to be constructed by hand, a boundary element
method can operate with a boundary representation (BRep) only.
Moreover, it decreases the dimension of the problem: where
FEM requires volumetric elements, BEM operates on the surface
only, thus the elements are 2D elements in the sense of the
reference domain.

Table III. Average magnetic field gradient and inhomogeneity factor
before and after the optimization.

CST (3D, original geometry)
CST (3D, optimized geometry)
Improvement
GeoPDEs (2D, original geometry)
GeoPDEs (2D, optimized geometry)
Improvement

τav
−237 T/m
−266 T/m
12.2%
−240 T/m
−282 T/m
17.5%


0.0503
0.0201
60.0%
0.0477
0.0122
74.4%

It should be noted that the indirect representations are entirely
insensitive to the choice of interior and exterior domain. Since
the BEM operates on the boundary of the chosen domain only,
this makes a boundary element approach an excellent choice to
solve exterior problems, since neither bounding boxes nor DOF
intensive exterior discretizations are required.
However, a BEM approach poses different challenges. Due to
the global nature of the integral term the matrices arising from
(33) are densely populated. This turns out to be not as bad
as it sounds: due to the decreased dimensionality fewer degrees of freedom are required, thus the linear systems become
much smaller than in the case of FEM. Due to the approach
via B-splines, this effect becomes even more notable. Even
for real-world problems the arising linear systems are almost of
sizes such that they can be solved in a dense representation by a
desktop computer. Moreover, efficient compression techniques
to deal with the discrete linear systems are available. Such techniques are well understood and ready for industrial applications.
An introduction to and comparison of compression methods can
be found in [62].
An introduction into the realm of isogeometric BEM in the
context of the Helmholtz equation which elaborates on everything mentioned above and explains the strength and drawbacks
of the combination of IGA and BEM is given by [63]; a first
computational approach for the Maxwell case was also recently
proposed [64].
V. Conclusions
Isogeometric Analysis can be interpreted as a Finite Element
Method that generalizes the set of basis functions from polynomials to B-splines and rational B-splines. This choice guarantees
several advantages, from the exact parametrization of geometries
defined via CAD to a higher accuracy per-degrees-of-freedom.
It also allows for solution fields with higher smoothness. These
reasons have made IGA a successful topic in recent years, particularly in the setting of mechanics and fluid simulation.
In this review article we show the benefits of IGA in the context of computational electromagnetics by presenting several 2and 3-dimensional simulation examples. Its application in the
context of shape optimization and uncertainty quantification is
of particular interest given the possibility of deforming the domain and the mesh at the same time by means of a few control
points. Finally, IGA shows great potential when combined with
the boundary element approach since it eliminates the necessity
of constructing trivariate parameterizations.

Fig. 15. Boundary mesh and potential evaluation according to (31) in the interior
(Laplace problem), see also [63].

lence Initiative’ of the German Federal and State Governments
and the Graduate School of CE at TU Darmstadt.
References
[1] M. Clemens, S. Schöps, C. Cimala, N. Gödel, S. Runke, and
D. Schmidthäusler. "Aspects of Coupled Problems in Computational Electromagnetics Formulations", ICS Newsletter,
3, 2012.
[2] T.J.R. Hughes, J.A. Cottrell and Y. Bazilevs, "Isogeometric
analysis: CAD, finite elements, NURBS, exact geometry
and mesh refinement", Computer Methods in Applied Mechanics and Engineering, 194, 4135-4195, 2005.
[3] J.A. Cottrell, T.J.R. Hughes and Y. Bazilevs, "Isogeometric
Analysis: Toward Integration of CAD and FEA", John Wiley
& Sons Ltd, 2009.
[4] H. Al Akhras, T. Elguedj, A. Gravouil and M. Rochette,
"Isogeometric analysis-suitable trivariate NURBS models
from standard B-Rep models", Computer Methods in Applied Mechanics and Engineering, 307, 256-274, 2016.
[5] P. Monk, "Finite Element Methods for Maxwell’s Equations", Oxford: Oxford University Press, 2003.
[6] J.-M. Jin, "The Finite Element Method in Electromagnetics,
3rd Edition", Wiley-IEEE Press, 2014,
[7] W.A. Wall, M.A. Frenzel and C. Cyron, "Isogeometric structural shape optimization", Computer Methods in Applied
Mechanics and Engineering, 197(33-40), 2976-2988, 2008.
[8] X. Qian, "Full analytical sensitivities in NURBS based isogeometric shape optimization", Computer Methods in Applied Mechanics and Engineering, 199(29–32), 2059-2071,
2010.
[9] T.J.R. Hughes, A. Reali and G. Sangalli, "Duality and unified
analysis of discrete approximations in structural dynamics
and wave propagation: Comparison of p-method finite elements with k-method NURBS", Computer Methods in Applied Mechanics and Engineering, 197(49-50), 4104-4124,
2008.

Acknowledgment
This work is supported by the German BMBF in the context
of the SIMUROM project (grant nr. 05M13RDA), by the DFG
(grant nr. SCHO1562/3-1 and KU1553/4-1) and by the ’Excel-

[10] A. Bartezzaghi, L. Dedè and A. Quarteroni, "Isogeometric Analysis of high order Partial Differential Equations on
surfaces", Computer Methods in Applied Mechanics and
Engineering, 295, 446-469, 2015.

[11] B.S. Hosseini, M. Möller and S. Turek, "Isogeometric
Analysis of the Navier–Stokes equations with Taylor–Hood
B-spline elements", Applied Mathematics and Computation,
267, 264-281, 2015.

[28] A.F. Sarmiento, A.M.A. Cortes, D.A. Garcia, L. Dalcin, N. Collier and V.M. Calo, "PetIGA-MF: a multifield high-performance toolbox for structure-preserving Bsplines spaces", arXiv:1602.08727, 2016.

[12] A. Buffa, J. Rivas, G. Sangalli and R. Vázquez, "Isogeometric Discrete Differential Forms in Three Dimensions",
SIAM Journal on Numerical Analysis, 49(2), 818-844, 2011.

[29] V.P. Nguyen, C. Anitescu, S.P. Bordas and T. Rabczuk,
"Isogeometric analysis: An overview and computer implementation aspects", Mathematics and Computers in Simulation, 117, 89-116, 2015.

[13] A. Buffa, G. Sangalli and R Vázquez, "Isogeometric analysis in electromagnetics: B-splines approximation", Computer Methods in Applied Mechanics and Engineering, 199,
1143-1152, 2010.
[14] A. Buffa, G. Sangalli and R. Vázquez. "Isogeometric methods for computational electromagnetics: B-spline and Tspline discretizations", Journal of Computational Physics,
257, Part B, 1291-1320, 2013.
[15] C. Blanc, and C. Schlick, "Accurate parametrization of
conics by NURBS", IEEE Computer Graphics and Applications, 16(6), 64-71, 1996.
[16] L. Piegl and W. Tiller, "The NURBS Book", 2nd ed.
Springer, 1997.
[17] C. de Boor, "A Practical Guide to Splines", Applied Mathematical Sciences, 27, Springer, 2001.
[18] S.J. Salon, "Finite element analysis of electrical machines",
Boston USA: Kluwer academic publishers, 1995.
[19] C. de Falco, A. Reali, and R. Vázquez, "GeoPDEs: a
research tool for isogeometric analysis of PDEs", Advances
in Engineering Software, 42(12), 1020-1034, 2011.
[20] R. Vázquez, "A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0",
Computers & Mathematics with Applications, 72(3), 523554, 2016.
[21] Mathworks, MATLAB Getting Started Guide, 2009.
[22] J.W. Eaton, D. Bateman, S. Hauberg and R. Wehbring,
"The GNU Octave 4.0 Reference Manual 1/2: Free Your
Numbers", Samurai Media Limited, 2015.
[23] J. Ahrens and A. Henderson, The ParaView Guide, Kitware
Inc., 2015.
[24] E.M. Garau and R. Vázquez, "Algorithms for the implementation of adaptive isogeometric methods using hierarchical B-splines", to appear in Applied Numerical Mathematics, DOI: 10.1016/j.apnum.2017.08.006
[25] M.S. Pauletti, M. Martinelli, N. Cavallini and P. Antolin,
"Igatools: An isogeometric analysis library", SIAM Journal
on Scientific Computing, 37(4), C465-C496, 2015.
[26] B. Juettler, U. Langer, A. Mantzaflaris, S. Moore and W.
Zulehner, "Geometry + simulation modules: Implementing
isogeometric analysis", Proceedings in Applied Mathematics and Mechanics, 14 (1), 961-962, special Issue: 85th Annual Meeting of the Int. Assoc. of Appl. Math. and Mech.
(GAMM), Erlangen, 2014.
[27] L. Dalcin, N. Collier, P. Vignal, A.M.A. Côrtes and V.M.
Calo, "PetIGA: A framework for high-performance isogeometric analysis", Computer Methods in Applied Mechanics
and Engineering, 308, 151-181, 2016.

[30] N.D. Manh, A. Evgrafov, J. Gravesen and D. Lahaye, "Isogeometric shape optimization of magnetic density separators", COMPEL: The International Journal for Computation
and Mathematics in Electrical and Electronic Engineering,
33(4), 1416-1433, 2014.
[31] S.W. Lee, J. Lee and S. Cho, "Isogeometric Shape Optimization of Ferromagnetic Materials in Magnetic Actuators,
IEEE Transactions on Magnetics, 52(2), 1-8, 2016.
[32] T. Sekine, K. T. Miura and H. Asai, "Method of moments
based on isogeometric analysis for electrostatic field simulations of curved multiconductor transmission lines" 2014
IEEE 23rd Conference on Electrical Performance of Electronic Packaging and Systems, Portland, OR, 195-198, 2014.
[33] T. Sekine, K. T. Miura and H. Asai, "Electrostatic field
simulations of curved conductors by using method of moments based on isogeometric analysis", 2014 International
Symposium on Electromagnetic Compatibility, Gothenburg,
192-195, 2014.
[34] J. Li, D. Dault, R. Zhao, B. Liu, Y. Tong and B. Shanker,
"Isogeometric analysis of integral equations using subdivision", 2015 IEEE Antennas and Propagation Society International Symposium, APS 2015 - Proceedings, Institute of
Electrical and Electronics Engineers Inc., 153-154, 2015.
[35] R.N. Simpson, S.P.A. Bordas, H. Lian and J.
Trevelyan, "An Isogeometric Boundary Element Method
for elastostatic analysis: 2D implementation aspects",
arXiv:1302.5305, 2013.
[36] G. Devanz, M. Luong and A. Mosnier, "Numerical simulations of dynamic Lorentz detuning of SC cavities", EPAC
(2002): Proceedings of the 8th European Particle Accelerator Conference, Paris, France, 2002.
[37] H. Gassot, "Mechanical Stability of the RF superconducting Cavities", EPAC (2002): Proceedings of the 8th European Particle Accelerator Conference, Paris, France, 2002.
[38] E. Zaplatin et al., "Lorentz Force Detuning Analysis for
Low-Loss, Reentrant, and Half-Reentrant Superconducting
RF Cavities", Proceedings of LINAC 2006, Knoxville, Tennessee USA, 2006.
[39] U. Schreiber, and U. van Rienen, "Coupled Calculation
of Electromagnetic Fields and Mechanical Deformation",
Scientific Computing in Electrical Engineering, Springer,
Berlin, Heidelberg, 63-68, 2006.
[40] MpCCI - Mesh-based parallel Code Coupling Interface.
Version 1.3, PALLAS GmbH, Hermülheimer Str. 10, D50321 Brühl, Version 1.3, 2002.
[41] CST STUDIO SUITE, Computer Simulation Technology
AG, 2016.

[42] R. Niekamp, and E. Stein, "An object oriented approach
for parallel 2– and 3–dimensional adaptive nonlinear FiniteElement-computations", International Journal of Computer
and Structures, 80, 317-328, 2002.
[43] E. Gjonaj, W. Ackermann, T. Lau, T. Weiland, and M.
Dohlus, "Coupler kicks in the third harmonic module for
the XFEL", Particle accelerator, Proceedings, 23rd Conference, 9, 4-8, 2009.
[44] J. Corno, C. de Falco, H. De Gersem, S. Schöps, "Isogeometric Simulation of Lorentz Detuning in Superconducting
Accelerator Cavities", Computer Physics Communications,
201, 1-7, 2016.
[45] B. Aune et al., "Superconducting TESLA cavities", Physical Review Special Topics-Accelerators and Beams, 3(9),
2000.
[46] J. Corno, "Numerical Methods for the Estimation of the
Impact of Geometric Uncertainties on the Performance of
Electromagnetic Devices", PhD Thesis, submitted, 2017.
[47] H. Padamsee, J. Knobloch and T. Hays, "RF Superconductivity for Accelerators", Wiley, 2008.
[48] G.R. Kreps, J. Sekutowicz and D. Proch, "Tuning of the
TESLA Superconducting Cavities and the Measurement of
Higher Order Mode Damping", Proc. XVth Conference on
Charged Particle Accelerators, Dubna, Russia. 1996.
[49] Z. Bontinck, J. Corno, H. De Gersem, S. Schöps, "Isogeometric Analysis and Harmonic Stator-Rotor Coupling for
Simulating Electric Machines", Computer Methods in Applied Mechanics and Engineering, submitted, 2017.
[50] D. Howe and Z.Q. Zhu, "The influence of finite element discretisation on the prediction of cogging torque in permanent
magnet excited motors", IEEE Transactions on Magnetics,
28(2), 1080-1083, 1992.

[58] T. Roggen, H. De Gersem, B. Masschaele, W. Ackermann,
S. Franke and T. Weiland, "Embedding finite element results
for accelerator components in a moment approach beam
dynamics code", International Particle Accelerator Conference (IPAC 2011), 2011.
[59] I. Cortes Garcia, S. Schöps, L. Bortot, M. Maciejewski,
M. Prioli, A. Fernandez Navarro, B. Auchmann and A. Verweij, "Optimized Field/Circuit Coupling for the Simulation
of Quenches in Superconducting Magnets", IEEE Journal
on Multiscale and Multiphysics Computational Techniques,
2(1), 97-104, 2017.
[60] I. Babuška, F. Nobile, and R. Tempone, "Worst case scenario analysis for elliptic problems with uncertainty", Numerische Mathematik, 101(2), 185-219, 2005.
[61] A. Pels, Z. Bontinck, J. Corno, H. De Gersem and S.
Schöps, "Optimization of a Stern-Gerlach Magnet by Magnetic Field-Circuit Coupling and Isogeometric Analysis",
IEEE Transactions on Magnetics, 51(12), 2015.
[62] H. Harbrecht and M. Peters, "Comparison of fast boundary
element methods on parametric surfaces", Computer Methods in Applied Mechanics and Engineering, 261, 39-55,
2013.
[63] J. Dölz, H. Harbrecht, S. Kurz, S. Schöps and F. Wolf, "A
Fast Isogeometric BEM for the Three Dimensional Laplaceand Helmholtz Problems", submitted, arXiv:1708.09162,
2017.
[64] R. N. Simpson, Z. Liu, R. Vázquez and J.A. Evans, "An
isogeometric boundary element method for electromagnetic scattering with compatible B-spline discretizations",
arXiv:1704.07128, 2017.
Authors Name and Affiliation

[51] P. Bhat, Z. Bontinck, J. Corno, H. De Gersem, and S.
Schöps, "Modeling of a Permanent Magnet Synchronous
Machine Using Isogeometric Analysis", 18th International
Symposium on Electromagnetic Fields in Mechatronics,
Electrical and Electronic Engineering (ISEF 2017), 2017.

Zeger Bontinck, Jacopo Corno, Herbert De Gersem, Stefan
Kurz, Andreas Pels, Sebastian Schöps, Felix Wolf from the
Institut für Theorie Elektromagnetischer Felder and the Graduate School CE, Technische Universität Darmstadt, Darmstadt,
Germany.

[52] S. Kurz, J. Fetzer, G. Lehner and W.M. Rucker, "A novel
formulation for 3D eddy current problems with moving bodies using a Lagrangian description and BEM-FEM coupling", IEEE Transactions on Magnetics, 34(5), 3068-3073,
1998.

Carlo de Falco from MOX – Modellistica e Calcolo Scientifico,
Dipartimento di Matematica, Politecnico di Milano, Milano,
Italy.

[53] S. Clénet, "Uncertainty Quantification in Computational
Electromagnetics: The stochastic approach", ICS Newsletter
(International Compumag Society), 13, 3-13, 2013.

Rafael Vázquez from the Institute of Mathematics, École
Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

[54] S. Russenschuck, "Field computation for accelerator magnets: analytical and numerical methods for electromagnetic
design and optimization", John Wiley and Sons, 2011.
[55] D. Xiu, "Numerical Methods for Stochastic Computations:
A Spectral Method Approach", Princeton University Press,
2010.
[56] M. Hinze, R. Pinnau, M. Ulbrich and S. Ulbrich, "Optimization with PDE Constraints". Springer, 2008.
[57] U. Römer, "Numerical Approximation of the Magnetoquasistatic Model with Uncertainties and its Application to
Magnet Design", PhD Thesis, Technische Universität Darmstadt, Springer, 2015.

Jürgen Dölz from the Departement Mathematik und Informatik,
Universität Basel, Basel, Switzerland.

Ulrich Römer from the Department of Mechanical Engineering
of Technische Universität Braunschweig, Braunschweig,
Germany.
Corresponding author:
tu-darmstadt.de)

Sebastian Schöps (schoeps@gsc.

