Isogeometric Methods for Computational Electromagnetics: B-spline and T-spline discretizations

Isogeometric Methods for Computational Electromagnetics:
B-spline and T-spline discretizations

A. Buffa G. Sangalli R. Vázquez Istituto di Matematica Applicata e Tecnologie Informatiche ’E. Magenes’ del CNR
via Ferrata 1, 27100, Pavia, Italy
Dipartimento di Matematica, Università di Pavia, via Ferrata 1, 27100, Pavia, Italy

In this paper we introduce methods for electromagnetic wave propagation, based on splines and on T-splines. We define spline spaces which form a De Rham complex and, following the isogeometric paradigm, we map them on domains which are (piecewise) spline or NURBS geometries. We analyse their geometric structure, as related to the connectivity of the underlying mesh, and we give a physical interpretation of the fields degrees-of-freedom through the concept of control fields. The theory is then extended to the case of meshes with T-junctions, leveraging on the recent theory of T-splines. The use of T-splines enhance our spline methods with local refinement capability and numerical tests show the efficiency and the accuracy of the techniques we propose.

Maxwell equations, De Rham diagram, Exact Sequences, Isogeometric Methods, Splines, T-splines.


1 Introduction

Electromagnetic field computations and, more generally, the numerical discretization of equations enjoying a relevant geometric structure, is one of the most interesting challenge of numerical analysis for PDEs and several results have been obtained in the last decade. Indeed, only for Galerkin methods, three Acta Numerica overview papers have been published: by Hipmair HIP02a (), by Arnold, Falk and Winther AFW06 (), and by Boffi Boff10 (), addressing different aspects of the problem.

On the one hand, discrete schemes have to preserve the geometric structure of the underlying PDEs in order to avoid spurious behaviors, instability or non-physical solutions (see e.g., the pioneering paper BFGP ()). For electromagnetics, as it is clear from the references above, numerical schemes have to be related with a discrete De Rham complex. On the other hand, especially in view of high frequency computations, numerical schemes need to be efficient and accurate. This requires many features, and among others it requires adaptivity, or at least local mesh refinement capability, in order to capture the strong singularities of the electromagnetic field, possibly driven by a-posteriori error estimator as, e.g., in BrSc08 ().

In this paper we present and analyse discretization techniques for electromagnetic fields based on splines and generalizations of splines, as NURBS (Piegl ()) or T-splines (Sederberg_Zheng () or below). Our work originates from IsoGeometric Analysis (IGA), IGA-book (). Isogeometric analysis has been introduced in 2005 by Hughes and co-authors in the seminal paper HCB05 () to solve structural mechanic problems directly on the geometry output by a CAD system, and has set the paradigm to use splines, NURBS or their generalization as generating functions for the construction of Galerkin spaces. This idea has been proved to be extremely effective and IGA is spreading very fast across different scientific communities: structural mechanics (see e.g., IGA-book (), MR2443159 (), Cottrell_Reali_Bazilevs_Hughes (), Echter2010374 (), Lipton_Evans_Bazilevs (), rank12 (), MR2835758 ()), geometric modeling (see e.g., DQS10 (), Martin_Cohen_Kirby (), CoMa10 () and also Zhang2012 ()) and numerical analysis (see e.g., BBCHS06 (), Beirao_Cho_Sangalli () , BRSV11 (), nwidth-iga (), Vuong_giannelli_juttler_simeon (), BePa12 (), kleissieti ()).

In this paper we present the recent advances in the use of the isogeometric paradigm and spline-based methods for electromagnetic wave computations. This research has started with the two papers Buffa_Sangalli_Vazquez () and BRSV11 () and can likely be considered as still in infancy (see also ratnani () for the applications of this results). This paper aims at showing the potential impact of these techniques in the electromagnetic community by addressing several aspects: from the geometric structure of the proposed methods, to local refinement strategies.

We introduce the spline complex studied in BRSV11 () (see (38) and (39)) and we present its properties: we construct canonical bases so that the matrices representing differential operators are the incidence matrices of the underlying meshes, and this enlightens the relation between the spline complex and the geometry of the underlying meshes. We show that for different choices of the degree of splines, the spline complex is isomorphic to the co-chain complex or to the chain complex of the underlying mesh. Besides this interesting fact, we also introduce the concept of control fields in analogy to control points which are ubiquitous in spline theory (see e.g., cohen2001geometric () or DeBoor ()) which provide the correct physical interpretation of degrees of freedom. Finally, we extend the results of BRSV11 () to multi-patch geometries, i.e., geometries which are piecewise spline or NURBS mapping of the unit cube. We refer the reader to kleissieti () for a detailed description of this class of geometries.

The second major contribution in this paper is a step towards adaptivity for spline-based methods. Leveraging on the recent work on T-splines, we design a two dimensional T-spline complex where meshes with T-junctions can be used to allow for adaptivity. T-splines are the most attractive way to break the tensor product structure of splines while keeping their structure and their accuracy. T-splines have been introduced in Sederberg_Zheng () and Sederberg_Cardon () and their use as a fundamental tool to enhance isogeometric analysis with adaptivity has been proposed in Bazilervs_Calo_Cottrell_Evans (). A series of papers has followed Buffa_Cho_Sangalli (), Scott_Li_Sederberg_Hughes (), Wang2011477 (), Beirao_Buffa_Cho_Sangalli (), together with the relevant class of Analysis Suitable T-splines LZSHS12 (), BBCS12 () which we use in our construction. The two dimensional T-spline complex is used to treat three dimensional problems with symmetry. We should also mention that the definition and use of T-splines in three dimension are not yet well understood, but object of an intensive study. Their use will allow, on a longer time perspective, to design full adaptive algorithms, on very general geometries parametrized on totally unstructured meshes. We refer the reader to Scott () for a monograph on the modern use of T-splines in geometric modeling.

Finally, we should remark that the spline spaces we study in this paper have a wide domain of applications and can be applied successfully to the discretization of other problems than electromagnetics. In fact, they can be used to solve the Darcy flows equations or more generally the Hodge laplacian operator as detailed in AFW06 () and AFW-2 (). Moreover, thanks to the regularity of spline spaces, their use in fluids is very promising. In the paper Buffa_deFalco_Sangalli () they are used for the first time to solve the Stokes equations, in EvHu12-4 () the Stokes eigenvalue problem is addressed, and in the sequence of three papers EvHu12 (), EvHu12-2 (), EvHu12-3 () they are applied to solve Stokes and Brinkman equations, steady and unsteady Navier-Stokes equations, providing impressive results.

The outline of the paper is the following. In Section 2 we set up the notation for the problems we address, in Section 3 we present known results about splines and NURBS in a self-contained way; in Section 4 we present the spline complex and all the related results while in Section 5 we introduce the T-spline complex and analyse its properties. Finally, in Section 6 we present numerical results: the first ones are two and three dimensional, academic tests aiming at demonstrating the validity of the proposed approach. As a last example, we compute the propagation in a waveguide with geometric inhomogeneity, on a three dimensional locally refined mesh.

2 Notation

In this section we present the notation that we need to describe the time-harmonic Maxwell problem. Let be a bounded Lipschitz domain. We denote by the space of complex square integrable functions on , endowed with standard norm , and by their vectorial counterparts. The Hilbert space contains functions of such that their first order derivatives also belong to . We denote by the subspace of functions with homogeneous boundary condition. We will also make use of the space , constituted by all functions in such that their curl also belongs to , and , the space of functions in such that their divergence belongs to . Moreover, we denote by (resp. ) the subspace of (resp. ) of fields with vanishing tangential (resp. normal) component.

For the sake of simplicity, we assume that the domain , referred to as physical domain in the following, is bounded Lipschitz and simply connected, and that its boundary is connected. We also assume that it is defined through a continuously differentiable parametrization with continuously differentiable inverse which we denote as , where will be referred to as the parametric domain. Further assumptions on the geometrical mapping will be given later.

Some notation will be borrowed from the context of differential forms: first of all, we define the spaces

Since the parametrization and its inverse are smooth, we can define the pullbacks that relate these spaces as (see (HIP02a, , Sect. 2.2)):


where is the Jacobian matrix of the mapping . Then, due to the curl and divergence conserving properties of and , respectively (see (Monk, , Sect. 3.9), for instance), the following commuting De Rham diagram is satisfied (see (HIP02a, , Sect. 2.2)):


We are also interested in spaces with boundary conditions, denoted with the subindex ,

for which the De Rham diagram reads


which also expresses the integral preserving property of .

Remark 2.1

As it is well known, the exactness of the sequences (2) and (3) relies on the assumption that (and ) has a trivial topology. All what we develop in this paper applies in principle also to the case of arbitrary topology but we do not present all the details here.

3 Preliminaries on splines and NURBS

We give here a brief overview on B-splines and, in the spirit of Bazilervs_Calo_Cottrell_Evans (), we also introduce some concepts that will be needed in the definition of T-splines. For more details on B-splines we refer the reader to HCB05 (); Piegl ().

3.1 Univariate B-splines

3.1.1 Knot vector and B-spline functions, refinement, spline derivatives

Given two positive integers and , we say that is a -open knot vector if

where repeated knots are allowed and denote by the multiplicity of the knot . We assume for all internal knots.

From the knot vector , B-spline functions of degree are defined following the well-known Cox-DeBoor recursive formula: we start with piecewise constants ():


and for the B-spline functions are defined by the recursion


This gives a set of B-splines that form a basis of the space of splines, that is, piecewise polynomials of degree with continuous derivatives at the internal knots , for . We denote this univariate spline space by


An example of some B-splines is given in Figure 1.

Figure 1: Example of B-splines of degree (left) and (right).

Notice that the B-spline function is supported in the interval , and in fact its definition only depends on the knots within that interval. For this reason, we define the local knot vector , and we will sometimes denote .

In the context of splines, three kinds of refinement are possible, as explained in HCB05 ():

  1. -refinement which corresponds to successive application of the Cox-DeBoor formula (4)–(5). Regularity is raised together with the degree: therefore, the spaces (6) are not nested under -refinement but, at each step (degree and regularity elevation), the dimension of the space increases by . The name -refinement has been introduced in HCB05 ();

  2. -refinement which corresponds to mesh refinement and is obtained by knot insertion. Let be the knot vector after inserting the knot in . Then, the new B-spline functions are constructed as follows:


    where , if , , if and, for . When is equal to or or to both, the knot insertion corresponds to reduction of the inter-element regularity at .

  3. -refinement which corresponds to the degree raising with fixed interelement regularity, and generates a sequence of nested spaces.

Assuming the maximum multiplicity of the internal knots is less than or equal to , i.e., the B-spline functions are at least continuous, the derivative of the B-spline is a spline as well. Indeed, it belongs to the spline space , where is a ()-open knot vector. Obviously, the regularity of splines in is one less than the regularity in .

In the following we assume that and . The domain of definition of the spline functions is the one-dimensional parametric domain. On it, the knot vector induces a partition of the interval that we denote by . Precisely, we define as the set of the knot spans , , that can also be empty due to knot multiplicity greater than . Empty intervals still play a role in the definition of B-splines and are graphically represented as points close one to the other, as proposed in Buffa_Cho_Sangalli (). Note that in this representation of , the number of lines is the knot multiplicity with one exception: for each boundary knot (at 0 or 1) of an open knot vector in we represent only a multiplicity of , which is lines if is odd, and lines if is even (see Figure 1). The reason for this construction of will be motivated in the next section.

Finally, it is worth noting the relationship between the space and the space of derivatives , and their respective meshes and . The meshes and may differ only as regards the number of points at the boundary. Indeed, according to the definition above, if is odd both meshes coincide, and if is even the number of elements of with respect to is reduced by two, one on each side.

3.1.2 Anchors and Greville sites

In this section we present the concept of anchors and of Greville sites as points in the parametric space which may be associated to each B-spline function. Greville sites, which are also known as knot averages, are classical and can be found for instance in DeBoor (), while the concept of anchors has been introduced recently in Bazilervs_Calo_Cottrell_Evans ().

Since splines are not interpolatory, the association of functions to points (or, as we will see, other geometric entities) is somehow more arbitrary than with Lagrangian finite elements. Anchors and Greville sites are two different choices, and we present here both. Anchors are very useful when dealing with non-tensor product extensions of splines as T-meshes, while Greville sites (and related geometric entities) carry

degrees of freedom in a more natural way.

Given a B-spline function , and its local knot vector , we set: if is odd, the anchor associated with is the central knot of . If is even, the anchor associated with is chosen to be the midpoint of the central knot span of , namely: . The position of the anchors for degrees and are represented in Figure 1.

Note that obviously the correspondence between anchors and B-splines functions is one to one, but different anchors may lie at the same position . A remedy to this abuse of notation, at the cost of more complex definition, is proposed in BBSV12 () where the use of both an index and a parametric domain is proposed.

The set of anchors is denoted as . When is odd anchors are located at all knots of the partition (which may be repeated), while when is even anchors are located at midpoints of all elements in (including the ones of zero area). Indeed, this fact is the reason for the definition of in particular as regards to the multiplicity of boundary knots.

Most often, we will use anchors to index functions and local knot vectors. Namely, for an anchor , and will denote the corresponding local knot vector and B-spline function, respectively. When no confusion occurs, the subindex may be removed.

Remark 3.1

The B-splines are, in general, not interpolatory at the anchor , while they are interpolatory at knots having multiplicity . This always happens at and , and happens in the interior of the parametric domain where the basis is continuous, i.e., at knots with multiplicity . See e.g., Figure 1(left).

Given , and for some , then the Greville site is defined as:


Unlike anchor positions, Greville sites are all different one from the other, when the multiplicities verify and thus B-splines are all continuous. The Greville sites induce a partition of the unit interval, referred as Greville mesh and denoted . These concepts are ubiquitous in spline theory and geometry representation. Greville sites are intimately related to control points and control polygon whose properties we briefly recall in the next section.

3.1.3 B-spline curves

A B-spline curve in is defined by a parametrization in the interval , in the form


where are called the control points. Control points are in a one-to-one correspondence with B-spline basis functions. The piecewise linear interpolation of the control points gives the control polygon . See Figure 2 for an example.

The control points have an important role not only in the definition of the spline parametrization (15), but also in the visualization and interaction with spline geometries within CAD softwares. Indeed, it is common in CAD softwares to represent, together with the parametrized curve , the control points and the associated control polygon . Typically, the CAD user defines or interacts with the control points in order to input and modify the geometry. Since the B-splines are not in general interpolatory (recall Remark 3.1), then the control polygon differs from , but it is “close” to it. Precisely, converges to under -refinement. This convergence is proved, e.g., in DeBoor () and discussed here below.

We introduce the usual Lagrangian basis for piecewise linear polynomials on the Greville mesh , denoted by :


The control polygon is then parametrized by the mapping defined by


that is, and share the same control points. When is smooth enough, the following approximation estimate holds (see, e.g., (DeBoor, , Ch. XI)):


denoting the mesh-size. In other words, approximates up to an error under -refinement. A graphical representation of this convergence can be seen in Figure 2.

Figure 2: B-spline curve and its control polygon (left), and the same curve after one step of -refinement (right).

3.2 Multivariate B-splines

Multivariate B-splines are defined from univariate B-splines by tensor product, see for instance Piegl (); DeBoor (). Anchors are defined in a similar way. We give here a quick overview.

3.2.1 Knot vectors, B-spline functions, anchors, Greville sites

Let be the space dimensions (in practical cases, ). Assume , the degree and the -open knot vector are given, for . These knot vectors define a tensor product mesh in the parametric domain where, as in Section 3.1.2, we have to take into account the knot multiplicity. The multiplicity of a knot vector in is represented graphically by lines () or planes () one close to the other, while the boundary is treated exactly as in one dimension.

The set of anchors is defined on as the Cartesian product . Considering, for example, the trivariate case () and recalling the definitions from Section 3.1.2 for the univariate case, we have that: if all are odd the anchors lie at the vertices of the mesh; if both and are odd and is even, then the anchors are middle-points of the vertical edges of ; if both and are even and is odd, then the anchors are centers of the horizontal faces of ; if all are even the anchors lie at the center of the elements of , and so on. As in the univariate case, the anchors may be located at the center of zero length edges or zero area faces or empty elements, according to knot repetition. Also, the computation of the local knot vectors for each anchor follows from the univariate case. Given an anchor , we have that its coordinates are . The three local knot vectors (one in each coordinate direction) corresponding to are defined as for and the B-spline associated to is constructed by tensor product as:


with .

The B-spline functions (13) span the space (or simply ), which is the space of piecewise polynomials of degree in the direction on , whose continuity at the internal mesh plane is , being the multiplicity of in the knot vector .

To each anchor (or, equivalently, to each B-spline function (13)) we also associate a Greville site in the natural way, that is


where each is defined as in (8), from the local knot vector . Connecting adjacent Greville sites, we obtain the Greville mesh , which is a regular tensor product mesh with all elements of positive volume.

3.2.2 Spline and NURBS geometries, multi-patch domains

Analogously to spline curves, a trivariate single-patch spline parametrization of the domain is defined as a linear combination of B-splines,


where are called control points. In a similar way, it is also possible to define bivariate spline domains in or surfaces in , which are commonly used in CAD (see, e.g., Piegl (); cohen2001geometric ()).

The control points have again the same important role in the visualization and interaction with geometries within CAD softwares. Now, the concept of control polygon generalizes to the control mesh , that is, the mesh connecting the control points. Figure 3 shows an example geometry, with its control points and control mesh. The control mesh defines a polyhedral domain, denoted , which is an approximation of : again, the control domain converges to under -refinement.

Figure 3: Representation of a geometry (green), with its control points (blue) and control mesh (black lines) for splines of degree .

This is stated as in the univariate case: we introduce the usual Lagrangian basis for piecewise trilinear polynomials on the tridimensional Greville mesh , still denoted by , for the sake of brevity,

The control mesh is the image of the Greville mesh through the piecewise trilinear mapping ,


which is a parametrization of , since

When is smooth enough, as for (12), we have


The control mesh plays a fundamental role in structural mechanics applications where the unknowns are sought as displacements of control points. In our work, we will show how this interpretation can be used also in other contexts.

Remark 3.2

When (and all anchors have multiplicity one) the Greville sites coincide with the anchor representations, i.e., , and , , that is, and coincide.

In CAD and isogeometric analysis the geometry is often parametrized by Non Uniform Rational B-splines (NURBS). NURBS are generated from projective transformations of splines (see Piegl ()). A trivariate single-patch NURBS parametrization of the domain is a function defined as quotient of linear combination of B-splines,


where are the NURBS control points and the positive NURBS weights.

In order to enhance flexibility and allow for more complex geometries, the definition of tensor-product spline and NURBS parametrized domain can be generalized to domains that are union of images of cubes; precisely


where the are referred to as patches, and are assumed to be disjoint. Each patch has its own parametrization , defined on its own spline or NURBS space. The whole is then referred to as a multi-patch domain. For the construction of discrete fields on a multi-patch domain we will introduce in Section 4.4 suitable conformity assumptions. These will restrict the framework to configurations where it is easy to implement the proper continuity of the fields at the patches interface.

In this paper, is assumed to be parametrized either by spline or NURBS functions but the unknown fields are always constructed by splines. This means that, in case of NURBS geometries, we leave the isoparametric concept which is a fundamental assumption for isogeometric methods in the context of continuum mechanics (see IGA-book ()).

4 The spline complex

This section is devoted to present the spline spaces that are compatible with the De Rham complex. The definition of the spaces is taken from BRSV11 (), and is given in three dimensions (though the same construction is generalizable to arbitrary dimension). We first recall the construction on the parametric domain , and then the discrete spaces on the physical domain are obtained by the push-forward mapping associated to (1). As shown in BRSV11 (), it is also possible to complement these spaces with commuting and continuous projectors, in the setting of the so called Finite Element Exterior Calculus (see AFW06 ()), however this issue is not discussed here. Instead, we discuss the selection of a suitable basis for the implementation of the proposed spaces, and the meaning of the associated degrees-of-freedom. We will see that the proposed spline spaces are a natural high-order extension of classical low-order Nédélec hexahedral finite elements of the first family (see NED80 ()), obtained in this setting for degree , and that in a natural way they are related to cochain or chain complexes of the mesh where they are defined.

4.1 Complex on the parametric domain

We recall the following property of univariate splines, from Section 3.1.1: the derivative of a (continuous) spline is a spline, and in particular


where is the -open knot vector that coincides with the -open knot vector except for the boundary knot repetitions. Moreover, we have that the derivative of the B-spline associated to an anchor in is a linear combination of the B-splines associated to the previous and next adjacent anchors and in (only one adjacent anchor for the first and last ); precisely, denoting by the local knot vector (formed by knots) of and by the local knot vectors (formed by knots) of , and with the general notation of Section 3.1.1, we have


where are the length of the support of the -degree B-splines , that is, the difference of the last and first knots in the local knot vectors . An example is given in Figure 4. When is the first (resp., last) anchor, (21) holds with (resp., ). This is a well known property of B-splines (see Piegl (); DeBoor ()) that also suggests the following scaling of the basis functions of and


The scaling in (23) gives the Curry-Schoenberg B-splines (see (DeBoor, , Ch. IX)), that have been already used in isogeometric analysis in ratnani (). Indeed, with the bases (22) and (23), the matrix associated to the operator is the edge-vertex incidence matrix related to the mesh , when is odd, or the vertex-edge incidence matrix related to , when is even. We recall that also contains zero length edges and repeated vertices.


Figure 4: Derivative of the spline associated to the anchor as a linear combination of the splines associated to and .

The observations above are the key ingredients of the trivariate construction. Following BRSV11 (), and using the notation of Section 3.2, we introduce the following discrete spaces on the parametric domain


From (20), they form a De Rham complex:


Moreover, we have the following result.

Theorem 4.1

The sequence (25) is exact.


This result has been already presented in BRSV11 (). We present an alternative proof, that will be useful Section 5.3.

We have to show that in (25) it holds


In particular, we have to prove the inclusion in (26)–(29), since the other inclusion is trivial in all cases. It is also trivial that

Let , then define


it is easy to check that when , and that ; then

Consider , and define such that


as before, it is easy to check that and that ; then

In order to complete the proof we need to show that

which is implied by


To count dimensions recall from Section 3.1.1 that , ; then from Section 3.2 and from (24) we get


Then, by (26)–(27),

and by (29)