Recent Advances of Isogeometric Analysis in Computational Electromagnetics

09/18/2017 ∙ by Zeger Bontinck, et al. ∙ Technische Universität Darmstadt 0

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 Non-Uniform 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.



There are no comments yet.


page 5

page 8

page 10

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 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 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 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 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 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 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 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 pde, it is possible to inherit these properties and directly use the cad design as the computational domain, without the need of a meshing step.

Figure 1: A nurbs curve in is the projection of a B-spline curve in , this allows for the definition of circles and other conic sections.

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 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.

2 Isogeometric Discretization of Maxwell’s Equations

The differential form of Maxwell’s equations is given by


with and the electric field strength and electric flux density, and the magnetic field strength and magnetic flux density, the electric charge density and the current density. The system is completed by the material relations


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

and are the electric polarization and the magnetization respectively.

Figure 2:

B-spline basis functions of degree 1 and 2 on open, uniform knot vectors (

on top and at the bottom).

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 4. 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 ) with weakly defined curl in ; this space is commonly denoted by . 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 2.2). For classical polynomial fem it is 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].

2.1 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 be the prescribed degree and


be a vector that partitions into elements (). Then, the Cox-de Boor’s formula [16] defines one-dimensional B-spline basis functions with as


with . An example of B-spline basis of degree are shown in Fig. 2. One can notice that the support of a B-spline of degree is always knot spans and, as a consequence, each -th degree function has 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.

Figure 3: Continuity can be controlled in single points (see the black circle) by knot repetition (here: ).

*nurbs of degree are defined as rational B-splines


with a weighting parameter associated with the -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 be the knot vectors, the degrees and the number of basis functions (with ), a trivariate B-spline is given by


where and is a multi-index in the set


We will denote by the multivariate B-spline space spanned by the basis functions .

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

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 one-dimensional basis functions and control points , , a curve parameterization is given by:


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 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 increases, the curve is pulled closer to the control point , 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 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 such that each one of them is the image of the reference space through a map of the type given in (9). An example of such a situation is depicted in Fig. 4 on the right.

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.

2.2 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 that can be exactly parametrized with a mapping of the type in (8), i.e.


with the reference domain . We will denote by the Jacobian of the transformation.

We define the following pull-back functions


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:


Here we have also introduced the usual Sobolev space of functions with square integrable gradient and the space of functions with weak divergence in .

To obtain a conforming discretization of , 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


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


To define the spaces in the physical domain, we use the pull-backs (10), i.e. a conforming discretization of is given by


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 , although inside each patch the discretization remains highly smooth.

3 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 4 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 2, 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 and the mass matrix 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 4. 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].

1% Load Geometry
2geo = geo_load (’geo_file.mat’);
Listing 1: A simple example of how to solve Maxwell’s eigenvalue problem in GeoPDEs

4 Applications

Figure 5: Lorentz detuning in a TESLA cavity cell: In green the design geometry, in color the displaced walls. The magnitude of the displacement is enhanced by a factor for visibility.

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-dimensional problems that require the more complicated curl- and 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 bem 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.

4.1 Radio Frequency Cavities

To achieve acceleration of the particle bunches in particle accelerator 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 fem 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 (, ), with no charges or currents. The first order system (1) can then be rewritten as a second order equation for the electric field only


The solution of problem (18) is a sequence of eigenmodes which represents the excitable modes in the cavity. The quantities are the resonance frequencies .

Figure 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 . On the right B-spline basis functions of degree . In this case the ratio is .

Equation (18) can be discretized using B-splines belonging to the space (17) in order to obtain a generalized eigenvalue problem


with and the curl-curl and mass matrix respectively, and 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].

Figure 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.

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


for the deformation , with , the Lamé constants of niobium and the pressure applied. The symbol denotes the symmetric gradient operator. The electromagnetic problem (18) couples into the mechanical problem by the radiation pressure


where and are the field peak values, is the outside normal to the cavity and denotes the complex conjugate operator.

As mentioned above the computed deformation 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 pill-box 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 TESLA cavity [45]. This deformation is typically of the order of tenths of , nevertheless, this leads to a measurable frequency shift in the order of hundredths of that needs to be addressed during operation. The computed shift is approximately [46] which is in good agreement with the literature.

Figure 8: Absolute value of the longitudinal electric field in the untuned (top) and tuned (bottom) TESLA cavity. The computed value for the field flatness improve from and to and .

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 as the longitudinal axis of the cavity and we denote by the peak value of in the -th cell. The field flatness is typically measured by two quantities:



and std denote the expected value and the standard deviation operators. Both

and are typically required to be for a well tuned cavity.

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 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 . 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 .

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.

4.2 Electric Machine Simulation

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 , with , one obtains


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


where is the space-dependent reluctivity, and represents the current excitations due to the presence of coils and/or permanent magnets.

In Fig. 9 the geometry of a 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 stator-rotor coupling method to overcome this issue and to allow for easy handling of the machine rotation.

Figure 9: Computational domain of a pmsm model as discussed in [49, 51].

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].

Figure 10: Spectrum of the first 32 modes of the EMF of the PMSM, cf. [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 4.4 we briefly introduce the bem setting and show that it can be readily applied in conjunction with the iga concept.

4.3 Accelerator Magnets

Figure 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].

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 in local coordinates, where for the initial shape we have set and , respectively. In a stochastic setting, uncertain interfaces can be modeled by choosing or

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 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 -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 . A total piecewise constant current of 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 is constructed in the same space as the mapping . The multipole coefficients are evaluated at a reference radius of .

We focus on the uncertainty in the quadrupole gradient, defined as where is the -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 -smoothness, are made for the shape uncertainty and, therefore, we resort to a worst-case analysis. The maximum deviation of from its design value 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 , see Table 1

. There, the error on each level is estimated with respect to a fine discretization consisting of

total DOF, denoted as . Additionally, the numerical computation of the shape gradient is verified. To this end in Table 1 the maximum deviation with respect to a finite difference gradient computation, denoted , is given. Both estimated errors are found to be sufficiently small.

Parameter level level level level
Table 1: 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 is computed with DOF.

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 2 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 and by directly solving the optimization problem, denoted . Here, for the latter case MATLAB’s fmincon routine is used to carry out sequential quadratic programming. We observe a difference of about , 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.

num. free CP / / rel. diff.
Table 2: Quadrupole error by linearization and error estimated by means of MATLAB’s fmincon. indicates the worst-case scenario computed by Taylor expansion, the one obtained by solving the optimization problem. Perturbation magnitude of .
Figure 12: 3D model of one half of the Stern-Gerlach magnet with Rabi-type pole tips (modeled with CST EM STUDIO [41]), see [61].

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, 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


where is the non-linear reluctivity and 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 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.

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

For the optimization, the quantities of interest in the beam area of the magnet are the average magnetic field gradient and the inhomogeneity which are given by




where is the magnetic field gradient. The goal function used for the optimization process is


where is a free parameter, and , 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 controls the preference of one over the other quantity. The overall optimization problem can be written as


where 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 . This way, the magnetic flux density and even the average magnetic field gradient 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 3. 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].

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

4.4 IGA and the Boundary Element Method

The framework of Isogeometric Analysis is not limited to FEM. In recent years, Isogeometric Boundary Element Methods (BEM) gained a lot of attraction in research. The idea of bem is the representation of the solution 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


is given by


on where denotes some unknown density on the boundary. This density can be found by solving the arising variational problem: find in such that

where denotes an appropriate function space on the boundary and is the restriction of to the boundary. Choosing a suitable 2D spline space  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 such that

After inserting the canonical basis of spline space we receive the linear system


where encodes the discretization of the Dirichlet data given in (30), and encodes the desired density function. Plugging the discrete density back into (31), one can evaluate the solution of (30) in any given point , 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.

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.

CST (3D, original geometry) 0.0503
CST (3D, optimized geometry) 0.0201
GeoPDEs (2D, original geometry) 0.0477
GeoPDEs (2D, optimized geometry) 0.0122
Table 3: Average magnetic field gradient and inhomogeneity factor before and after the optimization.

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].

5 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 2- and 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.

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


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 ’Excellence Initiative’ of the German Federal and State Governments and the Graduate School of CE at TU Darmstadt.


  • [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.
  • [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.
  • [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.
  • [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 T-spline 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), 523-554, 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.
  • [28] A.F. Sarmiento, A.M.A. Cortes, D.A. Garcia, L. Dalcin, N. Collier and V.M. Calo, "PetIGA-MF: a multi-field high-performance toolbox for structure-preserving B-splines spaces", arXiv:1602.08727, 2016.
  • [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.
  • [30] N.D. Manh, A. Evgrafov, J. Gravesen and D. Lahaye, "Iso-geometric 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, D-50321 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 Finite-Element-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.
  • [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.
  • [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.
  • [53] S. Clénet, "Uncertainty Quantification in Computational Electromagnetics: The stochastic approach", ICS Newsletter (International Compumag Society), 13, 3-13, 2013.
  • [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.
  • [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 Laplace- and 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

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.
Carlo de Falco from MOX – Modellistica e Calcolo Scientifico, Dipartimento di Matematica, Politecnico di Milano, Milano, Italy.
Jürgen Dölz from the Departement Mathematik und Informatik, Universität Basel, Basel, Switzerland.
Rafael Vázquez from the Institute of Mathematics, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.
Ulrich Römer from the Department of Mechanical Engineering of Technische Universität Braunschweig, Braunschweig, Germany.
Corresponding author: Sebastian Schöps (