1 Introduction↩︎

Subdivision surfaces have emerged as a promising alternative to standard CAD technologies in recent years. They overcome the topological limitations of multi-patch NURBS surfaces, while still providing a smooth surface representation. Besides being a powerful tool for geometric modeling, subdivision methods naturally extend to application in isogeometric analysis (IGA). First showcased in [1], the Loop subdivision scheme was utilized for discretization of the Kirchhoff-Love thin-shell formulation, satisfying the regularity requirements. While electromagnetic problems typically only require continuity, the integration of IGA with subdivision surfaces may still be advantageous for the simulation of various problems, in particular where smoothness of the numerical solution for the electromagnetic fields is a desirable property, e.g. particle tracking.

Most relevant electromagnetic formulations require not just the scalar space \(\Hgrad\), but also the vectorial spaces \(\Hcurl\) and \(\Hdiv\) (or variants of them) of a de Rham sequence. These spaces need special treatment when it comes to discretization, requiring vectorial edge and face basis functions, respectively [2]. The classical subdivision schemes (e.g. Catmull-Clark [3], Loop [4]) however generate scalar \(\Hgrad\) limit functions, which are not applicable for a consistent \(\Hcurl\) discretization.

There are specialized edge subdivision schemes in the literature, first introduced by [5] [5], trying to generalize subdivision to the discrete exterior calculus setting, by enforcing commutativity between subdivision and the exterior derivatives. The proposed edge subdivision schemes, based on the Loop and Catmull-Clark methods, are further investigated in [6], where application in the context of numerical analysis, in particular for IGA, is investigated. Recently, [7], [7] provide a rigorous mathematical framework for the generalization to \(k\)-form subdivision schemes, and even apply the proposed methods to solve a two-dimensional Maxwell eigenvalue problem. However, the works above do not use the subdivision limit functions for the IGA discretization and hence do not obtain increased inter-element regularity for the numerical solution. ³ To the authors’ best knowledge it is still unknown, how the limit functions from these subdivision schemes perform and how to evaluate them exactly.

A second challenge is the extension to three-dimensional problems. An isogeometric approach using the boundary element method (BEM) is generally well suited for the application of the subdivision limit functions, since discretization of the boundary data requires only bivariate basis functions, which subdivision surfaces naturally provide. This approach is for example investigated in [8], [9]. A volumetric discretization by the finite element method (FEM) however, requires trivariate basis functions, making subdivision surfaces inapplicable. There do exist extensions of well known subdivision algorithms to subdivision solids, and they have been applied in the context of IGA [10], but their analytical properties are even less understood than their two-dimensional counterparts [11]. Also, to our knowledge, there exist to this date no generalization of edge subdivision schemes, required for an \(\Hcurl\) discretization, in the trivariate case.

Considering the above mentioned challenges, this paper investigates the following approach: For the reduction from 3D to 2D, we exploit the symmetry of an axisymmetric problem setup. This is a reasonable assumption for various problems of practical relevance, such as simulation of accelerator components. For the Maxwell eigenvalue problem, the approach results in a problem formulation in \(\HgradAxi\), which can be discretized with standard subdivision limit functions, without the need for a specialized edge subdivision scheme.

The paper is structured as follows: First, in 2 the Maxwell eigenvalue problem is formulated and the simplification under axisymmetry is performed. In 3 the basics of Catmull-Clark subdivision, and especially the construction of the basis functions, are introduced. Discretization of the problem formulation using subdivision IGA is then performed in 4. All numerical results are presented in 5.

2 Axisymmetric Maxwell eigenvalue problem↩︎

Macroscopic electromagnetic phenomena are governed by Maxwell’s equations. Under certain physical assumptions, the equations can be simplified to boundary value problems solvable using numerical methods. The Maxwell eigenvalue problem we consider in this section, is an important formulation describing electromagnetic wave propagation in for example cavities and waveguides [12].

Consider a sufficiently smooth domain \(\Omega\subset\real^3\), constant materials \(\varepsilon, \mu > 0\) and zero sources. The three-dimensional time-harmonic Faraday’s and Ampère’s laws are then given by the two equations \[\begin{align} \curl\vb{E} &= -i \omega\mu\vb{H} \tag{1} \\ \curl\vb{H} &= i \omega\varepsilon\vb{E} \tag{2}, \end{align}\] with the eigenfrequency \(\omega\) and imaginary unit \(i\) [12]. Inserting 2 into 1 yields the \(\vb{H}\)-based Helmholtz wave equation \[\curl\curl\vb{H} = \omega^2\mu\varepsilon \vb{H} . \label{eq:helmholtz}\tag{3}\] To numerically compute solutions of 3 , a variational formulation is typically used, and for well-posedness of the problem, appropriate boundary conditions must be imposed. We consider homogeneous Dirichlet and homogeneous Neumann boundary conditions here. For notational simplicity, the boundary conditions are not explicitly expressed in the function space notation. The weak \(\vb{H}\)-formulation of the Maxwell eigenvalue problem for \((\lambda, \vb{H}) \in \real_{>0}\times\Hcurl\) is then given by \[\forall \vb{H}'\in\Hcurl: (\curl \vb{H}, \curl \vb{H}') = \lambda (\vb{H}, \vb{H}'), \label{eq:maxwell95ev95weak}\tag{4}\] where \(\vb{H}'\) is a test function [2]. We denote by \(\lambda = \omega^2 \varepsilon\mu\) the eigenvalue and by \((\cdot,\cdot)\) the inner product on \(L^2(\Omega)\) and \((L^2(\Omega))^3\).

The following analysis is restricted to the axisymmetric eigenvalue problem. In particular we assume that the domain \(\Omega\) is a solid of revolution of an axisymmetric slice \[\require{physics} S := \qty{ (\rho,z) \in\real_{\geq 0}\times\real : (\rho,0,z)\in\Omega } .\] Under this assumption the derivatives \(\partial_\varphi\) vanish and formulation 4 simplifies significantly. The above function spaces are then to be replaced by two-dimensional suitably weighted versions of them (indicated by a \(\rho\)-subscript). Considering magnetic fields with only a \(\varphi\)-component, i.e. \(\vb H = u \vb{e}_\varphi\), the formulation can be reduced to \[\forall v \in\HgradAxi[S]: (\rot u, \rot v)_\rho = \lambda (u, v)_\rho, \label{eq:maxwell95ev95weak95axi}\tag{5}\] where \(\rot: \HgradAxi[S]\to\HdivAxi[S]\) is the rotated gradient and \[\require{physics} (u,v)_\rho = \iint_S u v \,\rho\dd{\rho}\dd{z}\] is the two-dimensional \(\rho\)-weighted inner product on \(\LtwoAxi[S]\) and \(\require{physics} \qty(\LtwoAxi[S])^2\).

For any given eigenpair \((\lambda,u)\), the electric field can be computed by Ampère’s law as \[\HrotAxi[S]\ni\vb{E} = \frac{-i}{\omega\varepsilon} \rot u,\] with \(\rotScal: \HrotAxi[S]\to\LtwoAxi[S]\) being the scalar curl.

Note that this ansatz can be generalized by expanding the solution in terms of Fourier modes, see e.g. [13]. For simplicity, this is not investigated here.

3 Catmull-Clark subdivision↩︎

There exists a variety of different subdivision algorithms, all with different geometric and numerical properties. For a general introduction to subdivision, we refer to [14]. One of the most well known and conceptionally easy to understand subdivision schemes is Catmull-Clark subdivision. A Catmull-Clark subdivision surface is an almost-everywhere \(C^2\) surface, defined on a possibly unstructured quadrilateral mesh. Restricted to tensor-product patches, the surface coincides with classic cubic B-spline surfaces [14].

The Catmull-Clark subdivision algorithm is originally formulated in [3] and is understood as a refinement algorithm of a coarse control mesh. The iterative application of the subdivision stencil to the initial mesh will then, in the limit, converge to a smooth surface. Similar to the well-known B-spline surfaces, this limit surface can also be represented by a linear combination of basis functions with the initial control points. This point of view is crucial for the application of subdivision to isogeometric analysis.

To define the Catmull-Clark basis functions, first consider the parametric domain \[\require{physics} \widehat{S} := [0,1]^2 \times \qty{1,\dots,N}\] of \(N\) cells. All cells are equipped with adjacency information, by topologically identifying the corners and edges of neighboring cells. The mesh defined by these neighborhood relations is allowed to be unstructured. In that case, the valence (number of outgoing edges) of a mesh vertex is in the following denoted by \(n\) and equals \(4\) for regular interior vertices. Vertices of any other valence are called extraordinary.

Similar to classical B-spline theory, the subdivision basis functions on each cell \(\Sigma\in\widehat{S}\) have support extending onto the neighboring cells. While these basis functions are defined by the limit of the subdivision process, they can also be evaluated parametrically, as shown in [15]. In the following, for the construction of a parametric representation of the basis functions, we need to differentiate between regular and irregular cells, i.e. cells without or with an extraordinary vertex, respectively. Cells at the boundary of the mesh also need special treatment and are usually handled by extrapolation of control points. We do not discuss this case here; details can be found in [16].

3.1 Local regular basis functions↩︎

For a cell without an extraordinary vertex, consider the four uniform cubic B-splines \(N_i: [0,1]\to\real\), defined by \[\begin{align} N_1(\xi) &= \frac{1}{6} \left(1 - 3\xi + 3\xi^2 - \xi^3 \right) \\ N_2(\xi) &= \frac{1}{6} \left(4 - 6\xi^2 + 3\xi^3 \right) \\ N_3(\xi) &= \frac{1}{6} \left(1 + 3\xi + 3\xi^2 - 3\xi^3 \right) \\ N_4(\xi) &= \frac{1}{6} \xi^3. \end{align}\] The Catmull-Clark basis functions on this regular cell then coincide with the \(4\times4\) tensor-product B-splines \[\widehat{b}_{ij}: [0,1]^2 \to\real \quad (\xi, \eta) \mapsto N_i(\xi) N_j(\eta) .\] On each regular cell \(\Sigma\), these local basis functions are typically sorted lexicographically and denoted with a single linear index as \(\widehat{b}_{\Sigma,i}\). Each of these basis functions is associated with a vertex in the \(3\times3\) grid of neighboring cells. The 16 vertices of this grid and the basis functions along both parametric directions are schematically depicted in 1.

3.2 Local irregular basis functions↩︎

Handling irregular cells is more involved than their regular counterpart, and lies at the heart of each subdivision scheme. Because for such an irregular region, no structured \(3\times3\) grid of cells can be found, the tensor-product B-splines cannot meaningfully be defined. Instead, the idea is to subdivide repeatedly, until the parameters \((\xi, \eta)\) lie inside a regular cell and then evaluate the regular bi-cubic B-splines. Mathematically, this is achieved as follows: Consider an infinite family of tiles \(\Sigma_k^\ell \subset [0,1]^2\) partitioning the unit square, as depicted in 2. The upper index \(\ell > 0\) denotes the level of subdivision and the lower index \(\require{physics} k \in\qty{1,2,3}\) denotes the position inside the refined square (bottom right, top right, top left). For \((\xi, \eta) \in\Sigma_k^\ell\) consider the transformations \(t_k^\ell: \Sigma_k^\ell \to [0,1]^2\), defined by \[\begin{align} t_1^\ell(\xi, \eta) &= (2^\ell \xi - 1, 2^\ell \eta) \\ t_2^\ell(\xi, \eta) &= (2^\ell \xi - 1, 2^\ell \eta - 1) \\ t_3^\ell(\xi, \eta) &= (2^\ell \xi, 2^\ell \eta - 1). \end{align}\] These transformations simply map the parametric coordinates \((\xi, \eta)\) from a refined tile to the standard unit square in an affine way. This mapping is required to evaluate the B-splines on each tile.

Next, consider the extended subdivision matrices \(\vb A\) and \(\bar{\vb A}\), as defined in [15]. The matrices handle the subdivision process, by mapping initial control points to control points of the refined mesh. This is illustrated in 3, where the vertices of the initial mesh are shown in black, and the refined ones in gray. The dark gray vertices correspond to the matrix \(\vb A\) and the light gray ones to \(\bar{\vb A}\).

The \(2n+8\) local basis functions \(\widehat{b}_{\Sigma,i}: [0,1]^2 \to\real\) for an irregular cell \(\Sigma\), collected in a row-vector \(\widehat{\vb b}_{\Sigma} = (\widehat{b}_{\Sigma,i})_i\), are then defined by restriction onto each tile \(\Sigma_k^\ell\) as \[\require{physics} \widehat{\vb b}_{\Sigma_k^\ell} := \qty(\widehat{\vb b}_k \circ t_k^\ell) \overline{\vb A} \vb A^{k-1}. \label{eq:catmark95basis95irregular}\tag{6}\] For each \(k \in\{1,2,3\}\), the row-vector \(\widehat{\vb b}_k: [0,1]^2 \to\real^{2n+17}\) in 6 is a collection of the bi-cubic B-splines, sorted to match the ordering of the vertices in the subdivision process (the light gray ones in 3). As for the regular case, each irregular basis function \(\widehat{b}_{\Sigma,i}\) is associated with one of the \(2n+8\) vertices in the neighborhood of the irregular cell. For a more detailed explanation of the definition and evaluation of the irregular basis functions, see [15].

3.3 Surface parametrization↩︎

Combining the constructions for the regular and irregular basis functions, a parametrization for the subdivision surface can be given. Associating a control point \(\vb{c}_i \in\real^3\) with each of the \(2n+8\) vertices of a cell \(\Sigma\), each local patch of the subdivision surface \(S\) can be parametrized by \[\Phi: \Sigma\to S \quad (\xi, \eta) \mapsto \sum_{i=1}^{2n+8} \widehat{b}_{\Sigma,i}(\xi, \eta) \vb{c}_i .\] Note that for the construction of the global parametrization \(\Phi:\widehat{S}\to S\), control points of vertices identified by the neighborhood relation are required to match.

4 Subdivision-based isogeometric discretization↩︎

In this chapter, we introduce the discretization of the problem formulation by isogeometric analysis, i.e. by utilizing the Catmull-Clark basis functions for the discretization of both the domain and the solution fields.

4.1 Discrete function spaces↩︎

The Catmull-Clark limit functions, defined on each local cell in 3, are suitable basis functions for an isogeometric Galerkin discretization of the weak forms in 5 . On each cell, let \[b_{\Sigma,i}: \Phi(\Sigma)\to\real \quad \vb{p} \mapsto \widehat{b}_{\Sigma,i} (\Phi^{-1}(\vb{p}))\] be the local basis functions on a physical element, assuming that the pullback \(\Phi^{-1}\) exists. For some arbitrary but fixed global ordering of the vertices in \(\widehat{S}\) (where vertices identified by the adjacency relation have the same index), the global subdivision basis functions are defined by association with a vertex each and restriction onto the elements. Denoting these functions by \(b_i: S\to\real\), the subdivision function space is defined as \[\require{physics} \mathcal{S}_h := \Span\qty{b_i},\] with \(h\) denoting the maximum element size. For mathematical details on global regularity and linear independence, see [18]–[20].

Construction of the subdivision function space in the physical domain enables the discretization of 5 . Consider the discrete representation of the solution \(u_h = \sum_i u_i b_i \in\mathcal{S}_h\). The formulation is then given by \[\forall i,j: \sum_i u_i \underbrace{(\rot b_i, \rot b_j)_\rho}_{=: k_{ij}} = \lambda \sum_i u_i \underbrace{(b_i, b_j)_\rho}_{=: m_{ij}}, \label{eq:maxwell95ev95iga95axi}\tag{7}\] which yields the generalized eigenvalue problem \[\vb K \vb u = \lambda \vb M \vb{u},\] with the stiffness matrix \(\vb K = (k_{ij})_{ij}\), mass matrix \(\vb M = (m_{ij})_{ij}\) and solution vector \(\vb u = (u_i)_i\).

4.2 Numerical integration↩︎

For the computation of the matrix elements \(m_{ij}\) and \(k_{ij}\), numerical quadrature is employed to approximate the value of the integrals with sufficient accuracy. The standard for FE or B-spline and NURBS based IGA is the use of Gaussian quadrature on every element [2], [21]. Since the basis functions for these methods are smooth in the interior of each element, Gaussian quadrature has optimal approximation power. This is however not the case for the subdivision basis functions: Due to the partitioning of the surface into infinitely many tiles near an extraordinary vertex, the irregular limit functions are only \(C^2\) inside the element. In order to compute the integrals accurately, a hierarchical quadrature scheme is employed. The basic idea is, to perform Gaussian quadrature not on the entire element at once, but on a finite number of refined tiles, up until some fixed level \(L \geq \ell\). For more details, see [22].

5 Numerical results↩︎

The performance of the method proposed in this paper is evaluated in this section, by numerically solving various model problems.

5.1 Convergence analysis↩︎

To numerically verify the discretization scheme, a convergence study is performed. The first problem to investigate convergence is the cylindrical pillbox cavity, which has closed-form solutions for the eigenmodes and eigenfrequencies [12]. A pillbox cavity of radius \(R\) and height \(d\) is depicted in 4. The eigenfrequencies of the transverse magnetic (TM) modes in the cavity are given by \[\omega_{mnp} = \frac{1}{\sqrt{\mu\varepsilon}} \sqrt{ \frac{x_{mn}^2}{R^2} + \frac{p^2 \pi^2}{d^2} }\] with \(x_{mn}\) being the \(n\)-th root of the Bessel function \(J_m\).

In the following, convergence for the computation of the TM\(_{010}\) mode is analyzed, for the parameters \(R = \SI{35}{\milli\meter}\) and \(d = \SI{100}{\milli\meter}\). The eigenfrequency of this mode equals \(f \approx\SI{3.278}{\giga\hertz}\) and the fields are given by \[\require{physics} \begin{align} H_\varphi(\rho,z) &= -i \sqrt{\frac{\varepsilon}{\mu}} E_0 J_1\qty(\frac{x_{01} \rho}{R}) \\ E_z(\rho,z) &= E_0 J_0\qty(\frac{x_{01} \rho}{R}) \end{align}\] with all other components equaling zero. The (non-unique) value \(E_0\) denotes the field amplitude. For the geometry discretization of the axisymmetric slice \([0,R]\times[0,d]\), a purposefully unstructured mesh is utilized, to highlight the approximation properties under the existence of extraordinary vertices. Gaussian quadrature of degree \(p=4\) is used for the numerical integration on regular elements. On irregular elements \(L=10\) levels of subdivision are used for the hierarchical quadrature scheme.

The results of the convergence study using Catmull-Clark elements are depicted in 5. The error in the computed eigenfrequency, and the \(\LtwoAxi[S]\) and \(\HgradAxi[S]\) errors of the eigenmode are plotted against the square root \(\sqrt{n_\mathrm{dof}}\) of the number of degrees of freedom (DOFs). As indicated by the slope triangles, the observed rates are approximately 1.5 in \(\HgradAxi[S]\) and 2.5 in \(\LtwoAxi[S]\). These rates agree with observations in the literature for a Catmull-Clark discretization of the \(H^1\)-Laplace problem, see e.g. [11], [17]. The study experimentally confirms that the axisymmetric problem exhibits the same convergence rates as the standard 2D problem and that Catmull-Clark basis functions are hence suitable for axisymmetric discretizations. The observed rate of 3.0 for the eigenfrequency error is also consistent with the literature, as we would expect a doubling of the \(\HgradAxi[S]\)-rate [24].

As a more complex test case, we consider the eigenmodes of a single-cell TESLA cavity. For the experiment, the geometry parameters from [25] for a midcup are used. The shape of the cavity is depicted in 6.

In contrast to the pillbox test case, the geometry of the TESLA cavity is not exactly representable by a Catmull-Clark subdivision surface, as the 2D cross section contains circular and ellipsoidal arcs. Therefore, after each refinement step, the boundary of the subdivision surface is fitted to the exact geometry. For simplicity the fitting is performed by cubic spline quasi-interpolation \[\vb{c}_i = \frac{1}{6}\left(-\vb{p}_{i-1} + 8\vb{p}_i - \vb{p}_{i+1} \right) ,\] where \(\vb{c}_i\) are the control points of the Catmull-Clark boundary curve and \(\vb{p}_i\) samples of the exact boundary (see [26]). More sophisticated methods, like \(L^2\)-projection or least-squares, are also conceivable, but are not investigated in this paper.

The convergence of our Catmull-Clark discretization is analyzed for the five lowest TM\(_{0np}\) modes. The utilized quadrature scheme is the same as for the pillbox test case and the mesh also contains extraordinary vertices. As there are no closed-form solutions for the eigenfrequencies, the values used in the error computation are obtained numerically using the open-source IGA framework GeoPDEs [27] with very high accuracy (using cubic B-splines with \(\sim\num{100000}\) DOFs).

The results are shown in 7. For all analyzed eigenmodes, the trend for the frequency error decay is roughly the same, differing only in a constant. Similar to the pillbox test case, the observed rates are all approximately 3, validating the applicability of the method for complex geometries.

5.2 Longitudinal field evaluation↩︎

As an illustration of the increased smoothness of the eigenmodes, we showcase the graph of the electric field inside a 9-cell TESLA cavity with respect to the longitudinal \(z\)-direction. The geometry parameters for the cavity cells are again taken from [25]. The radial component \(E_\rho\) of the field is evaluated along the horizontal lines \(\rho = R_\text{iris}/2\) and \(\rho = R_\text{eq}/2\) and plotted in 8 and 9, respectively. The plots showcase a comparison between lowest order FE and Catmull-Clark elements. Both discretizations use the same mesh with extraordinary vertices.

The results visually confirm that the Catmull-Clark discretization method produces accurate results for the electric field, as the solution qualitatively agrees with the standard finite element one. However, it can also be seen, that the solutions noticeably differ towards the peaks of the field: For lowest order \(C^0\) elements, the electric field becomes increasingly noisy, while the Catmull-Clark solution does not show overshoots. These unphysical jumps at element interfaces can be attributed to the discretization method, since \(C^0\) basis functions cannot guarantee a continuous solution for the derivative (here, the electric field). Because Catmull-Clark basis functions are \(C^2\) across elements, they guarantee continuity even of the derivatives. Using higher order \(C^0\) basis functions (e.g. Lagrange elements) is also expected to significantly suppress the observed numerical noise, resulting in a better field quality. However, this approach cannot exactly suppress the jumps in the derivative and requires more degrees of freedom on the same mesh.

6 Conclusion↩︎

In this work, we have presented a subdivision-based isogeometric discretization of the axisymmetric Maxwell eigenvalue problem. The discretization is based on the Catmull-Clark subdivision scheme, utilizing the subdivision basis functions for both the representation of the geometry and for the discretization of the axisymmetric unknowns. For TM modes, the resulting magnetic field is \(C^2\)-continuous except at extraordinary vertices, and the electric field is \(C^1\)-continuous.

Performance of the method has been demonstrated by various numerical experiments. A convergence analysis for the computation of TM modes in a pillbox cavity showed that the rates in \(L^2\) and \(H^1\) agree with observed rates for the two-dimensional Poisson problem in the literature (2.5 and 1.5 respectively). The eigenvalues of the pillbox and single-cell TESLA cavity also converged with the expected (doubled) rate.

The quality in terms of smoothness of the numerical fields of a 9-cell TESLA cavity has been investigated and compared to classical lowest order finite elements. The results show that the subdivision-based discretization yields significantly smoother fields, overcoming the interface discontinuities of the finite element solution, while requiring the same amount of degrees of freedom.

References↩︎