Finite Element Time-Domain Body-of-Revolution Maxwell Solver based on Discrete Exterior Calculus
FFinite Element Time-Domain Body-of-Revolution MaxwellSolver based on Discrete Exterior Calculus
Dong-Yeop Na a , Ben-Hur V. Borges b , Fernando L. Teixeira a a ElectroScience Laboratory and department of Electrical and Computer Engineering, The Ohio StateUniversity, Columbus, OH 43212, USA b Electrical and Computer Engineering Department, University of S˜ao Paulo, S˜ao Carlos, SP13560-970, Brazil
Abstract
We present a finite-element time-domain (FETD) Maxwell solver for the analysis ofbody-of-revolution (BOR) geometries based on discrete exterior calculus (DEC) of dif-ferential forms and transformation optics (TO) concepts. We explore TO principles tomap the original 3-D BOR problem to a 2-D one in the meridian ρz -plane based on aCartesian coordinate system where the cylindrical metric is fully embedded into the con-stitutive properties of an effective inhomogeneous and anisotropic medium that fills thedomain. The proposed solver uses a (TE φ , TM φ ) field decomposition and an appropriateset of DEC-based basis functions on an irregular grid discretizing the meridian plane.A symplectic time discretization based on a leap-frog scheme is applied to obtain thefull-discrete marching-on-time algorithm. We validate the algorithm by comparing thenumerical results against analytical solutions for resonant fields in cylindrical cavitiesand against pseudo-analytical solutions for fields radiated by cylindrically symmetric an-tennas in layered media. We also illustrate the application of the algorithm for a particle-in-cell (PIC) simulation of beam-wave interactions inside a high-power backward-waveoscillator. Keywords: body-of-revolution, finite-element time-domain, Maxwell equations,discrete exterior calculus, transformation optics.
1. Introduction
The solution of Maxwell’s equations in circularly symmetric or body-of-revolution(BOR) geometries is important for a plethora of applications involving analysis anddesign of microwave devices (e.g. cavity resonators, coaxial cables, waveguides, antennas,high-power amplifiers, etc.) [1, 2, 3, 4, 5, 6, 7, 8, 9], electromagnetic scattering [10, 11,12, 13], metamaterials [14], and exploration geophysics [15, 16, 17, 18, 19, 20, 21], toname a few. Azimuthal field variations in BOR problems can be described by Fourier
Email addresses: [email protected] (Dong-Yeop Na), [email protected] (Ben-Hur V. Borges), [email protected] (Fernando L. Teixeira)
Preprint submitted to Elsevier September 26, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] S e p odal decomposition, with the modal field solutions reduced to a two-dimensional (2-D) problem in the meridian ρz -plane. Frequency-domain finite element (FE) Maxwellsolvers for BOR problems have been developed in the past by discretizing the second-order vector wave equation using edge elements for either the electric or the magneticfield [6, 7, 12, 14, 22] which avoids some of the pitfalls encountered when using scalarelements [10].It is highly desirable to develop BOR FE solvers in the time domain as well. Time-domain FE solvers are better suited for simulating broadband problems, for capturingtransient processes such as those involved in beam-wave interactions [23, 24, 25], and forhandling non-linear problems. However, the use of the second-order vector wave equationas a starting point for a time-domain FE formulation, as done in frequency-domainMaxwell FE solvers, is inadequate. This is because the vector wave equation admitssolutions of the form t ∇ φ , which are not original solutions of Maxwell’s equations and,even if not excited by (properly set) initial conditions, may emerge in the course of thesimulation due to round-off errors and pollute the results for long integration times [26].To avoid this problem, a mixed (basis) FE solver based directly on the first-order shouldbe adopted in the time domain [27, 28, 29, 30].In this paper, we present a mixed FE BOR solver for time-domain Maxwell’s curlequations based on transformation optics (TO) [31, 32, 33, 34, 35, 36] and discretizationprinciples based on the discrete exterior calculus (DEC) of differential forms [23, 27,37, 38, 39, 40, 41, 42, 43, 44]. We explore TO principles to map the original three-dimensional (3-D) BOR problem to an equivalent problem on the 2-D meridian planewhere the resulting metric is not the cylindrical one but instead the Cartesian one (i.e.,with no radial factors present). The cylindrical metric becomes fully embedded intothe constitutive properties of an effective (artificial) inhomogeneous anisotropic mediumthat fills the entire domain. In this way, a Cartesian 2-D FE code can be retrofittedto this problem with no modifications necessary except to accommodate the presence ofanisotropic media. Similar ideas have been explored in the past but restricted to thefrequency-domain finite-difference (FD) context and to structured grids only [45]. In theFE context considered here, DEC principles are used to discretize Maxwell’s equationson unstructured (irregular) grids using discrete differential (Whitney) forms [33, 37, 40,46, 47]. Unstructured grids permits a more flexible representation of irregular geometriesand reduce the need for geometrical defeaturing. In addition to the above advantages,the proposed formalism facilitates treatment of the coordinate singularity on the axis ofsymmetry ( z axis) because it does not require any modification of the basis functionsfor ρ = 0 (otherwise necessary in prior BOR FE solvers [6, 12, 22]). As detailed in theAppendix, the DEC formalism also facilitates implementation of perfectly matched layers(PML) to truncate the outer boundaries. We validate the algorithm against analyticalsolutions for resonant fields in cylindrical cavities and against pseudo-analytical solutionsfor the radiated fields by cylindrically symmetric antennas in layered media. We alsoillustrate the application of the algorithm to the simulation of wave-beam interactions ina high-power microwave backward-wave oscillator (BWO).2 igure 1: Depiction of an axisymmetric structure.
2. Formulation
Consider a BOR object with symmetry axis along z , such as the waveguide structuredepicted in Fig. 1. It is well known that the vector operators (gradient, curl, and di-vergence) in cylindrical coordinates have additional metric scaling factors not present inCartesian coordinates. However, by exploiting TO concepts [31, 32, 38], we can map thecylindrical-system Maxwell’s curl equations to a Cartesian-like equations where the met-ric factors are embedded into artificial constitutive tensors. For convenience we denotethese calculations under the generic banner of TO but some of these ideas actually pre-date TO per se. They can be traced to earlier applications involving Maxwell’s equationsin BOR geometries and to Weitzenbock identities involving differential forms of differentdegrees [48] in cylindrical (polar) coordinates.Starting from Maxwell’s equations in cylindrical coordinates, and considering artificialanisotropic permittivity and permeability tensors ¯¯ (cid:15) (cid:48) and ¯¯ µ (cid:48) of the form¯¯ (cid:15) (cid:48) = ¯¯ (cid:15) · ¯¯R (cid:15) = ¯¯ (cid:15) · ρ ρ −
00 0 ρ , (1)¯¯ µ (cid:48) = ¯¯ µ · ¯¯R µ = ¯¯ µ · ρ − ρ
00 0 ρ − , (2)where the constitutive parameters of the original medium are given by¯¯ (cid:15) = (cid:15) ρ (cid:15) φ
00 0 (cid:15) z , ¯¯ µ = µ ρ µ φ
00 0 µ z . E (cid:48) = ¯¯R E · E = ρ
00 0 1 · E , (3) D (cid:48) = ¯¯R D · D = ρ ρ · D , (4) B (cid:48) = ¯¯R B · B = ρ ρ · B , (5) H (cid:48) = ¯¯R H · H = ρ
00 0 1 · H , (6)we can rewrite the resulting Maxwell’s curl equations as ∇ (cid:48) × E (cid:48) = − ∂ B (cid:48) ∂t , (7) ∇ (cid:48) × H (cid:48) = ∂ D (cid:48) ∂t , (8) D (cid:48) = ¯¯ (cid:15) (cid:48) · E (cid:48) , (9) B (cid:48) = ¯¯ µ (cid:48) · H (cid:48) , (10)with ∇ (cid:48) × A (cid:48) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ ρ ˆ φ ˆ z ∂∂ρ ∂∂φ ∂∂z A (cid:48) ρ A (cid:48) φ A (cid:48) z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (11)The modified curl operator in the equivalent (primed) system seen in (11) is devoid ofany radial scaling and thus locally isomorphic to the Cartesian curl operator.4 a) (b)Figure 2: (2+1) setup for fields on (a) primal and (b) dual meshes at the meridian plane. The verticalaxis is ρ and the horizontal axis is z . We decompose the fields into two sets: TE φ - and TM φ -polarized fields, correspondingto { E (cid:48) ρ , E (cid:48) z , B (cid:48) φ } and { E (cid:48) φ , B (cid:48) ρ , B (cid:48) z } , respectively. In what follows, we use superscripts (cid:107) or ⊥ to denote fields transverse or normal to the 2-D meridian plane. The TE φ fieldcomponents can be expressed as E (cid:48)(cid:107) and B (cid:48)⊥ and the TM φ as E (cid:48)⊥ and B (cid:48)(cid:107) . In the DECcontext, the electric field intensity, the magnetic flux density, the electric flux density,and the magnetic field intensity are likewise represented as 1-, 2-, 2-, and 1-forms on the3-D Euclidean space, respectively [38]. For present analysis based on the meridian plane(a 2-D manifold), E (cid:107) is transverse to the plane and still is represented as a 1-form. Onthe other hand, E ⊥ should be represented as a 0-form since it is a point-based quantityon this manifold. Likewise, although B ⊥ is a 2-form in 3-D, B (cid:107) is represented as a 1-formon the 2-D meridian plane (see Fig. 2).
1- and 2-forms correspond to physical quantities naturally associated to line and surface integrals,respectively. .3. Mixed FE time-domain BOR solver We factor the transverse (i.e. ρ and z ) and normal (i.e. φ ) variations of the polarization-decomposed Maxwell fields on the 2-D meridian plane as E (cid:48) ( ρ, φ, z, t ) = M φ (cid:88) m = − M φ E (cid:48)(cid:107) m ( ρ, z, t ) Φ m ( φ ) + M φ (cid:88) m = − M φ E (cid:48)⊥ m ( ρ, z, t ) Ψ m ( φ ) , (12) B (cid:48) ( ρ, φ, z, t ) = M φ (cid:88) m = − M φ B (cid:48)⊥ m ( ρ, z, t ) Φ m ( φ ) + M φ (cid:88) m = − M φ B (cid:48)(cid:107) m ( ρ, z, t ) Ψ m ( φ ) , (13)where M φ is the maximum order of the Fourier harmonics considered andΦ m ( φ ) = cos ( mφ ) , for m < , for m = 0sin ( mφ ) , for m > , (14)Ψ m ( φ ) = sin ( mφ ) , for m < , for m = 0cos ( mφ ) , for m > . (15)Substituting (12) and (13) into (7), by using the orthogonality between modes, i.e. (cid:90) π Φ m ( φ ) Φ n ( φ ) dφ = C m δ mn , (16) (cid:90) π Ψ m ( φ ) Ψ n ( φ ) dφ = C m δ mn , (17)where C m = π for m (cid:54) = 0 and C = 2 π , we obtain the modal Faraday’s law as ∇ (cid:48)(cid:107) × E (cid:48)(cid:107) m ( ρ, z, t ) = − ∂ B (cid:48)⊥ m ( ρ, z, t ) ∂t , (18) ∇ (cid:48)(cid:107) × E (cid:48)⊥ m ( ρ, z, t ) = − ∂ B (cid:48)(cid:107) m ( ρ, z, t ) ∂t + | m | E (cid:48)(cid:107) m ( ρ, z, t ) × ˆ φ, (19)for m = − M φ , ..., M φ , where ∇ (cid:48)(cid:107) = ˆ ρ∂/∂ρ + ˆ z∂/∂z .We discretize (18) and (19) on the meridian plane using an unstructured mesh basedon simplicial (triangular) cells and by expanding the fields in a mixed basis as scalar orvector proxies of discrete differential forms (Whitney forms) [27, 38, 43]. In particular,the TE φ field is expanded as E (cid:48)(cid:107) m ( ρ, z, t ) = N (cid:88) j =1 E (cid:107) j,m ( t ) W (1) j ( ρ, z ) , (20) B (cid:48)⊥ m ( ρ, z, t ) = N (cid:88) k =1 B ⊥ k,m ( t ) W (2) k ( ρ, z ) , (21)where W ( p ) q is the vector proxy of a Whitney p -form w ( p ) q [24] associated with the q -th p -cell ( p = 0 , , N p is the6otal number of p -cells on the grid. The expressions for the Whitney forms and theirproxies are provided in Appendix A. Likewise, the TM φ field is represented as E (cid:48)⊥ m ( ρ, z, t ) = N (cid:88) i =1 E ⊥ i,m ( t ) ˆ φ W (0) i ( ρ, z ) , (22) B (cid:48)(cid:107) m ( ρ, z, t )= N (cid:88) j =1 B (cid:107) j,m ( t ) W (RWG) j ( ρ, z ) . (23)In what follows, we denote W (1) j × ˆ φ = W (RWG) j , since this expression recovers the so-called Rao-Wilton-Glisson (RWG) functions [49, 50] . Note that we use dummy indexsubscripts i , j , and k to indicate the i -th node, j -th edge, and k -th face, respectively.The various basis functions above are depicted in Fig. 3, see also [52, 53].By substituting (20) and (21) into (18), and (22) and (23) into (19), we obtain thefollowing equations N (cid:88) j =1 E (cid:107) j,m ( t ) (cid:16) ∇ (cid:48)(cid:107) × W (1) j (cid:17) = − ∂∂t N (cid:88) k =1 B ⊥ k,m ( t ) W (2) k (24) N (cid:88) i =1 E ⊥ i,m ( t ) ∇ (cid:48)(cid:107) W (0) i = − ∂∂t N (cid:88) j =1 B (cid:107) j,m ( t ) W (1) j + | m | N (cid:88) j =1 E (cid:107) j,m ( t ) W (1) j , (25)for m = − M φ , ..., M φ and where we have used the fact that ∇ (cid:48)(cid:107) × (cid:16) ˆ φ W (0) i (cid:17) = (cid:16) ∇ (cid:48)(cid:107) W (0) i (cid:17) × ˆ φ . The equations above can be recast using the exterior calculus of differential forms as N (cid:88) j =1 E (cid:107) j,m ( t ) (cid:16) d (cid:48)(cid:107) w (1) j (cid:17) = − ∂∂t N (cid:88) k =1 B ⊥ k,m ( t ) w (2) k , (26) N (cid:88) i =1 E ⊥ i,m ( t ) (cid:16) d (cid:48)(cid:107) w (0) i (cid:17) = − ∂∂t N (cid:88) j =1 B (cid:107) j,m ( t ) w (1) j + | m | N (cid:88) j =1 E (cid:107) j,m ( t ) w (1) j , (27)where d (cid:48)(cid:107) = dρ ∂/∂ρ + dz ∂/∂z is the exterior derivative on the meridian plane.Applying DEC principles, (26) can be paired to the 2-cells of the mesh and (27) tothe 1-cells of the mesh (see Appendix B) so that, by invoking the generalized Stokes’theorem [27, 38, 40, 43, 44] (see Appendix C), the exterior derivative can be replaced byincidence operators on the mesh (see also Appendix D). Next, by discretizing the timederivatives using central-differences in a staggered manner (leap-frog time discretization)we obtain the following update equations for Faraday’s law (cid:2) B ⊥ m (cid:3) n + = (cid:2) B ⊥ m (cid:3) n − − ∆ t [ D curl ] · (cid:104) E (cid:107) m (cid:105) n , (28) (cid:104) B (cid:107) m (cid:105) n + = (cid:104) B (cid:107) m (cid:105) n − − ∆ t (cid:16) [ D grad ] · (cid:2) E ⊥ m (cid:3) n − | m | (cid:104) E (cid:107) m (cid:105) n (cid:17) , (29) In other words, W (RWG) j is the Hodge dual of W (1) j in 2-D [40, 43, 51]. a) (b)(c) (d)Figure 3: Vector proxies of various degrees of Whitney forms on the mesh: (a) W (1) j , (b) W (2) k , (c) W (0) i ,and (d) W (RWG) j . Note that t j is a unit vector tangential to j − th edge and parallel to its direction and n k is a unit vector normal to k − th face. where ∆ t is a time step increment and the superscript n indicates the time-step index.[ D curl ] and [ D grad ] are N × N and N × N incidence matrices, respectively, that encodethe curl and the gradient operators on the FE mesh with elements in the set {− , , } (see Appendix D). The field unknowns are represented by the column vectors (cid:2) B ⊥ m (cid:3) = (cid:2) B ⊥ m, , ..., B ⊥ m,N (cid:3) T , (cid:104) E (cid:107) m (cid:105) = (cid:104) E (cid:107) m, , ..., E (cid:107) m,N (cid:105) T , (cid:104) B (cid:107) m (cid:105) = (cid:104) B (cid:107) m, , ..., B (cid:107) m,N (cid:105) T , and (cid:2) E ⊥ m (cid:3) = (cid:2) E ⊥ m, , ..., E ⊥ m,N (cid:3) T .We proceed along similar lines for Ampere’s law by expressing the D (cid:48) and H (cid:48) fieldsas D (cid:48) ( ρ, φ, z, t ) = M φ (cid:88) m =0 D (cid:48)(cid:107) m ( ρ, z, t ) Φ m ( φ ) + M φ (cid:88) m =0 D (cid:48)⊥ m ( ρ, z, t ) Ψ m ( φ ) , (30) H (cid:48) ( ρ, φ, z, t ) = M φ (cid:88) m =0 H (cid:48)⊥ m ( ρ, z, t ) Φ m ( φ ) + M φ (cid:88) m =0 H (cid:48)(cid:107) m ( ρ, z, t ) Ψ m ( φ ) . (31)8fter substituting (30) and (31) to (8), applying trigonometric orthogonality to theresulting equations, and matching the field components, we arrive at ∇ (cid:48)(cid:107) × H (cid:48)(cid:107) m ( ρ, z, t ) = ∂ D (cid:48)⊥ m ( ρ, z, t ) ∂t , (32) ∇ (cid:48)(cid:107) × H (cid:48)⊥ m ( ρ, z, t )= ∂ D (cid:48)(cid:107) m ( ρ, z, t ) ∂t − | m | H (cid:48)(cid:107) m ( ρ, z, t ) × ˆ φ. (33)As before, we discretize (32) and (33) on the 2-D meridian plane, the important differencebeing that the discretization for D (cid:48) and H (cid:48) is on the dual mesh [33, 38, 43, 51], as opposedto the FE (primal) mesh as done for E (cid:48) and B (cid:48) . In this way, we obtain D (cid:48)(cid:107) m ( ρ, z, t ) = ˜ N (cid:88) j =1 D (cid:107) j,m ( t ) ˜ W (RWG) j ( ρ, z ) , (34) H (cid:48)⊥ m ( ρ, z, t ) = ˜ N (cid:88) i =1 H ⊥ i,m ( t ) ˆ φ ˜W (0) i ( ρ, z ) , (35) D (cid:48)⊥ m ( ρ, z, t ) = ˜ N (cid:88) k =1 D ⊥ k,m ( t ) ˜ W (2) k ( ρ, z ) , (36) H (cid:48)(cid:107) m ( ρ, z, t ) = ˜ N (cid:88) j =1 H (cid:107) j,m ( t ) ˜ W (1) j ( ρ, z ) . (37)where we use the tilde˜to denote quantities associated with the dual mesh. Similar to thediscrete counterparts of Faraday’s law, by substituting (34) and (35) into (32) and (36)and (37) into (33) and by applying DEC principles and a leap-frog time discretization tothe resulting equations, we obtain the discrete representations of Ampere’s law as (cid:2) D ⊥ m (cid:3) n +1 = (cid:2) D ⊥ m (cid:3) n + ∆ t (cid:104) ˜ D curl (cid:105) · (cid:104) H (cid:107) m (cid:105) n + , (38) (cid:104) D (cid:107) m (cid:105) n +1 = (cid:104) D (cid:107) m (cid:105) n + ∆ t (cid:18)(cid:104) ˜ D grad (cid:105) · (cid:2) H ⊥ m (cid:3) n + − | m | (cid:104) H (cid:107) m (cid:105) n + (cid:19) , (39)where (cid:104) ˜ D curl (cid:105) and (cid:104) ˜ D grad (cid:105) are incidence matrices on the dual mesh, with sizes ˜ N × ˜ N and ˜ N × ˜ N , respectively. As before, (cid:2) H ⊥ m (cid:3) , (cid:104) D (cid:107) m (cid:105) , (cid:104) H (cid:107) m (cid:105) , and (cid:2) D ⊥ m (cid:3) are column vectorscontaining the degrees of freedom of the modal fields.We use the (discrete) Hodge star operator [33, 38, 43, 51] to convert the discreteAmpere’s law from the dual mesh to the primal mesh. In this way,[ (cid:63) (cid:15) ] → · (cid:2) E ⊥ m (cid:3) n +1 = [ (cid:63) (cid:15) ] → · (cid:2) E ⊥ m (cid:3) n + ∆ t (cid:18) [ D grad ] T · (cid:2) (cid:63) µ − (cid:3) → · (cid:104) B (cid:107) m (cid:105) n + (cid:19) , (40)[ (cid:63) (cid:15) ] → · (cid:104) E (cid:107) m (cid:105) n +1 = [ (cid:63) (cid:15) ] → · (cid:104) E (cid:107) m (cid:105) n +∆ t (cid:18) [ D curl ] T · (cid:2) (cid:63) µ − (cid:3) → · (cid:2) B ⊥ m (cid:3) n + − | m | (cid:2) (cid:63) µ − (cid:3) → · (cid:104) B (cid:107) m (cid:105) n + (cid:19) , (41)9here (cid:104) ˜ D curl (cid:105) = [ D grad ] T , (cid:104) ˜ D grad (cid:105) = [ D curl ] T and the discrete Hodge matrix elementsare given by [38, 43, 47][ (cid:63) (cid:15) ] → J,j = (cid:90) Ω ( (cid:15) ρ ) w (1) J ∧ (cid:63) (cid:16) w (1) j (cid:17) = (cid:90) Ω ( (cid:15) ρ ) W (1) J · W (1) j dV (cid:124) (cid:123)(cid:122) (cid:125) vector proxy representation , (42) (cid:2) (cid:63) µ − (cid:3) → K,k = (cid:90) Ω (cid:0) µ − ρ (cid:1) w (2) K ∧ (cid:63) (cid:16) w (2) k (cid:17) = (cid:90) Ω (cid:0) µ − ρ (cid:1) W (2) K · W (2) k dV (cid:124) (cid:123)(cid:122) (cid:125) vector proxy rep. , (43)[ (cid:63) (cid:15) ] → I,i = (cid:90) Ω (cid:0) (cid:15) ρ − (cid:1) w (0) I ∧ (cid:63) (cid:16) w (0) i (cid:17) = (cid:90) Ω (cid:0) (cid:15) ρ − (cid:1) (cid:104) W (0) I ˆ φ (cid:105) · (cid:104) W (0) i ˆ φ (cid:105) dV (cid:124) (cid:123)(cid:122) (cid:125) vector proxy rep. , (44) (cid:2) (cid:63) µ − (cid:3) → J,j = (cid:90) Ω ( µ ρ ) − w (RWG) J ∧ (cid:63) (cid:16) w (RWG) j (cid:17) = (cid:90) Ω ( µ ρ ) − (cid:104) W (1) J × ˆ φ (cid:105) · (cid:104) W (1) j × ˆ φ (cid:105) dV (cid:124) (cid:123)(cid:122) (cid:125) vector proxy rep. , (45)where Ω is the (compact) spatial support of the Whitney forms, and the ρ , ρ − factorsresult from the use of the TO in the mapping, as discussed before, where they enteras modifiers of constitutive properties rather than differential operator factors. Thediscrete Hodge matrices defined in (42), (43), (44), and (45) are instantiations of the(discrete) Galerkin-Hodge operator. It should be emphasized that the Galerkin-Hodgeoperator is not a natural consequence of DEC. The Galerkin-Hodge operator was origi-nally proposed in [54]. It satisfies a number of built-in properties for stability in arbitrarysimplicial meshes as discussed, for example, in references [43],[55],[56],[57]. In particular,the Galerkin-Hodge operator enforces standard local energy positivity [42].The field updates in (40) and (41) call for sparse linear solvers due to the presenceof the matrices [ (cid:63) (cid:15) ] → and [ (cid:63) (cid:15) ] → . From (42) and (44), it is seen that [ (cid:63) (cid:15) ] → and[ (cid:63) (cid:15) ] → are diagonally dominant and symmetric positive definite matrices; consequently,the linear solve can be performed very quickly. Nevertheless, this needs to be repeatedat every time step. The linear solve can be obviated by computing a sparse approximateinverse (SPAI) of [ (cid:63) (cid:15) ] → and [ (cid:63) (cid:15) ] → prior to the start of the time updating procedure.This strategy is discussed in [25] and [41]. The present algorithm is explicit and henceconditionally stable. The stability conditions are discussed in Appendix G. For BOR problems where the line ρ = 0 (symmetry axis) is part of the solutiondomain (for example, in hollow waveguides), it becomes necessary to treat the fieldbehavior there by means of appropriate boundary conditions. The boundary conditionsat ρ = 0 are mode-dependent and should account for the cylindrical coordinate systemsingularity and the related degeneracy of the ˆ ρ and ˆ φ unit vectors there. When m = 0,there is no field variation along azimuth and, in the absence of charges at ρ = 0, bothazimuthal and radial field components are zero at ρ = 0. On the other hand, the axialfield component should be zero for m (cid:54) = 0 [58] since the axial direction is invariant with10 a)(b)(c)(d)Figure 4: Field boundary conditions on the primal mesh for the TE φ field with (a) perfect magneticconductor ( m = 0) and (b) perfect electric conductor ( m (cid:54) = 0) and for the TM φ field with (c) perfectmagnetic conductor ( m (cid:54) = 0) and (d) perfect electric conductor ( m = 0). Dashed lines indicate Dirichletboundary condition, for example edges on the z axis representing a perfect electric conductor boundaryfor TE φ field in (b), or nodes on the z axis representing a perfect electric conductor boundary for theTM φ field in (d). φ and a field dependency of the form cos ( mφ ) or sin ( mφ ) with m (cid:54) = 0 wouldimply a multivalued result at ρ = 0 due to the coordinate degeneracy there. As a result,when m = 0, the boundary ρ = 0 can be represented as a perfect electric conductorfor the TE φ field and as a perfect magnetic conductor for the TM φ field. Conversely,when m (cid:54) = 0, the ρ = 0 boundary can be represented as a perfect magnetic conductorfor the TE φ field and as a perfect electric conductor for the TM φ field. A homogeneousNeumann boundary condition for the electric field can be used to represent the perfectmagnetic conductor case and a homogeneous Dirichlet boundary condition for the perfectelectric conductor case. Implementation of such boundary conditions on the primal meshis illustrated in Fig. 4. Dashed lines in Fig. 4b and 4d denote the Dirichlet boundaryimplementation: along the z axis, the perfect electric conductor condition is enforcedon grid edges for the TE φ case and on grid nodes for the TM φ case. Likewise, Fig. 4aand 4c illustrate application of the Neumann boundary condition: along the z axis, theperfect magnetic conductor condition is enforced on grid edges for the TE φ case and ongrid nodes for the TM φ case.Using the boundary conditions described above, the present FETD-BOR Maxwellsolver does not require any modifications in the basis functions on the grid cells adjacentto the z axis, unlike prior FE-BOR Maxwell solvers.
3. Numerical Examples
In order to validate present FETD-BOR Maxwell solver, we first consider a cylindricalcavity and compare the resonance frequency results to the analytical predictions. Then,we illustrate two practical examples of devices based on BOR geometries: logging-while-drilling sensors used for Earth formation resistivity profiling in geophysical explorationand relativistic BWO for high-power microwave applications.
We simulate the eigenfrequencies of a hollow cylindrical cavity with metallic wallsusing the present FETD-BOR Maxwell solver, and compare the results to analytic pre-dictions. The cavity has radius a = 0 . h = 1 m, as depicted in Fig. 5.Magnetic and electric dipole current sources M ( r , t ) and J ( r , t ) oriented along φ andexcited by broadband Gaussian-modulated pulses are placed at arbitrary locations insidethe cavity r s = ( ρ s , φ s , z s ), so that M ( r , t ) , J ( r , t ) = ˆ φ G ( t ) δ ( r − r s ) == ˆ φ G ( t ) δ (cid:16) r (cid:107) − r (cid:107) s (cid:17) π + 2 π M φ (cid:88) m =1 cos [ m ( φ − φ s )] (46)where G ( t ) = e − [( t − t g ) / (2 σ g )] sin [2 πf g ( t − t g )] with t g = 20 ns, σ g = 1 . f g = 300MHz, and r (cid:107) = ρ ˆ ρ + z ˆ z . We use Fourier series expansion to describe δ ( φ − φ s ) in(46) in order to match the modal field expansion used before. A total of four dipolesources (electric and magnetic currents) are used to excite a rich gamut of eigenmodes,as illustrated in Fig. 5. The meridian plane of the cylindrical cavity is discretized byan unstructured grid with 4 ,
045 nodes, 11 ,
939 edges, and 7 ,
895 faces (seen as the ρz igure 5: Schematic view of the simulated cylindrical cavity with perfect electric conductor (PEC) walls.The cavity dimensions are a = 0 . h = 1 m.Table 1: Maximum time-step intervals for various cases in the simulation of cylindrical metallic cavity. m = 0 m (cid:54) = 0TE φ -pol. TM φ -pol. m = 1 m = 2 m = 3 m = 4∆ t max [ps] 10.009 10.249 10.009 6.4792 4.5545 3.4843plane for φ = 180 o in Fig. 7). The metallic boundaries are treated as perfect electricconductors. In this case, the maximum azimuthal modal order M φ was set equal to 4 toinvestigate the field solution up to this order. Higher order modes can be included bysimply increasing M φ . This is straightforward since azimuthal modal fields with differentorders are orthogonal to each other. From the stability analysis in Appendix G, themaximum time-step intervals for various cases are presented in Table 1. Here we chose∆ t = 1 ps for the simulations and used a total of 1 × time steps to provide sufficientlynarrow resonance peaks. By recording the time history of the electric field values atarbitrary locations inside the cavity and performing a Fourier transform, we obtain theeigenfrequencies as peaks in the Fourier spectrum. Fig. 6 shows the normalized spectralamplitude as a function of frequency. The black solid line is the result obtained byusing present FETD-BOR Maxwell solver. The red dashed and blue solid lines indicateanalytic predictions for the eigenfrequencies of the TE mnp and TM mnp modes in this13 igure 6: Normalized spectral amplitude for E , showing the eigenfrequencies of the cavity. Blacksolid lines correspond to the present FETD-BOR result. Red solid and blue dashed lines are analyticpredictions for the TE mnp and TM mnp eigenfrequencies, respectively. cavity, respectively. The analytic expressions for the eigenfrequencies are given by f TE mnp = 2 cπ (cid:114) χ (cid:48) mn + (cid:16) pπh (cid:17) , for m = 0 , , ..., n = 1 , , ..., p = 1 , , ... , (47) f TM mnp = 2 cπ (cid:114) χ mn + (cid:16) pπh (cid:17) , for m = 0 , , ..., n = 1 , , ..., p = 0 , , ... , (48)where c is speed of light, χ mn and χ (cid:48) mn are the roots of the equations J m ( aχ mn ) = 0and J (cid:48) m ( aχ (cid:48) mn ) = 0, respectively, with J m ( · ) being the Bessel function of first kind and J (cid:48) m ( · ) its derivative with respect to the argument. It is clear from Fig. 6 that there isa great agreement between the simulated and analytic eigenfrequencies. Table 2 showsthe relative error between the simulated f s and analytical f a frequencies. The relativeerror is below 0 .
03 % in all cases, indicating the accuracy of the proposed field solver.To illustrate the field behavior, Figs. 7 and 8 show snapshots for electric field intensityand magnetic flux density distribution inside the cavity on four ρz planes with φ =0 o , φ = 90 o , 180 o , 270 o and two ρφ planes with z = 0 . . . µ s, 1 . µ s, 1 . µ s, and 1 . µ s. Due to the location of thedipole sources, the transient fields produced include many eigenmodes, and are basicallyasymmetric. It can be seen that the (tangential or normal) boundary conditions on theouter perfect electric conductor walls for electric field intensity and magnetic flux densityare well satisfied. Moreover, the correct field distribution along the symmetry axis is wellreproduced by the chosen boundary conditions at ρ = 0, without any spurious artifacts.14 a) (b)(c) (d)Figure 7: Transient snapshots for E z inside the cylindrical cavity at (a) 1 . µ s], (b) 1 . µ s], (c)1 . µ s], and (d) 1 . µ s]. a) (b)(c) (d)Figure 8: Transient snapshots for B z inside the cylindrical cavity at (a) 1 . µ s], (b) 1 . µ s], (c)1 . µ s], and (d) 1 . µ s]. able 2: Eigenfrequencies for the cylindrical cavity and normalized errors between numerical and analyticresults. Resonant modes f a [MHz] | f a − f s | /f a ×
100 [%]TM . . × − TE . . × − TM . . × − TE . . × − TE . . × − TM . . × − TM . . × − TE , TM . . × − TE . . × − TE . . × − TE , TM . . × − TE . . × − TM . . × − TE . . × − TM . . × − TM . . × − TM . . × − TE . . × − TE . . × − TE . . × − TM . . × − Logging-while-drilling sensors have BOR geometries and are routinely used for hy-drocarbon exploration [16, 17, 18, 19, 20, 21]. As the drilling process is performed, thesesensors record logs obtained by the measurements of fields produced by loop (multi-coil)antennas present in the sensor and reflected from the surrounding geological formation.Logging-while-drilling sensors are typically equipped with a series of transmitter and re-ceiver loop antennas that are wrapped around the outer diameter of a metallic mandrelattached to the bit drill [59, 60, 61, 62, 63, 64]. Fields produced by the transmitter coil(s)interact with the adjacent well-bore environment and are detected by a pair (or more)of receiver coils along the logging-while-drilling sensor at same axial distance from thetransmitter(s). Two types of measurements are typically used to determine the resis-tivity profiles of the adjacent formation. The first is the amplitude ratio (AR) betweenthe electromotive force (e.m.f.) excited at the two receiver coils and the second is theirphase difference (PD). In this section, we consider a prototypical concentric logging-while-drilling sensor generating a TM φ field distribution in the formation with m = 0 . Thelogging-while-drilling sensor depicted in Fig. 9 consists of a metallic cylindrical mandrelmodeled as a perfect electric conductor inside a concentric cylindrical borehole. Three Not only the geometry but also the field excitation is axisymmetric in this case. igure 9: Logging-while-drilling sensor problem geometry (from inner to outer features): metallic man-drel, transmit (Tx) and receive (Rx) coil antennas, mud-filled borehole, and adjacent geological forma-tion. loop antennas are used: one as transmitter and two as receivers. The borehole createdby the drilling process is filled with a lubricant fluid (mud). The three coil antennas aremoving downward in tandem as the drilling process occur.We consider two scenarios for the adjacent Earth formation, as shown in Fig. 10.In the first scenario, the borehole is filled with a low conductive (oil-based) fluid (mud)having σ = 0 . σ = 2 S/m, and the formation has three horizontal layers with different con-ductivities as shown. We compute the AR and PD as the set of coil antennas (sensor)moves downward. In both cases, the relative permittivity and permeability are assumedequal to one everywhere, and the transmitter coil radiates a 2 MHz signal. In the timedomain, this is implemented through a current signal along the transmitter coil given by I Tx ( t ) = r ( t ) sin ( ωt ), where r ( t ) = , t < . (cid:20) − cos (cid:18) ωt α (cid:19)(cid:21) , (cid:54) t < αT , t (cid:62) αT, (49)is a raised-cosine ramp function, T = 2 π/ω is the signal period, and α is the number ofsine wave cycles during the ramp duration αT . The use of ramp function mitigates highfrequency components otherwise produced by an abrupt turn-on at t = 0, and yieldsfaster convergence of AR and PD (after approximately one time period T ) [21]. We18 a) (b)Figure 10: Logging-while-drilling sensor responses. (a) First scenario: the conductivity of the adjacentgeological formation is varied. (b) Second scenario: the sensor moves downward through a boreholesurrounded by a geological formation with three horizontal layers. choose α = 0 . θ and amplitudes A using θ = tan − (cid:18) q sin ( ωt ) − q sin ( ωt ) q cos ( ωt ) − q cos ( ωt ) (cid:19) , (50) A = (cid:12)(cid:12)(cid:12)(cid:12) q sin ( ωt + θ ) (cid:12)(cid:12)(cid:12)(cid:12) , (51)where q and q are signals computed at times t and t , respectively [21]. Next, the ARand PD are calculated as AR = A Rx /A Rx , (52)PD = θ Rx − θ Rx . (53)The azimuthal electric current along the transmitter coil is modeled as a nodal currentdensity on the meridian plane and the metallic mandrel is regarded as perfect electricconductor. The FE domain is truncated by a PML to mimic an open domain. We use 8layers for the PML to yield a reflectance below −
50 dB [30].Fig. 11 shows results for the behavior of AR and PD versus the conductivity ona homogeneous formation. The results are compared against previous results obtainedby the finite-difference time-domain (FDTD) and the numerical mode matching (NMM)methods [21]. There is excellent agreement between the results. Results for the secondscenario are shown in Fig. 12, where PD is plotted as a function of the z position ofthe transmitter, z Tx , and compared against previous results obtained by the FDTD andNMM methods [21]. Again, an excellent agreement is obtained. As expected, the PD19s higher when the coil antennas are within high attenuation (high conductivity) layerand vice versa. The conductance profile and the corresponding axial extension of eachformation is shown in green color in Fig. 12. Fig. 13a − Fig. 13f show snapshots of theelectric field distributions for different z Tx to illustrate the field behavior. In this section, we consider a backward-wave oscillator (BWO) driven by energeticelectron beams in the relativistic regime designed to produce a high-power microwavesignal [65], as depicted in Fig. 14. The proposed FETD-BOR solver is incorporated into aPIC algorithm [66, 67, 68] to simulate the wave-plasma interaction in the device [9]. ThePIC algorithm is based on an unstructured grid and explained in detail in [9, 24, 25]. Forthis problem it suffices to consider the TE φ polarized field with m = 0. In a relativisticBWO, the energy of space-charge modes is converted into microwaves via Cerenkovradiation [69]. The BWO employs a slow-wave structure to produce such radiation [70].In present case, the BWO system consists of a cathode, an anode, a slow-wave structurewith sinusoidal corrugations, a beam collector, and a coaxial output port, as depictedin Fig. 15a. The electron beam is produced by an external voltage between the cathodeand anode. In the slow-wave structure, the space charge modes evolve to TM modes.The oscillation of the modal field leads to beam velocity modulation and a quasi-periodicbunching of the electron beam distribution. Lateral beam confinement is obtained by anexternally applied static magnetic field. Coherent RF signals are detected and extractedat the coaxial RF output port as illustrated in Fig. 15a. The outer radial boundary ofthe slow-wave structure is expressed as R ( z ) = ( A − B ) cos (2 πz/C ) + B where A and B are maximum and minimum radii, respectively, and C is the axial corrugation period.For X-band operation, we set A = 1 .
95 cm, B = 1 .
05 cm, and C = 1 .
67 cm for a beam-velocity v = 2 . × m/s. The coaxial RF output port is truncated by a PML [29, 30].The unstructured mesh has N = 2 , N = 8 , N = 5 , l ave =1.4468 mmwhere l ave is an average edge size. The time step is ∆ t = 0 . .
5. As typical in PIC simulations, we employ acoarse-graining of the phase space, and each “superparticle” in the simulation represents1 . × electrons. The resultant electron density n e yields a Debye length λ D = 20 . N D equals 5 . × particlesand hence a collisionless plasma assumption is valid in this case. The self-field evolutionand spectrum at the output port are shown in Figs. 16a and 16b, respectively. FromFig. 16a, it is seen that the field grows to an RF oscillation near 50 ns and saturates ataround 100 ns. The output signal has a peak at 8 .
27 GHz, as shown in Fig. 16b. Fig.15a shows the electron beam distribution at 43 . . mode is indeedpresent.
4. Conclusion
We presented a new finite-element time-domain (FETD) Maxwell solver for the anal-ysis of body-of-revolution (BOR) geometries. The proposed solver is based on discrete20xterior calculus (DEC) and transformation optics (TO) concepts. We explored TOprinciples to map the original 3-D problem from a cylindrical coordinate system to anequivalent problem on a 2-D (Cartesian-like) meridian ρz plane, where the cylindricalmetric is factored out from the differential operators and embedded on an effective (ar-tificial) inhomogeneous and anisotropic medium that fills the domain. This enables theuse of Cartesian 2-D FE code with no modifications necessary except to accommodatethe presence of anisotropic media. The spatial discretization is done on an unstructuredmesh on the 2-D meridian plane and effected by decomposing the fields into their TE φ and TM φ components and expanding each eigenmode into an appropriate set of (vectoror scalar) basis functions (Whitney forms) based on DEC principles. A leap-frog (sym-plectic) time-integrator is applied to the semi-discrete Maxwell curl equations and usedto obtain a fully discrete, marching-on-time evolution algorithm. Unlike prior solvers,the present FETD-BOR Maxwell solver does not require any modifications on the basisfunctions adjacent to the symmetry axis. Rather, the field behavior on the symmetryaxis can be simply implemented through properly selected homogeneous Dirichlet andNeumann applied to the eigenmodal expansion. Acknowledgments
This work was supported in part by National Science Foundation grant ECCS-1305838, Department of Defense, Defense Threat Reduction Agency grant HDTRA1-18-1-0050, Ohio Supercomputer Center grants PAS-0061 and PAS-0110, S˜ao Paulo StateResearch Foundation (FAPESP) grant 2015/50268-5, and the Ohio State University Pres-idential Fellowship program.The content of the information does not necessarily reflect the position or the policyof the U.S. federal government, and no official endorsement should be inferred.21 a)(b)Figure 11: Computed (a) AR and (b) PD (in deg.) by a logging-while-drilling sensor surrounded byhomogeneous geological formations with different conductivities. This corresponds to the first scenarioin Fig. 10. The results from the present algorithm are compared against FDTD and NMM results [21](see more details in the main text). igure 12: Computed PD (deg.) between the two receivers of the logging-while-drilling sensor versusthe z position of the transmitter coil antenna. This corresponds to the second scenario in Fig. 10. Theresults from the present algorithm are compared against FDTD and NMM results [21] (see more detailsin the main text). a) (b)(c) (d)(e) (f)Figure 13: Electric field distribution during the half period for z Tx = (a) −
50 inch, (b) −
25 inch, (c) 5inch, (d) 25 inch, (e) 50, and (f) 70 inch. Note that z Tx = 0 at the interface between first (5 S/m) andsecond (0 . igure 14: Relativistic backward-wave oscillator with a sinusoidally-corrugated slow-wave structuredriven by a relativistic electron beam. (a)(b)Figure 15: Snapshots for (a) the velocity-modulated electron beam at 43 . . ρ and horizontal axis is z . a)(b)Figure 16: Output signals from the BWO device in (a) time domain and (b) frequency domain. ppendix A. Whitney forms and pairing operations Whitney p -forms are canonical interpolants of discrete differential p -forms [ ? ]. Asexplained below, Whitney p -forms are naturally paired to the p -cells of the mesh, where p refers to the dimensionality, i.e. p = 0 refers to nodes, p = 1 to edges, p = 1 to facetsand so on [38]. On simplices (e.g. on triangular cells in 2-D or tetrahedral cells in 3-D),Whitney 0-, 1-, and 2-forms are expressed as [38, 52 ? ] w (0) i = λ i , (A.1) w (1) i = λ i a dλ i b − λ i b dλ i a , (A.2) w (2) i = 2 ( λ i a dλ i b ∧ dλ i c + λ i b dλ i c ∧ dλ i a + λ i c dλ i a ∧ dλ i b ) , (A.3)where d is the exterior derivative, ∧ is the exterior product, i a , i b , and i c denote the gridnodes belonging to the i -th p -cell for p = 1 or 2, and λ denotes the barycentric coordinateassociated to a given node.The corresponding vector proxies for Whitney 0-, 1-, and 2-forms write as [24, 38]W (0) i = λ i , (A.4) W (1) i = λ i a ∇ λ i b − λ i b ∇ λ i a , (A.5) W (2) i = 2 ( λ i a ∇ λ i b × ∇ λ i c + λ i b ∇ λ i c × ∇ λ i a + λ i c ∇ λ i a × ∇ λ i b ) , (A.6)One of the key properties of Whitney p -forms is that they admit a natural “pairing”with the p -cells of the mesh [38]. Computationally, the pairing operation between an i -th p -cell of the grid σ i ( p ) and a Whitney form w ( p ) j associated with the j -th p -cell is effectedby the integral below and yields [38, 43] (cid:68) σ i ( p ) , w ( p ) j (cid:69) = (cid:90) σ i ( p ) w ( p ) j = δ i,j , (A.7)where δ i,j is the Kronecker delta, for p = 0 , . . . , Appendix B. Generalized Stokes’ theorem
The generalized Stokes’ theorem of exterior calculus [38, 43, 40, 71, 72] states (cid:68) σ ( p +1) , dw ( p ) j (cid:69) = (cid:68)(cid:0) ∂σ ( p +1) (cid:1) ( p ) , w ( p ) j (cid:69) (B.1)where ∂ is the boundary operator that maps an (oriented) p -cell on the grid to the setof (oriented) ( p − ∂ = 0 and hence d = 0from (B.1). This latter identity is the exterior calculus counterpart of the vector calculusidentities ∇ × ∇ = and ∇ · ∇× = 0.The generalized Stokes’ theorem recovers Stokes’ and Gauss’ theorems of vector cal-culus for p = 1 ,
2, respectively, and the fundamental theorem of calculus for p = 0.27 ppendix C. Discrete Maxwell’s equations By pairing Faraday’s law for the TE φ field set in (26) with K -th 2-cells σ K (2) of theFE grid (primal mesh) and applying the generalized Stokes’ theorem, we obtain (cid:42) σ K (2) , N (cid:88) j =1 E (cid:107) j,m ( t ) (cid:104) d (cid:48)(cid:107) w (1) j (cid:105)(cid:43) = − (cid:42) σ K (2) , ∂∂t N (cid:88) k =1 B ⊥ k,m ( t ) w (2) k (cid:43) , (C.1) (cid:42)(cid:16) ∂σ K (2) (cid:17) (1) , N (cid:88) j =1 E (cid:107) j,m ( t ) w (1) j (cid:43) = − (cid:42) σ K (2) , ∂∂t N (cid:88) k =1 B ⊥ k,m ( t ) w (2) k (cid:43) . (C.2)Using (cid:16) ∂σ K (2) (cid:17) (1) = (cid:80) N j =1 C K,j σ j (1) , where C K,j is the incidence matrix associated to theexterior derivative applied to 1-forms (curl operator on the mesh), see Appendix D), weobtain [38, 43, 73, 74] N (cid:88) j =1 C K,j E (cid:107) j,m ( t ) = − ∂∂t B ⊥ K,m ( t ) , (C.3)for m = − M φ , ..., M φ . The elements of the incidence matrix take values in the set of {− , , } ,Likewise, pairing (27) with J -th 1-cells σ J (1) of the primal mesh gives (cid:42) σ J (1) , N (cid:88) i =1 E ⊥ i,m ( t ) (cid:104) d (cid:48)(cid:107) w (0) i (cid:105)(cid:43) − (cid:42) σ J (1) , | m | N (cid:88) j =1 E (cid:107) j,m ( t ) w (1) j (cid:43) = − (cid:42) σ J (1) , ∂∂t N (cid:88) j =1 B (cid:107) j,m ( t ) w (1) j (cid:43) , (C.4)and applying generalized Stokes’ theorem to the left-hand side of (C.4) yields (cid:42)(cid:16) ∂σ J (1) (cid:17) (0) , N (cid:88) i =1 E ⊥ i,m ( t ) w (0) i (cid:43) − (cid:42) σ J (1) , | m | N (cid:88) j =1 E (cid:107) j,m ( t ) w (1) j (cid:43) = − (cid:42) σ J (1) , ∂∂t N (cid:88) j =1 B (cid:107) j,m ( t ) w (1) j (cid:43) , (C.5)Similarly to before, we can write (cid:16) ∂σ J (1) (cid:17) (0) = (cid:80) N i =1 G J,i σ i (1) , where G J,i is the incidencematrix associated to the exterior derivative applied to 0-forms (gradient operator on themesh), and obtain N (cid:88) i =1 G J,i E ⊥ i,m ( t ) − | m | E (cid:107) J,m ( t ) = − ∂∂t B (cid:107) J,m ( t ) , (C.6)An analogous procedure can be used to obtain the discrete rendering of Ampere’s lawfor on the dual mesh. 28 igure D.17: Example (primal) unstructured mesh. Appendix D. Incidence Matrices
Incidence matrices can be used to represent on a mesh the discrete exterior deriva-tive or, equivalently, the grad, curl, and div operators distilled from their metric struc-ture [38, 43, 72]. Since, from (B.1), the discrete exterior derivative can be seen as thedual of the boundary operator, incidence matrices encode the relationship between eachoriented p -cell of the mesh and its boundary oriented ( p − D curl ], of size N × N , there arethree edges wrapping face number 6: edges 8, 9, and 20. As a result, [ D curl ] , = 1,[ D curl ] , = −
1, and [ D curl ] , = 1. The sign is determined by comparing the intrinsicorientation of each edge with the curl in Fig. D.17: if they are opposite, the elementis −
1, otherwise it is +1. Furthermore, [ D curl ] ,j = 0 for all other j − th edges. This isrepresented in Fig. D.18a, which shows the entire [ D curl ] for this mesh. A curl orien-tation on each face is supposed to follow the intrinsic orientation of the first local edge(i.e. an edge with the smallest index among three edges for the face). Likewise, if weconsider [ D grad ], of size N × N , there are two nodes connected to edge 10: nodes 4 and5. The corresponding elements are [ D grad ] , = − D grad ] , = 1. The element forthe diverging node with the gradient (the intrinsic edge orientation) in Fig. D.17 is − ppendix E. Discrete Hodge Matrix A (discrete) Hodge star operator encodes all metric information and is used to trans-fer information between the primal and dual meshes [38, 40, 47, 51, 75]. Here, we usea Galerkin-Hodge construction [40, 41, 54, 75], which leads to symmetric positive defi-nite matrices and enables energy-conserving discretizations with standard local energypositivity in arbitrary simplicial meshes [43]. As noted in Section 2, the Galerkin-Hodgeoperator is not a natural consequence of DEC [57].The Hodge operator also incorporates the constitutive properties (permittivity andpermeability) of the background medium [30]. Inhomogeneous and anisotropic media canbe easily dealt with by incorporating piecewise constant permittivity and permeabilityover each cell, for example. In the present FETD-BOR solver, the elements of the Hodgematrices including the radial scaling factor from the cylindrical metric are assembled byadding the contributions from all cells as:[ (cid:63) (cid:15) ] → J,j = N (cid:88) k =1 (cid:90) Ω k ( (cid:15) k ρ k ) W (1) J · W (1) j dV, (E.1) (cid:2) (cid:63) µ − (cid:3) → K,k = N (cid:88) k =1 (cid:90) Ω k (cid:0) µ − k ρ k (cid:1) W (2) K · W (2) k dV, (E.2)[ (cid:63) (cid:15) ] → I,i = N (cid:88) k =1 (cid:90) Ω k (cid:0) (cid:15) k ρ − k (cid:1) (cid:104) W (0) I ˆ φ (cid:105) · (cid:104) W (0) i ˆ φ (cid:105) dV, (E.3) (cid:2) (cid:63) µ − (cid:3) → J,j = N (cid:88) k =1 (cid:90) Ω k (cid:0) µ − k ρ − k (cid:1) (cid:104) W (1) J × ˆ φ (cid:105) · (cid:104) W (1) j × ˆ φ (cid:105) dV, (E.4)where Ω k is the area of the k − th cell, and ρ k = (cid:80) i =1 ρ k i / ρ k i is ρ coordinateof i − th node touching k − th face and for simplicity we have assumed isotropic mediaassuming permittivity and permeability values (cid:15) k and µ k , resp., on cell k . Since Whitneyforms have compact support, we can express the global discrete Hodge matrix as a sumof local matrices (excluding element-wise permittivity and permeability information) forthe K -th face as [ T ] → K = ∆ K / /
12 1 / /
12 1 / / /
12 1 /
12 1 / , (E.5)[ T ] → K = ∆ K T → T → T → T → T → T → T → T → T → , (E.6)[ T ] → K = 4∆ K ( ∇ λ × ∇ λ ) · ˆ φ, (E.7)30here ∆ K is the area of K -th face and T → = ∇ λ · ∇ λ ∇ λ · ∇ λ − ∇ λ · ∇ λ , (E.8) T → = ∇ λ · ∇ λ − ∇ λ · ∇ λ − ∇ λ · ∇ λ , (E.9) T → = ∇ λ · ∇ λ − ∇ λ · ∇ λ ∇ λ · ∇ λ , (E.10) T → = T → , (E.11) T → = ∇ λ · ∇ λ ∇ λ · ∇ λ ∇ λ · ∇ λ , (E.12) T → = ∇ λ · ∇ λ ∇ λ · ∇ λ ∇ λ · ∇ λ , (E.13) T → = T → , (E.14) T → = T → , (E.15) T → = ∇ λ · ∇ λ ∇ λ · ∇ λ ∇ λ · ∇ λ . (E.16)Due to the local support of the Whitney forms, the above Hodge matrices are verysparse (and diagonally dominant). Their sparsity patterns for the mesh in Fig. D.17 areprovided in Fig. E.19. The number of non-zero elements per row (or column) in theseHodge matrices is invariant with respect to the mesh size, so the sparsity increases forlarger meshes. Appendix F. Cartesian-like PML implementation
A perfectly matched layer (PML) is used to absorb outgoing waves in FE simulations,enabling analysis of open-domain problems [76, 77]. As described before, in the presentFETD-BOR the spatial discretization is performed in the meridian plane mapped ontoa Cartesian domain with the cylindrical metric factor transferred to the constitutiverelations. The resulting constitutive relations correspond to a medium that is inhomo-geneous and doubly anisotropic. As such, a Cartesian PML implementation extendedto such media can be used. Such formulation exists [78] and is adapted here to theFETD-BOR case as follows.In the 2-D Cartesian plane, the PML can be effected as an analytic continuation onthe spatial variables to complex space [77, 78], given by u → ˜ u = (cid:82) u s u ( u (cid:48) ) du (cid:48) where s u ( u (cid:48) ) is a complex stretching variable and u stands for ρ or z . This transformation canalso be expressed as r (cid:48) (cid:107) → ˜ r (cid:48) (cid:107) = ¯¯ Γ · r (cid:48) (cid:107) , (F.1)where ¯¯ Γ = ˆ ρ ˆ ρ (˜ ρ/ρ ) + ˆ z ˆ z (˜ z/z ). As before, the apostrophe (cid:48) in r (cid:48) (cid:107) denotes the transversecoordinates on the 2-D meridian plane. The modified nabla operator (posterior to theTO-based transformation and hence devoid of the 1 /ρ factor in the φ derivative) followingsuch analytical continuation is given by ∇ (cid:48) → ˜ ∇ (cid:48) = ˆ ρ s ρ ∂∂ρ + ˆ φ ∂∂φ + ˆ z s z ∂∂z , (F.2)31r simply ˜ ∇ (cid:48) = ¯¯ S · ∇ (cid:48) , (F.3)where ¯¯ S = ˆ ρ ˆ ρ (1 /s ρ ) + ˆ φ ˆ φ (1) + ˆ z ˆ z (1 /s z ). Following [78], since s u ( u ) and ∂/∂u (cid:48) commutewhen u (cid:54) = u (cid:48) and ¯¯ S is a diagonal tensor, the following identity holds for any vector a inthe Cartesian-like 2-D meridian plane: ∇ (cid:48) × (cid:16) ¯¯ S − · a (cid:17) = (cid:16) det¯¯ S (cid:17) − ¯¯ S · (cid:16) ¯¯ S · ∇ (cid:48) (cid:17) × a . (F.4)Applying this analytic continuation to (18), (19), (32), and (33) in the Fourier domain(with time convention of e jωt ) yields the modified Maxwell’s equations for each mode m as ˜ ∇ (cid:48)(cid:107) × E (cid:48)(cid:107) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) = − jω B (cid:48)⊥ cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) , (F.5)˜ ∇ (cid:48)(cid:107) × E (cid:48)⊥ cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) = − jω B (cid:48)(cid:107) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) + | m | E (cid:48)(cid:107) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) × ˆ φ, (F.6)˜ ∇ (cid:48)(cid:107) × H (cid:48)(cid:107) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) = jω D (cid:48)⊥ cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) , (F.7)˜ ∇ (cid:48)(cid:107) × H (cid:48)⊥ cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) = jω D (cid:48)(cid:107) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) − | m | H (cid:48)(cid:107) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) × ˆ φ, (F.8)with constitutive relations in analytic-continued complex space as D (cid:48) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) = ¯¯ (cid:15) (cid:48) ( ω ) · E (cid:48) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) , (F.9) B (cid:48) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) = ¯¯ µ (cid:48) ( ω ) · H (cid:48) cm (cid:16) ˜ r (cid:48) (cid:107) (cid:17) , (F.10)where the superscript c denotes non-Maxwellian (complex space) fields and ¯¯ (cid:15) (cid:48) and ¯¯ µ (cid:48) indicates constitutive parameters of the original medium incorporating the radial scalingfactors from the TO mapping. Next, using (F.1) and (F.3), we can revert (F.5) − (F.8)back to a real-valued spatial domain by writing (cid:16) ¯¯ S · ∇ (cid:48)(cid:107) (cid:17) × E (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) = − jω B (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) , (F.11) (cid:16) ¯¯ S · ∇ (cid:48)(cid:107) (cid:17) × E (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) = − jω B (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) − | m | ˆ φ × E (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) , (F.12) (cid:16) ¯¯ S · ∇ (cid:48)(cid:107) (cid:17) × H (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) = jω D (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) , (F.13) (cid:16) ¯¯ S · ∇ (cid:48)(cid:107) (cid:17) × H (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) = jω D (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) + | m | ˆ φ × H (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) . (F.14)32sing the identity (F.4), we can rewrite (F.11) − (F.14) as ∇ (cid:48)(cid:107) × (cid:104) ¯¯ S − · E (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:105) = − jω (cid:20)(cid:16) det¯¯ S (cid:17) − ¯¯ S · B (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:21) , (F.15) ∇ (cid:48)(cid:107) × (cid:104) ¯¯ S − · E (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:105) = − jω (cid:20)(cid:16) det¯¯ S (cid:17) − ¯¯ S · B (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:21) − | m | (cid:20)(cid:16) det¯¯ S (cid:17) − ¯¯ S · (cid:110) ˆ φ × E (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:111)(cid:21) , (F.16) ∇ (cid:48)(cid:107) × (cid:104) ¯¯ S − · H (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:105) = jω (cid:20)(cid:16) det¯¯ S (cid:17) − ¯¯ S · D (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:21) , (F.17) ∇ (cid:48)(cid:107) × (cid:104) ¯¯ S − · H (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:105) = jω (cid:20)(cid:16) det¯¯ S (cid:17) − ¯¯ S · D (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:21) + | m | (cid:20)(cid:16) det¯¯ S (cid:17) − ¯¯ S · (cid:110) ˆ φ × H (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:111)(cid:21) . (F.18)We can further verify the identity below (cid:16) det¯¯ S (cid:17) − ¯¯ S · (cid:110) ˆ φ × E (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:111) = ˆ φ × (cid:104) ¯¯ S − · E (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:105) , (F.19) (cid:16) det¯¯ S (cid:17) − ¯¯ S · (cid:110) ˆ φ × H (cid:48)(cid:107) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:111) = ˆ φ × (cid:104) ¯¯ S − · H (cid:48)⊥ cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17)(cid:105) . (F.20)and introduce a new set of fields defined as E (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) = ¯¯ S − · E (cid:48) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) , (F.21) H (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) = ¯¯ S − · H (cid:48) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) , (F.22) D (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) = (cid:16) det¯¯ S (cid:17) − ¯¯ S · D (cid:48) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) , (F.23) B (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) = (cid:16) det¯¯ S (cid:17) − ¯¯ S · B (cid:48) cm (cid:16) ¯¯ Γ · r (cid:48) (cid:107) (cid:17) , (F.24)so that, by substituting (F.21) − (F.24) back into (F.15) − (F.18), and utilizing the identi-ties (F.19) and (F.20), we finally obtain ∇ (cid:48)(cid:107) × E (cid:48)(cid:107) am (cid:16) r (cid:48) (cid:107) (cid:17) = − jω B (cid:48)⊥ am (cid:16) r (cid:48) (cid:107) (cid:17) , (F.25) ∇ (cid:48)(cid:107) × E (cid:48)⊥ am (cid:16) r (cid:48) (cid:107) (cid:17) = − jω B (cid:48)(cid:107) am (cid:16) r (cid:48) (cid:107) (cid:17) + | m | E (cid:48)(cid:107) am (cid:16) r (cid:48) (cid:107) (cid:17) × ˆ φ, (F.26) ∇ (cid:48)(cid:107) × H (cid:48)(cid:107) am (cid:16) r (cid:48) (cid:107) (cid:17) = jω D (cid:48)⊥ am (cid:16) r (cid:48) (cid:107) (cid:17) , (F.27) ∇ (cid:48)(cid:107) × H (cid:48)⊥ am (cid:16) r (cid:48) (cid:107) (cid:17) = jω D (cid:48)(cid:107) am (cid:16) r (cid:48) (cid:107) (cid:17) − | m | H (cid:48)(cid:107) am (cid:16) r (cid:48) (cid:107) (cid:17) × ˆ φ. (F.28)with D (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) = (cid:20)(cid:16) det¯¯ S (cid:17) − (cid:110) ¯¯ S · ¯¯ (cid:15) (cid:48) ( ω ) · ¯¯ S (cid:111)(cid:21) · E (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) , (F.29) B (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) = (cid:20)(cid:16) det¯¯ S (cid:17) − (cid:110) ¯¯ S · ¯¯ µ (cid:48) ( ω ) · ¯¯ S (cid:111)(cid:21) · H (cid:48) am (cid:16) r (cid:48) (cid:107) (cid:17) . (F.30)33he above expressions show that E (cid:48) am , H (cid:48) am , D (cid:48) am , and B (cid:48) am obey Maxwell’s equations inan equivalent PML medium with constitutive parameters given by¯¯ (cid:15) PML = (cid:20)(cid:16) det¯¯ S (cid:17) − (cid:110) ¯¯ S · ¯¯ (cid:15) (cid:48) ( ω ) · ¯¯ S (cid:111)(cid:21) , (F.31)¯¯ µ PML = (cid:20)(cid:16) det¯¯ S (cid:17) − (cid:110) ¯¯ S · ¯¯ µ (cid:48) ( ω ) · ¯¯ S (cid:111)(cid:21) . (F.32)As an example, consider a background medium with¯¯ (cid:15) ( ω ) = (cid:15) ρ ( ω ) 0 00 (cid:15) φ ( ω ) 00 0 (cid:15) z ( ω ) , (F.33)¯¯ µ ( ω ) = µ ρ ( ω ) 0 00 µ φ ( ω ) 00 0 µ z ( ω ) , (F.34)with (cid:15) ρ ( ω ) = (cid:15) φ ( ω ) = (cid:15) z ( ω ) = (cid:16) σ m jω(cid:15) (cid:17) , corresponding to a lossy, isotropic, homoge-neous medium. After the TO-based mapping, we obtain¯¯ (cid:15) (cid:48) ( ω ) = ¯¯ (cid:15) ( ω ) · ¯¯ R (cid:15) = (cid:15) ρ ( ω ) ρ (cid:15) φ ( ω ) ρ
00 0 (cid:15) z ( ω ) ρ , (F.35)¯¯ µ (cid:48) ( ω ) = ¯¯ µ ( ω ) · ¯¯ R µ = µ ρ ( ω ) ρ µ φ ( ω ) ρ
00 0 µ z ( ω ) ρ , (F.36)As a result, by using (F.31) and (F.32), the elements of the resulting PML constitutivetensor write as: (cid:15) PML ρ ( ω ) = (cid:15) (cid:18) σ m jω(cid:15) (cid:19) (cid:0) jω(cid:15) + σ PML ρ (cid:1) ( jω(cid:15) + σ PML z ) , (F.37) (cid:15) PML φ ( ω ) = (cid:15) (cid:18) σ m jω(cid:15) (cid:19) ( jω(cid:15) ) (cid:0) jω(cid:15) + σ PML ρ (cid:1) ( jω(cid:15) + σ PML z ) , (F.38) (cid:15) PML z ( ω ) = (cid:15) (cid:18) σ m jω(cid:15) (cid:19) (cid:0) jω(cid:15) + σ PML z (cid:1)(cid:0) jω(cid:15) + σ PML ρ (cid:1) , (F.39) µ PML ρ ( ω ) = µ (cid:0) jω(cid:15) + σ PML ρ (cid:1) ( jω(cid:15) + σ PML z ) , (F.40) µ PML φ ( ω ) = µ ( jω(cid:15) ) (cid:0) jω(cid:15) + σ PML ρ (cid:1) ( jω(cid:15) + σ PML z ) , (F.41) µ PML z ( ω ) = µ (cid:0) jω(cid:15) + σ PML z (cid:1)(cid:0) jω(cid:15) + σ PML ρ (cid:1) . (F.42)where σ PML ρ and σ PML z are the artificial PML conductivities along ρ and z respectively.The presence of jω factors in the above Fourier-domain elements produce modifica-tions in the corresponding field equations in the time-domain. These modifications are34mplemented using an auxiliary differential equation (ADE) approach as described in,e.g., [29, 30]. Appendix G. Stability Conditions
To determine the stability conditions, we express the field update in matrix form as¯ w n +1 = ¯¯ G · ¯ w n = (cid:16) ¯¯ I + ¯¯ T (cid:17) · ¯ w n (G.1)with ¯ w n = (cid:2) B ⊥ m (cid:3) n − (cid:104) B (cid:107) m (cid:105) n − (cid:2) E ⊥ m (cid:3) n (cid:104) E (cid:107) m (cid:105) n , ¯ w n +1 = (cid:2) B ⊥ m (cid:3) n + (cid:104) B (cid:107) m (cid:105) n + (cid:2) E ⊥ m (cid:3) n +1 (cid:104) E (cid:107) m (cid:105) n +1 , (G.2)and¯¯ T = ¯¯ N × N , ¯¯ N × N , ¯¯ N × N , − ∆ t [ D curl ]¯¯ N × N , ¯¯ N × N , − ∆ t [ D grad ] , ∆ t | m | ¯¯ I N × N ¯¯ N × N , ∆ t ¯¯ X TM φ , − ∆ t ¯¯ X TM φ · [ D grad ] , ∆ t | m | ¯¯ X TM φ ∆ t ¯¯ X TE φ , − ∆ t | m | ¯¯ A , − ∆ t | m | ¯¯ A · [ D grad ] , − ∆ t ¯¯ X TE φ · [ D curl ] − ∆ t | m | ¯¯ A , (G.3)where ¯¯ X TM φ = (cid:16) [ (cid:63) (cid:15) ] → (cid:17) − · [ D grad ] T · (cid:2) (cid:63) − µ (cid:3) → , (G.4)¯¯ X TE φ = (cid:16) [ (cid:63) (cid:15) ] → (cid:17) − · [ D curl ] T · (cid:2) (cid:63) − µ (cid:3) → , (G.5)¯¯ A = (cid:16) [ (cid:63) (cid:15) ] → (cid:17) − · (cid:2) (cid:63) − µ (cid:3) → . (G.6)A necessary condition for stability is (cid:12)(cid:12) λ ¯¯ G (cid:12)(cid:12) ≤ λ ¯¯ G of ¯¯ G [79].When m = 0, the field update equation becomes decoupled into two independentnumerical integrators for TE φ and TM φ fields. In this case, following [41], we can easilyobtain the stability criteria for both polarizations in closed form as∆ t TE φ ,m =0 ≤ (cid:114) max (cid:16) λ X TE φ · [ D curl ] (cid:17) , (G.7)∆ t TM φ ,m =0 ≤ (cid:114) max (cid:16) λ X TM φ · [ D grad ] (cid:17) , (G.8)where λ X TE φ · [ D curl ] and λ X TM φ · [ D grad ] denote the eigenvalues of X TE φ · [ D curl ] and X TM φ · [ D grad ] respectively. 35hen m (cid:54) = 0, we can simply represent ¯¯ G using 2 × X and [ D ] as¯¯ G = (cid:20) ¯¯ I ( N + N ) × ( N + N ) , − ∆ t [ D ]∆ t ¯¯ X , ¯¯ I ( N + N ) × ( N + N ) − ∆ t ¯¯ X · [ D ] (cid:21) (G.9)where ¯¯ X = (cid:20) ¯¯ N × N , ¯¯ X TM φ ¯¯ X TE φ , − | m | ¯¯ A (cid:21) , (G.10)and [ D ] = (cid:20) ¯¯ N × N , [ D curl ][ D grad ] − | m | ¯¯ I N × N (cid:21) . (G.11)Therefore, the stability condition is similarly obtained as∆ t m (cid:54) =0 ≤ (cid:114) max (cid:16) λ ¯¯ X · [ D ] (cid:17) (G.12)where λ ¯¯ X · [ D ] are the eigenvalues of ¯¯ X · [ D ]. Note that in this case the maximum timestep depends on the modal index magnitude | m | . ReferencesReferences [1] J.-M. Jin, The finite element method in electromagnetics, John Wiley & Sons, New Jersey, 2015.[2] J.-F. Lee, G. M. Wilkins, R. Mitra, Finite-element analysis of axisymmetric cavity resonator usinga hybrid edge element technique, IEEE Trans. Microw. Theory Techn. 41 (11) (1993) 1981–1987.[3] F. L. Teixeira, J. R. Bergmann, Moment-method analysis of circularly symmetric reflectors usingbandlimited basis functions, IEE Proc. - Microw. Antennas Prop. 144 (3) (1997) 179–183.[4] F. L. Teixeira, J. R. Bergmann, B-spline basis functions for moment-method analysis of axisym-metric reflector antennas, Microw. Opt. Tech. Lett. 14 (3) (1997) 188–191.[5] G. M. Wilkins, J. F. Lee, R. Mittra, Numerical modeling of axisymmetric coaxial waveguide dis-continuities, IEEE Trans. Microw. Theory Techn. 39 (8) (1991) 1323–1328.[6] A. D. Greenwood, J.-M. Jin, Finite-element analysis of complex axisymmetric radiating structures,IEEE Trans. Antennas Propag. 47 (8) (1999) 1260–1266.[7] X. Rui, J. Hu, Q. H. Liu, Higher order finite element method for inhomogeneous axisymmetricresonators, Progress In Electromagnetics Research B 21 (2010) 189–201.[8] W. Tierens, D. D. Zutter, BOR-FDTD subgridding based on finite element principles, Journal ofComputational Physics 230 (12) (2011) 4519 – 4535. doi:https://doi.org/10.1016/j.jcp.2011.02.028 .[9] D.-Y. Na, Y. A. Omelchenko, H. Moon, B.-H. V. Borges, F. L. Teixeira, Axisymmetric charge-conservative electromagnetic particle simulation algorithm on unstructured grids: Application tomicrowave vacuum electronic devices, J. Comp. Phys. 346 (2017) 295 – 317.[10] A. Khebir, J. D’Angelo, J. Joseph, A new finite element formulation for RF scattering by complexbodies of revolution, IEEE Transactions on Antennas and Propagation 41 (5) (1993) 534–541. doi:10.1109/8.222272 .[11] L. Medgyesi-Mitschang, J. Putnam, Electromagnetic scattering from axially inhomogeneous bodiesof revolution, IEEE Transactions on Antennas and Propagation 32 (8) (1984) 797–806. doi:10.1109/TAP.1984.1143430 .
12] A. D. Greenwood, J.-M. Jin, A novel efficient algorithm for scattering from a complex BOR usingmixed finite elements and cylindrical PML, IEEE Trans. Antennas Propagat. 47 (4) (1999) 620–629.[13] A. N. O’Donnell, R. J. Burkholder, High-frequency asymptotic solution for the electromagneticscattering from a small groove around a conical or cylindrical surface, IEEE Trans. Antennas Prop.61 (2) (2013) 1003–1008.[14] Y. B. Zhai, X. W. Ping, W. X. Jiang, T. J. Cui, Finite-element analysis of three-dimensionalaxisymmetric invisibility cloaks and other metamaterial devices, Commun. Comput. Phys. 8 (4)(2010) 823–834.[15] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, M. Paszynski, Simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp finite element method, SIAMJ. Appl. Math 66 (6) (2006) 2085–2106.[16] M. S. Novo, L. C. da Silva, F. L. Teixeira, Comparison of coupled-potentials and field-based finite-volume techniques for modeling of borehole EM tools, IEEE Geosci. Remote Sens. Lett. 5 (2) (2008)209–211.[17] M. S. Novo, L. C. da Silva, F. L. Teixeira, Three-dimensional finite-volume analysis of directionalresistivity logging sensors, IEEE Trans. Geosci. Remote Sens. 48 (2) (2010) 1151–1158.[18] D. Hong, W. F. Huang, H. Chen, Q. H. Liu, Novel and stable formulations for the response ofhorizontal-coil eccentric antennas in a cylindrically multilayered medium, IEEE Trans. AntennasPropag. 65 (4) (2017) 1967–1977.[19] S. Yang, D. Hong, W. F. Huang, Q. H. Liu, A stable analytic model for tilted-coil antennas in aconcentrically cylindrical multilayered anisotropic medium, IEEE Geosci. Remote Sens. Lett. 14 (4)(2017) 480–483.[20] Y. Fang, Z. Y. J. Dai, J. Zhou, Q. H. Liu, Through-casing hydraulic fracture evaluation by inductionlogging i: An efficient EM solver for fracture detection, IEEE Trans. Geosci. Remote Sens. 55 (2)(2017) 1179–1188.[21] Y.-K. Hue, F. L. Teixeira, L. S. Martin, M. S. Bittar, Three-dimensional simulation of eccentricLWD tool response in boreholes through dipping formations, IEEE Trans. Geosci. Remote Sens.43 (2) (2005) 257–268.[22] M. F. Wong, M. Prak, V. F. Hanna, Axisymmetric edge-based finite element formulation for bodiesof revolution: Application to dielectric resonators, IEEE MTT-S Digest (1995) 285–288.[23] F. L. Teixeira, Time-domain finite-difference and finite-element methods for Maxwell equationsin complex media, IEEE Trans. Antennas Propag. 56 (2008) 2150–2166. doi:10.1109/TAP.2008.926767 .[24] H. Moon, F. L. Teixeira, Y. A. Omelchenko, Exact charge-conserving scattergather algorithm forparticle-in-cell simulations on unstructured grids: A geometric perspective, Comput. Phys. Com-mun. 194 (2015) 43–53. doi:http://dx.doi.org/10.1016/j.cpc.2015.04.014 .[25] D.-Y. Na, H. Moon, Y. A. Omelchenko, F. L. Teixeira, Local, explicit, and charge-conservingelectromagnetic particle-in-cell algorithm on unstructured grids, IEEE Trans. Plasma Sci. 44 (2016)1353–1362. doi:10.1109/TPS.2016.2582143 .[26] R. A. Chilton, R. Lee, The discrete origin of FETD-newmark late time instability, and a correctionscheme, J. Comput. Phys. 224 (2007) 1293–1306.[27] B. He, F. L. Teixeira, On the degrees of freedom of lattice electrodynamics, Phys. Lett. A 336(2005) 1–7. doi:http://dx.doi.org/10.1016/j.physleta.2005.01.001 .[28] B. He, F. L. Teixeira, Sparse and explicit FETD via approximate inverse Hodge (mass) matrix,IEEE Microw. Wireless Compon. Lett. 16 (2006) 348–350.[29] B. Donderici, F. L. Teixeira, Conformal perfectly matched layer for the mixed finite-element time-domain method, IEEE Trans. Antennas Propag. 56 (4) (2008) 1017–1026.[30] B. Donderici, F. L. Teixeira, Mixed finite-element time-domain method for transient Maxwell equa-tions in doubly dispersive media, IEEE Trans. Microw. Theory Techn. 56 (1) (2008) 113–120. doi:10.1109/TMTT.2007.912217 .[31] F. L. Teixeira, W. C. Chew, Differential forms, metrics, and the reflectionless absorption of elec-tromagnetic waves, J. Electromagn. Waves Appl. 13 (1999) 665–686. doi:http://dx.doi.org/10.1163/156939399X01104 .[32] J. B. Pendry, D. Schurig, D. R. Smith, Controlling electromagnetic fields, Science 312 (2006) 1780–1782. doi:10.1126/science.1125907 .[33] B. He, F. L. Teixeira, Differential forms, Galerkin duality, and sparse inverse approximations infinite element solutions of Maxwell equations, IEEE Trans. Antennas Propag. 55 (2007) 1359–1368. doi:10.1109/TAP.2007.895619 .[34] J. A. Silva-Macedo, M. A. Romero, B.-H. V. Borges, An extended FDTD method for the analysis f electromagnetic field rotations and cloaking devices, Progress In Electromagnetics Research 87(2008) 183–196.[35] O. Ozgun, M. Kuzuoglu, Software metamaterials: Transformation media based multiscale tech-niques for computational electromagnetics, J. Comput. Phys. 236 (2013) 203–219.[36] O. Ozgun, M. Kuzuoglu, Cartesian grid mapper: Transformation media for modeling arbitrarycurved boundaries with Cartesian grids, IEEE Antennas Wireless Propag. Lett. 13 (2014) 1771–1774.[37] L. Kettunen, K. Forsman, A. Bossavit, Discrete spaces for div and curl-free fields, IEEE Transactionson Magnetics 34 (5) (1998) 2551–2554. doi:10.1109/20.717588 .[38] F. L. Teixeira, W. C. Chew, Lattice electromagnetic theory from a topological viewpoint, J. Math.Phys. 40 (1999) 169–187. doi:http://dx.doi.org/10.1063/1.532767 .[39] D. N. Arnold, R. S. Falk, R. Winther, Finite element exterior calculus, homological techniques, andapplications, Acta Numerica 15 (2006) 1–155.[40] J. Kangas, T. Tarhasaari, L. Kettunen, Reading Whitney and finite elements with hindsight, IEEETransactions on Magnetics 43 (4) (2007) 1157–1160. doi:10.1109/TMAG.2007.892276 .[41] J. Kim, F. L. Teixeira, Parallel and explicit finite-element time-domain method for Maxwell’s equa-tions, IEEE Trans. Antennas Propag. 59 (2011) 2350–2356. doi:10.1109/TAP.2011.2143682 .[42] F. L. Teixeira, Differential forms in lattice field theories: An overview, ISRN Math. Phys. 2013(2013) 16. doi:http://dx.doi.org/10.1155/2013/487270 .[43] F. L. Teixeira, Lattice Maxwell’s equations, Prog. Electromagn. Res. 148 (2014) 113–128. doi:10.2528/PIER14062904 .[44] S. C. Chen, W. C. Chew, Numerical electromagnetic frequency domain analysis with discrete ex-terior calculus, J. Comp. Phys. 350 (2017) 668 – 689. doi:https://doi.org/10.1016/j.jcp.2017.08.068 .[45] D. M. Shyroki, Efficient Cartesian-grid-based modeling of rotationally symmetric bodies, IEEETrans. Microw. Theory Techn. 55 (6) (2007) 1132–1138.[46] B. He, F. L. Teixeira, Mixed E-B finite elements for solving 1-D, 2-D, and 3-D time-harmonicMaxwell curl equations, IEEE Microw. Compon. Lett. 17 (5) (2007) 313–315.[47] B. He, F. L. Teixeira, Geometric finite element discretization of Maxwell equations in primal anddual spaces, Phys. Lett. A 349 (2006) 1–14. doi:http://dx.doi.org/10.1016/j.physleta.2005.09.002 .[48] P. R. Kotiuga, Weitzenbock identities and variational formulations in nanophotonics and micromag-netics, IEEE Transactions on Magnetics 43 (4) (2007) 1669–1672. doi:10.1109/TMAG.2007.892497 .[49] S. Rao, D. Wilton, A. Glisson, Electromagnetic scattering by surfaces of arbitrary shape, IEEETransactions on Antennas and Propagation 30 (3) (1982) 409–418. doi:10.1109/TAP.1982.1142818 .[50] K. F. Warnick, Numerical Analysis for Electromagnetic Integral Equations, Artech House, Boston,2008.[51] A. Gillette, C. Bajaj, Dual formulations of mixed finite element methods with applications,Computer-Aided Design 43 (10) (2011) 1213 – 1221, solid and Physical Modeling 2010. doi:https://doi.org/10.1016/j.cad.2011.06.017 .[52] A. Bossavit, Whitney forms: A class of finite elements for three-dimensional computations in elec-tromagnetism, IEE Proc., Part A: Phys. Sci., Meas. Instrum., Manage. Educ. 135 (1988) 493–500. doi:10.1049/ip-a-1.1988.0077 .[53] P. W. Gross, P. R. Kotiuga, Electromagnetic Theory and Computation: A Topological Approach,Cambridge Univ. Press, Cambridge, 2004.[54] J. Dodziuk, Finite-difference approach to the Hodge theory of harmonic forms, Am. J. Math. 98 (1)(1976) 79–104.URL [55] T. Tarhasaari, L. Kettunen, A. Bossavit, Some realizations of a discrete Hodge operator: a rein-terpretation of finite element techniques [for EM field analysis], IEEE Trans. Magn. 35 (3) (1999)1494–1497. doi:10.1109/20.767250 .[56] A. Bossavit, Computational electromagnetism and geometry (5): The Galerkin hodge, J. JapanSoc. Appl. Electromagn. Mech. 8 (2) (2000) 203–209.[57] P. R. Kotiuga, Theoretical limitations of discrete exterior calculus in the context of computationalelectromagnetics, IEEE Trans. Magn. 44 (6) (2008) 1162–1165. doi:10.1109/TMAG.2007.915998 .[58] J. E. Lebaric, D. Kajfez, Analysis of dielectric resonator cavities using the finite integration tech-nique, IEEE Trans. Microw. Theory Techn. 37 (11) (1989) 1740–1748.[59] H. Li, H. Wang, Investigation of eccentricity effects and depth of investigation of azimuthal resis-tivity LWD tools using 3d finite difference method, J. Petroleum Sci. Eng. 143 (2016) 211–225.
60] Z. Q. Zhang, Q. H. Liu, Simulation of induction-logging response using conjugate gradient methodwith nonuniform fast Fourier and fast Hankel transforms, Radio Sci. 36 (4) (2001) 599–608.[61] M. S. Novo, L. C. da Silva, F. L. Teixeira, A comparative analysis of Krylov solvers for three-dimensional simulations of borehole sensors, IEEE Geosci. Remote Sens. Lett. 8 (1) (2011) 98–102.[62] H. O. Lee, F. L. Teixeira, L. E. S. Martin, M. S. Bittar, Numerical modeling of eccentered LWDborehole sensors in dipping and fully anisotropic earth formations, IEEE Trans. Geosci. RemoteSens. 50 (3) (2012) 727–735.[63] G. S. Liu, F. L. Teixeira, G. J. Zhang, Analysis of directional logging tools in anisotropic andmultieccentric cylindrically-layered earth formations, IEEE Trans. Antennas Propag. 60 (1) (2012)318–327.[64] G. S. Rosa, J. R. Bergmann, F. L. Teixeira, A robust mode-matching algorithm for the analysis oftriaxial well-logging tools in anisotropic geophysical formations, IEEE Trans. Geosci. Remote Sens.55 (5) (2017) 2534–2545.[65] S. H. Gold, G. S. Nusinovich, Review of high-power microwave source research, Rev. Sci. Instrum.68 (1997) 3945–3974. doi:http://dx.doi.org/10.1063/1.1148382 .[66] J. M. Dawson, Particle simulation of plasmas, Rev. Mod. Phys. 55 (1983) 403–447. doi:10.1103/RevModPhys.55.403 .[67] R. W. Hockney, J. W. Eastwood, Computer Simulation Using Particles, CRC Press, New York,1988.[68] C. K. Birdsall, A. B. Langdon, Plasma Physics via Computer Simulation, CRC Press, New York,2004.[69] R. A. Cairns, A. D. R. Phelps, Generation and Application of High Power Microwaves, CRC Press,New York, 1997.[70] U. Chipengo, M. Zuboraj, N. K. Nahar, J. L. Volakis, A novel slow-wave structure for high-power-band backward wave oscillators with mode control, IEEE Trans. Plasma Sci. 43 (2015) 1879–1886. doi:10.1109/TPS.2015.2431647 .[71] S. S. Cairns, The generalized theorem of Stokes, Trans. Amer. Math. Soc. 40 (1936) 167–174.[72] L. Kettunen, K. Forsman, A. Bossavit, Gauging in Whitney spaces, IEEE Transactions on Magnetics35 (3) (1999) 1466–1469. doi:10.1109/20.767243 .[73] T. J. Hughes, W. K. Liu, T. K. Zimmermann, Lagrangian-Eulerian finite element formulation forincompressible viscous flows, Comput. Method Appl. M. 29 (1981) 329–349.[74] A. H. Guth, Existence proof of a nonconfining phase in four-dimensional U(1) lattice gauge theory,Phys. Rev. D 21 (1980) 2291–2307.[75] T. Tarhasaari, L. Kettunen, A. Bossavit, Some realizations of a discrete hodge operator: a rein-terpretation of finite element techniques [for EM field analysis], IEEE Transactions on Magnetics35 (3) (1999) 1494–1497. doi:10.1109/20.767250 .[76] J.-P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal ofComputational Physics 114 (2) (1994) 185 – 200. doi:https://doi.org/10.1006/jcph.1994.1159 .[77] F. L. Teixeira, W. C. Chew, Complex space approach to perfectly matched layers: a review and somenew developments, International Journal of Numerical Modelling: Electronic Networks, Devices andFields 13 (5) (2000) 441–455. doi:10.1002/1099-1204(200009/10)13:5<441::AID-JNM376>3.0.CO;2-J .[78] F. L. Teixeira, W. C. Chew, General closed-form PML constitutive tensors to match arbitrarybianisotropic and dispersive linear media, IEEE Microw. Guided Wave Lett. 8 (6) (1998) 223–225.[79] S. Wang, F. L. Teixeira, Some remarks on the stability of time-domain electromagnetic simulations,IEEE Trans. Antennas Propag. 52 (3) (2004) 895–898. a)(b)Figure D.18: Incidence matrices for (a) curl [ D curl ] and (b) gradient (cid:2) D grad (cid:3) operators for the mesh inFig. D.17. a) (b)(c) (d)Figure E.19: Sparsity patterns for discrete Hodge matrices corresponding to the toy mesh depicted inFig. D.17: (a) [ (cid:63) (cid:15) ] → , (b) [ (cid:63) (cid:15) ] → , (c) (cid:104) (cid:63) − µ (cid:105) → , and (d) (cid:104) (cid:63) µ − (cid:105) → ..