# The Time-Harmonic Discontinuous Galerkin Method as a Robust Forward Solver for Microwave Imaging Applications.

1. INTRODUCTIONIn the field of computational electromagnetics, the discontinuous Galerkin method (DGM) has been applied to time-domain (TD) simulations of Maxwell's equations [1, 2] and is usually formulated as a high-order version of the time-domain finite-volume method (FVM) [3]. The benefit of TD-DGMs over time-domain finite-element methods (FEM) is that DGM implementations result in a global mass matrix that is locally invertible, thereby facilitating the implementation of explicit time-stepping schemes [1,4,5]. More recently, DGM formulations have been applied to elliptic partial differential equations (PDEs) [6-14], including the time-harmonic Maxwell's equations. In the frequency domain, both DGM and FEM formulations result in large sparse systems of equations that must be solved to determine the electromagnetic fields of interest. In this context, DGM formulations have several attractive properties including the ability to support unstructured and even non-conformal meshes, as well as local order refinement where the approximation orders in neighbouring mesh elements may differ [1,9]. The DGM is usually applied directly to Maxwell's first-order curl equations, providing simultaneous simulation of electric and magnetic fields arising from both electric and magnetic media and sources [1,9]. Solving for both electric and magnetic fields requires additional computational time and memory compared to standard FEM implementations that are formulated for the vector wave equation [15].

While the number of publications targeting general DGM formulations for time- harmonic applications has increased in recent years, these works are largely focused on theoretical aspects of interior cavity problems [6,7], DGM formulations for Maxwell's curl equations [8,9,16] and time-harmonic fluid dynamics [10,11]. Existing work focuses on mathematical issues such as convergence [12-14]. To the best of our knowledge, no existing work presents DGM formulations for the time-harmonic electromagnetic vector wave equations with both electric and magnetic sources and locally varying constitutive parameters. Our desire to model both dielectric and magnetic materials and sources is motivated by biomedical microwave imaging (MWI) applications. PDE formulations of Maxwell's equations readily support inhomogeneous background materials and boundary conditions, and have been recently exploited in MWI algorithms by using an FEM forward solver [17,18]. Recent research suggests the use of magnetic contrast agents for cancer-detecting MWI [19] and so we consider forward solvers that support both magnetic and dielectric materials. Also, in MWI a scattered-field formulation is typically adopted implying that the presence of magnetic materials forces forward solvers to support magnetic, as well as electric, contrast sources. These considerations lead us naturally to DGM formulations.

Herein, we provide the details of the nodal DGM formulations for solving the time-harmonic vector wave equations when both electric and magnetic sources are present. In Section 2, the governing equations are provided. Section 3 reviews the standard nodal unstructured DGM (NUDG) applied directly to Maxwel's curl equations and extends the concepts to the vector wave equations for both electric and magnetic fields. Attention is given to how the electric and magnetic sources contribute to the numerical fluxes required for the wave equations. In addition, we formulate exact radiating boundary conditions (ERBCs) for the DGM. Section 4 presents numerical results where ERBCs are shown to improve simulation accuracy for DGM solutions when compared to the frequently used Silver-Muller absorbing boundary conditions [1,9], while the wave equation formulations reduce the computational cost. The numerical experiments are limited to two-dimensional problems applicable to tomographic MWI systems and demonstrate that the DGM formulations for both the curl equations and the wave equations are accurate for conducting, dielectric and magnetic scatterers. Finally, the mechanism of approximating complex media using the high-order expansions provided by the DGM basis is explored.

2. THE CONTINUOUS EQUATIONS

We consider Maxwell's curl equations and a time-harmonic scattered-field formulation in a domain [OMEGA] bounded by [GAMMA]:

[mathematical expression not reproducible] (1)

where [??] [member of] [OMEGA] is the Cartesian position vector, [omega] is the angular frequency, the medium is characterized by the complex scalar permittivity [epsilon]([??]) and permeability [mu]([??]), j = [square root of -1] and an exp (j[omega]t) time-dependence is assumed and suppressed. The electric source [mathematical expression not reproducible] and magnetic source [mathematical expression not reproducible] support the scattered electric and magnetic fields, denoted respectively by [[??].sup.sct] and [[??].sup.sct]. For an inhomogeneous background characterized by [[epsilon].sub.b]([??]) and [[mu].sub.b]([??]) the sources are

[mathematical expression not reproducible], (2)

where the incident fields [[??].sup.inc] and [[??].sup.inc] are supported by sources [mathematical expression not reproducible] and [mathematical expression not reproducible] according to

[mathematical expression not reproducible]. (3)

Both the incident and scattered fields are subject to appropriate boundary conditions on [GAMMA]. We assume that [GAMMA] is made up of perfect electric conductors (PECs), and surfaces where absorbing boundary conditions (ABCs) are applied, i.e., [GAMMA] = [[GAMMA].sup.PEC] [union] [[GAMMA].sup.ABC]. The DGM is most commonly applied to systems such as (1) as a direct discretization of Maxwell's curl equations. However, it may also be formulated for the vector wave equations for the electric and magnetic fields [20]:

[mathematical expression not reproducible] (4)

Note that the magnetic wave equation is obtained from the electric wave equation by duality: [mathematical expression not reproducible].

3. THE DISCONTINUOUS GALERKIN METHOD

To formulate the DGM for the curl equations (1) and the vector wave equations (4) we assume a partition [[union].sub.n][V.sub.n] [approximately equal to] [OMEGA], where [V.sub.n] denotes the nth finite-volume (element). On [V.sub.n] each field component is expanded as a polynomial of order p = p(n). The expansion order can vary with element index n. Choosing the basis as interpolating Lagrange polynomials corresponding to a set of M = M(p(n)) interpolating nodes results in NUDG methods [1]. Thus, for example, the z-component of the scattered electric field is expanded in [V.sub.n] as

[E.sub.z]([??]) = [M.summation over (m=1)] [E.sub.z]([[??].sup.n.sub.m]) [l.sup.n.sub.m] ([??]) + O([[absolute value of [V.sub.n].sup.p(n)+1]), [??] [member of] [V.sub.n]. (5)

Here, [l.sup.n.sub.m]([??]) is the mth Lagrange polynomial of order p(n) in [V.sub.n], {[[??].sup.n.sub.1], [[??].sup.n.sub.2], ... [[??].sup.n.sub.M]} is the set of nodal points and [absolute value of [V.sub.n]] denotes the element measure. The z- component of the electric field is uniquely defined in [V.sub.n] to order p(n) by the nodal coefficient vector [[E.bar].sup.n.sub.z] [member of] [C.sup.Mx1], where [([[E.bar].sup.n.sub.z]).sub.m] = [E.sub.z]([[??].sup.n.sub.m]) such that

[E.sub.z]([??]) = [([l.sup.n]([??])).sup.T] [[E.bar].sup.n.sub.z], [??] [member of] [V.sub.n], [([l.sup.n]([??])).sup.T] = [[l.sup.n.sub.1]([??]), [l.sup.n.sub.2]([??]), ..., [l.sup.n.sub.M]([??])] (6)

where T denotes transposition. Details on implementing this basis can be found in [1].

3.1. DGM Testing

Two approaches to producing DGM systems from Maxwell's curl equations are the 'weak' and 'strong' forms [1]. We consider the strong form, which for Maxwell's curl equations can be obtained by testing (1) with a scalar function [psi]([??]) and applying integration by parts twice over a finite volume [V.sub.n] [1]:

[mathematical expression not reproducible]. (7)

We have suppressed the spatial dependence for notational brevity. In these equations, [??] = [??](x) denotes the outward unit normal to the boundary [partial derivative][V.sub.n] of [V.sub.n]. We emphasize that high-order representations of [epsilon] = [epsilon]([??]) and [mu] = [mu]([??]) are possible as will be implemented in the sequel, thus the constitutive parameters cannot be pulled out of the integral. The tested equations (7) are local to [V.sub.n] with the exception of the boundary integral. The boundary values are referred to as (numerical) fluxes and the superscript [conjunction] is used to indicate that they should incorporate information from both sides of the boundary in a way that is consistent with the physical boundary conditions dictated by Maxwell's equations ([dagger]) [1]. Testing the electric vector wave equation in (4) with [psi]([??]) and applying integration by parts twice results in the strong form:

[mathematical expression not reproducible] (8)

where the equivalent tested magnetic vector wave equation can be obtained by duality.

3.2. The Local DGM Discrete System for Maxwell's Curl Equations

Substituting expansions of the form (5) for each field component and source term into the tested Maxwell's curl equations (7) and letting the test function [psi] range over all functions in the expansion set {[l.sup.n.sub.m]([??])}, (i.e., applying a Galerkin testing procedure), results in 6M equations in 6M unknowns on [V.sub.n]. For example the x-components of the curl equations become

[mathematical expression not reproducible]. (9)

Here, [M.sup.n] and [S.sup.n.sub.[zeta]] [member of] [R.sup.MxM] are the element mass and stiffness matrices, [[??].sup.E/H,n.sub.[zeta]] are the tested numerical fluxes, and [J.sup.sct,n.sub.[zeta]] and [M.sp.sct,n.sub.[zeta]] are the nodal values of the sources in [V.sub.n]. The M x M diagonal matrices [[epsilon].sup.n] and [[mu].sup.n] respectively contain [[[epsilon].bar].sup.n] and [[[mu].bar].sup.n] on their diagonals and permit locally varying constitutive parameters, the interpretation being that for the first term in (9) we are expanding [epsilon]([??])[E.sub.x]([??]) in the basis. The elemental matrices are computed as

[mathematical expression not reproducible]. (10)

Details for evaluating the elemental matrices are provided in [1]. The tested numerical fluxes couple the fields to adjacent elements and will be described in Section 3.4. The remaining equations can be obtained by cyclic permutations x [right arrow] y [right arrow] z [right arrow] x. Combining the six single-component equations into a single system for [V.sub.n] (suppressing the dependence on n) we arrive at the local system

[mathematical expression not reproducible] (11)

where [alpha] = [[alpha].sup.n] = j[omega][M.sup.n][[epsilon].sup.n], and where [beta] = [[beta].sup.n] = j[omega][M.sup.n][[mu].sup.n]. Collecting the equations for each element into a global system results in the discrete DGM formulation of Maxwell's curl equations.

From the three-dimensional local systems, we can obtain the two-dimensional transverse magnetic (TM) and transverse electric (TE) strong-form DGM systems. For example, by retaining only the rows (and corresponding columns) relating to [E.sub.z], [H.sub.x] and [H.sub.y] from the system (11) we obtain the z- invariant local TM system

[mathematical expression not reproducible] (12)

where we must interpret the entries of the mass and stiffness matrices as integrals over two-dimensional elements of two-dimensional Lagrange interpolating polynomials, a basis whose implementation is also discussed in detail in [1].

3.3. The Local DGM Discrete System for the Vector Wave Equations

To derive the local DGM system for the electric vector wave equation, we consider the x-component of the electric field equation given in (8). A Galerkin testing procedure on [V.sub.n] gives the local system

[mathematical expression not reproducible], (13)

where a derivation of the required operators is provided in the appendix. The matrices [D.sup.n.sub.[zeta]] are derivative operators such that the product [D.sup.n.sub.[zeta]][[E.bar].sup.n.sub.[zeta]] results in the coefficients of the derivative of [E.sub.[xi]]([??]) with respect to [zeta] in [V.sub.n]. The definition of the tested numerical flux [[??].sup.E,n.sub.x] is specific to this system of equations and should not be confused with the analogous term in (9); its computation will be discussed in Section 3.5. Letting [[gamma].sub.[zeta][xi]] = [(j[omega]).sup.2] [M.sup.n][e.sup.n] - ([S.sup.n.sub.[zeta]][([[mu].sup.n]).sup.- 1][D.sup.n.sub.[zeta]] + [S.sup.n.sub.[xi]][([[mu].sup.n]).sup.- 1][D.sup.n.sub.[xi]]) and [[delta].sub.[zeta][xi]] = [S.sup.n.sub.[zeta]][([[mu].sup.n]).sup.-1][D.sup.n.sub.[xi]], the local system on [V.sub.n] for the electric wave equation (with n suppressed) is

[mathematical expression not reproducible]. (14)

The corresponding system for the magnetic fields is

[mathematical expression not reproducible] (15)

where, if [mathematical expression not reproducible] and similar notation is applied to [[??].sub.[zeta][xi]]. The z-invariant Tm equations can be obtained as two decoupled systems, one from the electric field wave equation:

[[gamma].sub.yx][[E.bar].sup.sct.sub.z] + [[??].sup.E.sub.z] = -j[omega]M[[J.bar].sup.sct.sub.z] - (- [S.sub.y][[mu].sup.-1] [[M.bar].sup.sct.sub.x] + [S.sub.x][[mu].sup.-1] [[bar.M].sup.sct.sub.y]). (16)

and the other from the magnetic field wave equation:

[mathematical expression not reproducible] (17)

We emphasize that the presented theory permits the construction of the DGM discrete wave equation systems from operators readily available in the DGM discretization of Maxwell's curl equations.

3.4. Numerical Fluxes for the Curl Equations

Evaluation of the fluxes for the DGM formulation of Maxwell's curl equations can be found in [1]. Nevertheless, we repeat the equations herein to facilitate the exposition of evaluating the fluxes for the wave equations in the next section. The flux terms in Equation (11) come from the boundary integral in (7). To simplify the notation we introduce

[mathematical expression not reproducible] (18)

so we may write

[mathematical expression not reproducible] (19)

For the moment we consider only the evaluation of [[??].sup.E/H,n] and will return to the full evaluation of the above expression in Section 3.9. From [1], [[??].sup.E/H,n] satisfies

[mathematical expression not reproducible] (20)

where Z = [square root of [mu]/[epsilon]] = [Y.sup.-1] is the impedance of the medium, [[[[??].sup.sct]]] = [[??].sup.sct,-] - [[??].sup.sct,+], {{Z}} = [Z.sup.-] + [Z.sup.+], and negative superscripts denote quantities associated with [V.sub.n] while positive superscripts denote quantities associated with the neighbouring element at a given position on [partial derivative][V.sub.n] ([double dagger]). Thus [[??].sup.sct,-] [equivalent to] [[??].sup.sct,n].

The parameter [alpha] [member of] [0,1] is chosen as zero for a centered flux or 1 for an upwind flux. It follows from duality that for [??] [member of] [partial derivative][V.sub.n] [3, 21],

[mathematical expression not reproducible] (21)

where for notational brevity we have introduced the operators:

[mathematical expression not reproducible]. (22)

with N = N([??]) [member of] [R.sup.3x3] being the [??]x operator:

[mathematical expression not reproducible]. (23)

For TM and TE formulations care must be taken when implementing the operator N. In fact, two different [??]x operators are required due to the different number of components in the electric and magnetic fields. For z-invariant problems we define the following operators:

[mathematical expression not reproducible], (24)

such that, in the TM case, [N.sup.2] [??] becomes [N.sup.H.sub.TM] [N.sup.E.sub.TM] [E.sub.z] while [N.sup.2] [??] becomes [N.sup.E.sub.TM] [N.sup.H.sub.TM] [[H.sub.x] [H.sub.y]].sup.T]. In this way the relationship (21) can be computed for TM problems by constructing [A.sup.-/+] using the appropriate [N.sup.E/H.sub.TM] operators and limiting the field components to [E.sub.z], [H.sub.x] and [H.sub.y].

3.5. Numerical Fluxes for the Vector Wave Equations

For the discrete electric-field wave Equation (14) we define

[mathematical expression not reproducible] (25)

such that the required tested fluxes can be written as:

[mathematical expression not reproducible]. (26)

We seek to cast this in a form that is recognizable from the standard flux calculation (21). Again we only consider evaluating [[??].sup.E,n], while details on how to evaluate (26) are deferred to Section 3.9. As [mathematical expression not reproducible] it follows from (25) and (18) that

[mathematical expression not reproducible] (27)

where we have assumed that the boundary integrals of the terms involving the magnetic sources cancel ([section][section]). The right-hand side of this expression has previously been provided in (21), but some additional work is required because we do not carry the magnetic field in the electric wave equation. To eliminate the magnetic field from the expression we again use the fact that [mathematical expression not reproducible], which when combined with (27) and (21) produces

[mathematical expression not reproducible]. (28)

The curl of the electric field can be computed in a straightforward fashion using standard basis operators found both in [1] and the appendix. The terms that involve the magnetic sources [[??].sup.sct,+/-] will contribute, after testing, to the right-hand side of the global system of equations. The fluxes required for the magnetic vector wave equation can be derived in a similar manner. Letting

[mathematical expression not reproducible] (29)

the required tested fluxes are:

[mathematical expression not reproducible]. (30)

It follows that [mathematical expression not reproducible] resulting in

[mathematical expression not reproducible]. (31)

The relations given in (21), (28) and (31) are valid for interfaces between elements involving finite constitutive parameters. Boundary-specific fluxes can be computed from these by making assumptions regarding the neighbouring element as shown in the following sections.

3.6. Silver-Muller Boundary Conditions

For an element [V.sub.n] adjacent to an ABC, Silver-Muller (SM) boundary condition can be implemented by assuming that [Z.sup.+] = [Z.sup.-] and that [mathematical expression not reproducible] [1, 3]. Using the notation:

[mathematical expression not reproducible] (32)

the Silver-Muller absorbing boundary flux for the curl equations is simply:

[mathematical expression not reproducible] (33)

while for the electric and magnetic vector wave equations respectively we have:

[mathematical expression not reproducible] (34)

3.7. PEC Boundary Conditions

At a point on [partial derivative][V.sub.n] lying on a PEC surface, the numerical flux can be obtained for the total fields by assuming no contribution from across the boundary and taking the limit of the flux terms as [Z.sup.+] tends to 0. The scattered-field fluxes are then recovered by decomposing the total field into scattered and incident fields [3, 21]. For the curl equations we have:

[mathematical expression not reproducible] (35)

such that the curl equation flux is:

[mathematical expression not reproducible]. (36)

Here, the term involving the incident field must be shifted to the right-hand- side of the resulting global system of equations. Due to the definition of the operator [A.sup.-.sub.PEC] the magnetic fields can be taken as zero. With this in mind, the PEC fluxes for the electric and magnetic wave equations are respectively:

[mathematical expression not reproducible]. (37)

3.8. Exact Radiating Boundary Conditions

Silver-Muller (SM) boundary conditions are attractive for DGMs because they are simple to implement and offer relatively accurate results provided that the boundary condition is located sufficiently far from all scattering targets. Perfectly matched layers (PML) offer an alternative to the SM conditions that are generally more accurate but require tuning the spatial profile of an artificially lossy medium [4,22,23]. As an alternative, we use exact radiating boundary conditions (ERBCs). ERBCs use a Huygens' surface contained inside the computational domain to evaluate the fields at the computational boundary [[GAMMA].sup.ABC], and have been previously used for FEMs [24], TD-FVMs [25] and TD finite- difference methods [26, 27]. Our ERBC implementation follows the FVTD approach in [25]. On [[GAMMA].sup.ABC] the quantities [[??].sup.sct,+] and [[??].sup.sct,+] are computed using the Huygens' surface [[GAMMA].sup.ERBC] shown in Fig. 1, which is assumed to enclose all scattered-field sources. The required time-harmonic Huygens' integrals are available, for example, in [15]. The electric field on [[GAMMA].sup.ABC] can be computed from

[mathematical expression not reproducible] (38)

where [[GAMMA].sup.ERBC] is the integration surface. The magnetic field can be obtained by duality. Here [G.sub.b]([??], [??]') is the Green's function for the background medium external to [[GAMMA].sup.ERBC]. For the wave equations additional transformations are used to determine both the electric and magnetic fields. Evaluation of the Huygens' integrals is performed by interpolating the fields to a set of quadrature points on [[GAMMA].sup.ERBC] and evaluating the integrals to a desired order. We assume that [[GAMMA].sup.ABC] and [[GAMMA].sup.ERBC] are sufficiently separated that no singularity extraction is required. To account for variation in the Green's function we typically assume an integration rule whose order is twice that used in the DGM. In order to simplify the interpolation process we assume that the integration surface coincides with element boundaries. When implementing ERBCs for FVTD it was found that evaluating the integration surface values using elements inside the volume it encloses provided stable solutions [25], and we maintain that approach herein.

3.9. Implementation Details

To facilitate efficient construction of the global DGM system of equations, we provide implementation details related to evaluating the tested numerical fluxes (19), (26) and (30). We exploit the particular choice of nodal expansions presented in [1] for simplex elements, such that each simplex face maintains a subset of the nodal points that is commensurate with the expansion order p(n) for (D - 1)-dimensional functions along [partial derivative][V.sub.n]. Examples are shown (for 2D- geometries) in Fig. 2 for orders 1 and 2. Three nodes are required to capture a first-order function in D = 2 dimensions. On each of the three faces there are two points, sufficient to capture first-order behaviour in D = 1 dimensions along the face. Details for generating nodal points maintaining this property are provided in [1].

The tested fluxes for the curl equations (19) are decomposed into a sum of integrals along th D + 1 faces of a D-dimensional simplex [V.sub.n]. We use [partial derivative][V.sub.n,I] to denote the ith face of [V.sub.n] and [[PHI].sup.E,n,i.sub.[zeta]]([??]) to represent the restriction of [[PHI].sup.E,n.sub.[zeta]]([??]) to [??] [member of] [partial derivative][V.sub.n,i]. Along the face the nodal points may be expressed as [mathematical expression not reproducible], where [j.sub.i](k) maps the kth face point index to the nodal index in [V.sub.n]. Expanding each [[PHI].sup.E/H,n,i.sub.[zeta]]([??]) using a nodal basis along [partial derivative][V.sub.n,I] results in

[mathematical expression not reproducible] (39)

where [mathematical expression not reproducible] represents the nodal basis functions corresponding to points along the ith face of [V.sub.n] and where [[[PHI].bar].sup.E/H,n,i.sub.[zeta]] represents the expansion coefficients for the flux functions along [partial derivative][V.sub.n,i]. Computation of the element-face mass matrix, [M.sup.n,i], can be found in [1].

The evaluation of the curl-equation fluxes in (39), using (21) at a given point along [partial derivative][V.sub.n,i], is facilitated when the field components at that point are collected together. To this end assume that [V.sub.n] and [V.sub.n'] share a common edge [partial derivative][V.sub.n,i] and let [[E.bar].sup.sct,n',n.sub.[zeta]] denote the value of [E.sup.sct.sub.[zeta]]([??]) in [V.sub.n'] evaluated at the nodal points in [V.sub.n] that lie on the shared edge [partial derivative][V.sub.n,i], i.e.,

[mathematical expression not reproducible]

and let [mathematical expression not reproducible] denote the collection of all field components in [V.sub.n'] at a point along the edge, i.e.,

[mathematical expression not reproducible]

We may now introduce the compact notation:

[mathematical expression not reproducible] (40)

where concat is used to represent the concatenation of column vectors and the bold face type is used to emphasize that the above quantities are collections of nodal coefficients for multiple functions. Thus,

* [[U.bar].sup.sct,n] contains the coefficients for all field components on [V.sub.n] in the usual volumetric order expected by the elemental system (11).

* [[U.bar].sup.sct,n',n] contains the evaluation of the fields in [V.sub.n'] at the nodal points in [V.sub.n] that lie on [partial derivative][V.sub.n,i] in the usual volumetric order.

* [[??].sup.sct,n',n] is a re-ordering of [[U.bar].sup.sct,n',n] such that all field components at a single point are collected together.

It follows that

[[??].sup.sct,n',n] = P[[U.bar].sup.sct,n'n] = PV(n', n) [[U.bar].sup.sct,n]

where V(n', n) is a block interpolation operator that maps each field component expansion in [V.sub.n'] to the nodes in [V.sub.n] on the shared face, and P is a permutation matrix that then reorders the field components on [partial derivative][V.sub.n,i] to be compatible with numerical flux evaluations. An illustration is provided in Fig. 2. This expression can be evaluated when n' = n or when n' denotes the index of the neighbour to [V.sub.n] across [partial derivative][V.sub.n,i]. When n' = n or p(n) = p(n'), V(n',n) is simply a selection matrix. Otherwise, V(n',n) interpolates the order p(n') expansion of each field component in [V.sub.n'] to the [V.sub.n] nodes on [partial derivative][V.sub.n,i]. As a result, the coefficients of each component of the flux terms in (39) can be computed as:

[mathematical expression not reproducible] (41)

where [P.sup.-1] = [P.sup.T] puts the coefficients back in the usual volumetric order and where [A.sup.-/+] is a block- diagonal matrix whose blocks are composed of [A.sup.-/+] evaluated at the nodes in [V.sub.n] on [partial derivative][V.sub.n,i]. The tested numerical fluxes in (19) can now be evaluated directly, but we note that an additional interpolation is required to compute the constitutive parameters at the nodal points in [V.sub.n] on [partial derivative][V.sub.n,i] when the permittivity and/or permeability vary over [V.sub.n'].

To evaluate the numerical fluxes required for the electric wave Equation (26) we use

[mathematical expression not reproducible] (42)

where [P.sup.-1.sub.E] is a permutation matrix commensurate with the number of field components for the electric field. To evaluate the above fluxes we must compute the magnetic fields in [[U.bar].sup.sct,n] from the electric fields and magnetic sources. For the magnetic wave equation we have

[mathematical expression not reproducible] (43)

where [P.sup.-1.sub.H] is a permutation matrix commensurate with the number of field components for the magnetic field and where the electric field components of [[U.bar].sup.sct,n] must be determined from the magnetic fields and electric sources.

3.10. The Global DGM System

The global system for the curl-equations can be constructed by collecting together the elemental systems (11) resulting in:

K[[U.bar].sup.sct] = R[[U.bar].sup.inc] (44)

where K = K([omega], [epsilon], [[mu].bar]) contains the volumetric and flux matrices that multiply the scattered-field coefficients, and [[U.bar].sup.sct/inc] is the concatenation of [[U.bar].sup.sct/inc,n] over each mesh element. The matrix R = R([omega], [[epsilon].bar], [[mu].bar], [[[epsilon].bar].sub.b], [[[mu].bar].sub.b]) produces the right-hand-side volumetric and flux terms due to the incident fields and the difference between the constitutive parameters and the background, where e and /x denote the collection of complex constitutive coefficients over the entire mesh.

The discrete versions of the electric and magnetic wave equations can be written using similar notation:

[K.sub.E][[U.bar].sup.E,sct] = [R.sub.E][[U.bar].sup.inc], [K.sub.H][[U.bar].sup.H,sct] = [R.sub.H][[U.bar].sup.inc] (45)

where [[U.bar].sup.E/H,sct] contains the electric or magnetic field coefficients collected over all elements and where [K.sub.E/H] and [R.sub.E/H] are the appropriate operators for the electric or magnetic wave equation formulations. Note that in the case of general constitutive parameters, both the electric and magnetic incident fields contribute to the right-hand-side of all three systems.

4. NUMERICAL RESULTS

In this section we present simulation results to demonstrate the performance of the DGM applied to both Maxwell's curl equations and the vector wave equations for transverse magnetic field problems. We denote the Maxwell's curl equations formulation as DGM-TM-ME (or simply ME) and the wave equation formulations as DGM-TM-WE and DGM-TM-WH (simply WE or WH), for the electric and magnetic wave equations respectively. In all tests, the upwind parameter a = 1 was used. Due to the modest sizes of the problems considered, Matlab was used to implement the DGM formulations, using basis function and plotting routines provided in [1]. Simulations were executed on an iMac equipped with a single 2-core 3.06 GHz Intel Core i3 processor and 8 GB of 1333 MHz DDR3 memory. For all examples the systems (44) and/or (45) were solved using Matlab's sparse LU factorization routine.

4.1. Analytic Comparisons for PEC and Penetrable Circular Cylinders

In order to validate the DGM-TM-ME, DGM-TM-WE and DGM-TM-WH implementations we begin with TM planewave scattering from both PEC and penetrable cylinders. The incident field was selected to be an x-propagating planewave:

[mathematical expression not reproducible] (46)

where [k.sub.0] = 2[pi] f/[c.sub.0] is the free space wavenumber for frequency f, [Z.sub.0] is the impedance of free space and [c.sub.0] is the speed of light in vacuum. Analytic series solutions for PEC and penetrable cylinders can be found in [28]. For the penetrable cylinders we selected two configurations, an [epsilon]-target, having a relative permittivity of [[epsilon].sub.r] = 2 (with relative permeability [[mu].sub.r] = 1) and a [mu]-target having relative permeability [[mu].sub.r] = 1.5 (with releative permittivity [[epsilon].sub.r] = 1). The computational domains were generated using Gmsh [29] and are shown for a particular discretization in Fig. 3. The radii of the cylinder targets are 0.1061m corresponding to approximately 0.7 wavelengths in freespace at the simulation frequency of 2 GHz. The domain [OMEGA] is truncated at a radius of 0.4243 m by [[GAMMA].sup.ABC] upon which either a Silver- Muller or ERBC boundary condition was imposed. In the former case the surface [[GAMMA].sup.ERBC] was simply ignored. The PEC and penetrable meshes are respectively composed of 558 and 848 triangles, and at 2 GHz are approximately 5.5 wavelengths across. The largest edge length of any element in the meshes is 0.077 m, which is larger than a half wavelength. The cylinder boundaries are approximated using small elements to properly capture the geometries. Fig. 4 shows the real part of both the analytic and the DGM-TM-WH computed [E.sup.sct.sub.z] fields for the [epsilon]-target. The fields were computed using constant 4th order DGM formulations without ERBC mesh truncation. Visually, the solutions are indistinguishable from each other. They are also visually the same as the ME and WE formulations although their solutions are not shown. Rather than plot the fields produced for the PEC and [mu]-target cases (as they are also visually identical to the analytic solutions at 4th order), in Fig. 5 we plot the relative error between the analytic [E.sup.sct.sub.z] field for all three cylinder targets and the three DGM formulations as a function of a constant expansion order p(n) = p. Given a pth-order function f([??]) on an element [V.sub.n] it is straightforward to show that the inner product of f with itself satisfies <f, f>[v.sub.n] = [[f.bar].sup.h][M.sup.n][f.bar] where h denotes the conjugate transpose. Thus the relative [L.sup.2] error summarized in Fig. 5 is calculated as:

[mathematical expression not reproducible] (47)

where there are N total elements in the mesh and where [[E.bar].sup.sct,analytic,n.sub.z] represents the coefficients of the analytic solution expanded in the DGM basis on [V.sub.n], while [[E.bar].sup.sct,DGM,n.sub.z] contains the numerically computed DGM field coefficients. To ensure that the computed DGM solution is accurate throughout each element (and not just at the nodal points), the error quantifer (47) was not computed directly using the DGM coefficients over the simulation meshes. Instead, the error was computed by first interpolating the DGM solutions to refined meshes (not shown), consisting of N = 8, 062 and 10, 576 elements for the PEC and penetrable cylinder simulations respectively. Interpolating from the simulation meshes to the the refined meshes poses a problem as smaller elements in the refined meshes used for error calculations may fall outside of the simulation domains. To avoid having to extrapolate from the DGM solution space to the refined meshes, the latter were taken with a decreased outer radius of 0.3182 m and a slightly increased inner radius of 0.1071 m for the PEC case. Thus the errors shown in Fig. 5 are computed over a large sub-domain of [OMEGA].

For all three DGM formulations the error displays high-order convergence at low orders. The WE and WH formulations lag the DGM-TM-ME in terms of convergence due to the curl- curl term in the wave equations which requires an additional derivative computation when compared to the curl equation formulation. The 2nd-order solutions for both the WE and WH formulations are so poor (due to the curl-curl derivatives and the prescribed element sizes) that they have not been included in the figure. Further, WE outperforms WH at low orders due to the fact that that recovery of the electric field from the magnetic fields produced by the WH formulation requires an additional curl computation.

At 4th order, the fields produced by the three methods are practically indistinguishable. The benefits of the wave equation formulations can then be demonstrated by the performance summary in Table 1, which shows the reduction in degrees of freedom as well as the fill- and factor-times of the respective discrete operators for the [epsilon]-target. We note the increased factor times for the WE and WH formulations with ERBC over the ME formulation. In our experience this increase is problem-specific and the relatively small problem size considered herein may play a role. ERBCs improve the solution accuracy in all cases at the cost of the increased computational complexity required to relate all elements on the ERBC surface to all elements on the ABC surface. Mitigating this expense is possible, for example by using MLFMM [27], but is beyond the scope of this work.

For these examples, increasing the expansion order beyond p = 4 fails to improve the solution accuracy, a result that is attributed to the geometric modelling error in representing the cylinders [30]. One way to improve the results is to use h-refinement, i.e., to reduce the element size. While h-refinement can be a computational burden when the order p(n) is constant, the flexibility of DGM to support locally varying expansion orders can offset the additional cost. To validate the variable order formulations we consider the [epsilon]-target on a refined mesh consisting of 1, 778 elements and shown in Fig. 6. The number of elements within the cylinder results in an excessive oversampling of the fields and therefore we applied a locally varying order as a function of each element's electrical size. For this particular problem the order ranged from 3 through 5, resulting in a discrete system having 22, 514 degrees of freedom in the DGM-TM-WE case. By comparison, a constant 5th order WE formulation on this same mesh would require 37,338 degrees of freedom. The result of this h-refinement is also shown in Fig. 6, where we have supplemented the error convergence for the e-target shown in Fig. 5 with an additional set of data points for the variable-order solutions. The combination of ERBCs and more accurate geometric modelling achieves better error levels than the 6th order simulations on the original meshes. The refined mesh used here is an extreme case; the number of elements within the cylinder could be significantly reduced by only enforcing small element sizes in the vicinity of the cylinder boundary.

4.2. High-Order Approximations to Constitutive Parameters

An alternative geometric modelling approach for DGM formulations is to use high- order polynomial expansions for the constitutive parameters. To begin we will illustrate the concept with a continuous constitutive profile. We consider a gradient relative permittivity profile,

[mathematical expression not reproducible] (48)

in cylidrical coordinates ([rho], [phi], z) where the radius of the cylinder a = 0.1061m. We solved for the scattered fields from this gradient target using the incident fields defined in (46) and two discretization schemes. First, the refined mesh shown in Fig. 7(a), consisting of 3,966 elements, was used and a constant relative permittivity value is assigned to each element. The relative permittivity was prescribed as the value of the gradient profile at the centroid of each element. In this configuration, a 2nd order DGM-TM-ME solution was computed as a reference solution. Next, a coarse mesh consisting of 96 elements was adopted. On this mesh, shown in Fig. 7(b), a 10th-order basis expansion was assumed and the relative permittivity was prescribed directly at the underlying nodal points using Equation (48). The results of the simulations are summarized in Table 2 and show excellent agreement and significantly improved computational performance for the WE and WH formulations.

As a final elementary example we return to the scattering from the [mu]-target cylinder defined in Section 4.1. Now, however, we have approximated the exact geometry by evaluating it at the 10th order nodes of the 96 element mesh described in the previous example, as well as at the 4th-order nodes in a more refined mesh consisting of 2, 938 elements. The resulting constitutive approximations are shown in Fig. 8 and simulation results are provided in Table 3. There are two points of interest highlighted by the results. First, using high-order approximations to the constitutive parameters reduces the computational cost, demonstrated by the computational times on the 10th-order 96-element mesh. Second, and more importantly, the DGM-TM-WE formulation performs poorly. Approximating a step- function with high-order polynomials will introduce oscillations in the representation of the constitutive parameters. The DGM-TM-WE flux calculation in (28) results in a contour integral of the magnetic source vector (unlike DGM-TM-WH or DGM-TM-ME), which suggests that the effects of the oscillations along the contour corrupt the results beyond the effects of the oscillations over each volume. Although not included herein, we can confirm that for the e-target, it is the DGM-TM-WH formulation that suffers similarly.

4.3. A Practical Example of a Forward Solution for Biomedical Microwave Imaging Applications

As a final practical example, used to demonstrate the capabilities of a DGM forward solver for tomographic microwave imaging applications, we consider simulating scattering from a human breast model. The model, a realistic MRI-derived numerical phantom from the University of Wisconsin-Madison UWCEM Numerical Breast Phantom Repository [31,32], has been converted to an FEM model shown in Fig. 9. The FEM model consists of 100,584 triangles resulting in 50,482 degrees of freedom. The problem was also formulated using the DGM-TM-WE on a coarser mesh, made up of 568 triangles, using a variable-order basis with orders 5 through 7. The DGM constitutive parameters were obtained by a simple nearest-neighbour interpolation from the FEM mesh, resulting in a profile that is also shown in Fig. 9. Note that only the real part of the dielectric is shown, but that simulations include the imaginary part of the complex dielectric. The DGM-TM-WE model results in 13,628 degrees of freedom. Simulations were performed for a 1 GHz line source located at (x, y) = (0.09,0) meters in a freespace background medium. FEM and DGM simulation results are provided, for comparison, in Fig. 9. The DGM solution has 3% error with ERBCs (8% without) relative to the FEM solution, with fill and factor times not exceeding 2.5 seconds with ERBCs. The error was computed by first interpolating the DGM solution to the FEM mesh. This example demonstrates the flexibility of the proposed DGM formulations to handle geometries that are typically found in biomedical MWI applications. Improvements to the DGM model could be made by using automated local h-p refinement in regions where the constitutive parameters exhibit large gradients. However, we wish to emphasize the flexibility of the DGM formulation for efficiently representing very complicated constitutive structures without prior knowledge of their location or features. This flexibility is very attractive for imaging applications where we typically have severely limited knowledge of a target's constitutive properties.

5. CONCLUSION

The DGM formulated for Maxwell's curl equations offers flexibility that is desirable for microwave imaging applications including: local solution order refinement, support for locally varying constitutive parameters and support for both electric and magnetic constitutive parameters. To improve computational performance we have developed DGM formulations for the electric and magnetic vector wave equations when both electric and magnetic sources support the fields. The wave equation formulations also permit locally varying field order, locally varying constitutive parameters, and are compatible with exact radiating boundary conditions (ERBCs).

We have demonstrated through a number of examples that both the curl-equation DGM and the wave equation DGMs provide accurate solutions to transverse magnetic field problems in the presence of both electric and magnetic sources. ERBCs improve the overall accuracy of the scheme and provide a mechanism for validating computational models using Silver-Muller boundary conditions: if solutions with and without ERBCs do not differ substantially, Silver-Muller conditions are sufficient.

An important feature of the DGM formulations is that high-order constitutive parameters provide a mechanism for essentially decoupling the physical medium from the mesh used to support the fields. For imaging applications, this implies that we can reduce the number of degrees of freedom required to represent unknown physical parameters yet still maintain accuracy in the fields. Results have shown that when modelling discontinuous constitutive parameters using the DGM polynomial basis, one should be careful in choosing a formulation. For example, the accuracy of the inhomogeneous magnetic vector wave equation suffers for discontinuous dielectric targets when the mesh does not align itself with the discontinuity. We are currently investigating a hybrid formulation that applies, on an element-by-element basis, either the curl-equation formulation or a wave equation formulation. Details on such a formulation are forthcoming.

ACKNOWLEDGMENT

The authors wish to thank Anastasia Baran from the University of Manitoba for helpful discussions and for generating the FEM solution to the breast phantom problem. We also thank the Natural Sciences and Engineering Research Council of Canada for supporting this work.

APPENDIX A. NODAL DGM OPERATORS

The purpose of this appendix is to provide details for deriving the tested DGM formulations for the wave equations. We introduce the differentiation operator [D.sup.n.sub.[zeta]] defined such that [D.sup.n.sub.[zeta]][[bar.E].sup.n.sub.[xi]] gives the coefficients for the local function [partial derivative][E.sub.[xi]]/[partial derivative][zeta] in [V.sub.n]. It is straightforward to show that [1]

[mathematical expression not reproducible] (A1)

and it follows that the coefficients for the curl of the electric field can be computed as

[mathematical expression not reproducible] (A2)

This operator can be used to compute the flux terms in (28) that depend on the curl of [[??].sup.sct]. Similarly, the integrals [integral] [psi][nabla] x [[??].sup.sct]dv required for computing (7), where [psi] takes on all expansion functions [E.sup.n.sub.m], can be written as:

[mathematical expression not reproducible]. (A3)

where it can be shown that [M.sup.n][D.sup.n.sub.[zeta]] = [S.sup.n.sub.[zeta]] [1]. In (8), the stiffness terms for the electric field wave equation involve [integral] [psi][nabla] x ([[mu].sup.-1][nabla] x [[??].sup.sct])dv. It follows from (A2) that the coefficients for the curl-curl of the electric field in [V.sub.n] are (suppressing n on [D.sub.[zeta]]):

[mathematical expression not reproducible] (A4)

and it follows that the tested values can be computed as:

[mathematical expression not reproducible] (A5)

Expanding the right-hand-side of this expression gives the leading matrix in (14).

REFERENCES

[1.] Hesthaven, J. S. and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis and Applications, Springer, New York, 2008.

[2.] Hesthaven, J. S. and T. Warburton, "Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell's equations," J. Comput. Phys., Vol. 181, No. 1, 186-221, 2002.

[3.] Jeffrey, I., "Finite-volume simulations of Maxwell's equations on unstructured grids," Ph.D. Dissertation, University of Manitoba, 2011.

[4.] Liu, M., K. Sirenko, and H. Bagci, "An efficient discontinous Galerkin finite element method for highly accurate solution of Maxwell's equations," IEEE Trans. Antennas. Propag., Vol. 60, No. 8, 3992-3998, 2012.

[5.] Shi, Y. and C.-H. Liang, "Simulations of the left-handed medium using discontinuous Galerkin method based on the hybrid domains," Progress In Electromagnetics Research, Vol. 63, 171-191, 2006.

[6.] Buffa, A. and I. Perugia, "Discontinuous Galerkin approximation of the Maxwell eigenproblem," SIAM J. Numer. Anal. , Vol. 44, No. 5, 2198-2226, 2006.

[7.] Warburton, T. and M. Embree, "The role of the penalty in the local discontinuous Galerkin method for Maxwell's eigenvalue problem," Comput. Method Appl. Mech. Eng., Vol. 195, No. 25, 3205-3223, 2006.

[8.] Li, L., S. Lanteri, and R. Perrussel, "A hybridizable discontinuous Galerkin method for solving 3D time-harmonic Maxwell's equations," Numerical Mathematics and Advanced Applications, 119-128, Springer, Berlin Heidelberg, 2013.

[9.] Bouajaji, M. E. and S. Lanteri, "High order discontinuous Galerkin method for the solution of 2D time-harmonic Maxwell's equations," Appl. Math. Comput., Vol. 219, No. 13, 7241- 7251, 2013.

[10.] Arnold, D. N., F. Brezzi, B. Cockburn, and L. D. Marini, "Unified analysis of discontinuous Galerkin methods for elliptic problems," SIAM J. Numer. Anal., Vol. 39, No. 5, 1749-1779, 2002.

[11.] Gabard, G., "Discontinuous Galerkin methods with plane waves for time- harmonic problems," J. Comput. Phys., Vol. 225, No. 2, 1961-1984, 2007.

[12.] Perugia, I., D. Schotzau, and P. Monk, "Stabilized interior penalty methods for the time-harmonic Maxwell equations," Comput. Method Appl. Mech. Eng., Vol. 191, No. 41, 4675- 4697, 2002.

[13.] Lohrengel, S. and S. Nicaise, "A discontinuous Galerkin method on refined meshes for two-dimensional time-harmonic Maxwell equations in composite materials," J. Comput. Appl. Math., Vol. 206, No. 1, 27-54, 2007.

[14.] Hiptmair, R., A. Moiola, and I. Perugia, "Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations," Math. Comput., Vol. 82, No. 281, 247- 268, 2013.

[15.] Jin, J., The Finite Element Method in Electromagnetics, Wiley, New York, 2002.

[16.] Dolean, V., H. Fol, S. Lanteri, and R. Perrussel, "Solution of the time- harmonic Maxwell equations using discontinuous Galerkin methods," J. Comput. Appl. Math., Vol. 218, No. 2, 435-445, 2008.

[17.] Zakaria, A., I. Jeffrey, and J. LoVetri, "Full-vectorial parallel finite- element contrast source inversion method," Progress In Electromagnetics Research/, Vol. 142, 463-483, 2013.

[18.] Zakaria, A., A. Baran, and J. LoVetri, "Estimation and use of prior information in FEM-CSI for biomedical microwave tomography," IEEE Antennas Wireless Propag. Lett., Vol. 11, 1606-1609, 2012.

[19.] Bellizzi, G., O. M. Bucci, and I. Catapano, "Microwave cancer imaging exploiting magnetic nanoparticles as contrast agent," IEEE Trans. Biomed. Eng., Vol. 58, No. 9, 2528-2536, 2011.

[20.] Hanson, G. W. and A. B. Yakovlev, Operator Theory for Electromagnetics: An Introduction, Springer, New York, 2002.

[21.] Bonnet, P., X. Ferrieres, B. L. Michielsen, P. Klotz, and J. L. Roumiguires, Finite-volume Time Domain Method, 307-367, Academic Press, San Diego, 1999.

[22.] Sankaran, K., C. Fumeaux, and R. Vahldieck, "Cell-centered finite-volume-based perfectly matched layer for time-domain Maxwell system," IEEE Trans. Microw. Theory Techn., Vol. 54, No. 3, 1269-1276, 2006.

[23.] Dosopoulos, S. and J.-F. Lee, "Interior penalty discontinuous Galerkin finite element method for the time-dependent first order Maxwell's equations," IEEE Trans. Antennas Propag., Vol. 58, No. 12, 4085-4090, 2010.

[24.] Pearson, L., R. Whitaker, and L. Bahrmasel, "An exact radiation boundary condition for the finite-element solution of electromagnetic scattering on an open domain," IEEE Trans. Magn., Vol. 25, No. 4, 3046-3048, 1989.

[25.] Firsov, K. D. and J. LoVetri, "FVTD-integral equation hybrid for Maxwell's equations," Int. J. Numer. Model. El., Vol. 21, No. 1-2, 29-42, 2008.

[26.] Ziolkowski, R. W., N. K. Madsen, and R. C. Carpenter, "Three-dimensional computer modeling of electromagnetic fields: A global lookback lattice truncation scheme," J. Comput. Phys., Vol. 50, No. 3, 360-408, 1983.

[27.] Shanker, B., M. Lu, A. A. Ergin, and E. Michielssen, "Plane-wave time- domain accelerated radiation boundary kernels for FDTD analysis of 3-D electromagnetic phenomena," IEEE Trans. Antennas Propag., Vol. 53, No. 11, 3704-3716, 2005.

[28.] Harrington, R. F., Time-harmonic Electromagnetic Fields, 2nd Edition, Wiley-Interscience, New York, 2001.

[29.] Geuzaine, C. and J.-F. Remacle, "GMSH: A 3-D finite-element mesh generator with built-in pre- and post-processing facilities," Int. J. Numer. Meth. Eng., Vol. 79, No. 11, 1309- 1331, 2009.

[30.] Jeffrey, I., J. Aronsson, M. Shafieipour, and V. Okhmatovski, "Error controllable solutions of large-scale problems in electromagnetics: MLFMA-accelerated locally corrected Nystrom solutions of CFIE in 3D," IEEE Antennas Propag. Mag., Vol. 55, No. 3, 294-308, 2013.

[31.] Lazebnik, M., et al., "A large-scale study of the ultrawideband microwave dielectric properties of normal, benign and malignant breast tissues obtained from cancer surgeries," Phys. Med. Biol., Vol. 52, No. 20, 6093, 2007.

[32.] Burfeindt, M. J., et al., "MRI-derived 3-D-printed breast phantom for microwave breast imaging validation," IEEE Antennas Wireless Propag. Lett., Vol. 11, 1610-1613, 2012.

Ian Jeffrey *, Nicholas Geddert, Kevin Brown, and Joe LoVetri

Received 4 September 2015, Accepted 23 October 2015, Scheduled 6 November 2015

* Corresponding author: Ian Jeffrey (Ian.Jeffrey@umanitoba.ca).

The authors are with Department of Electrical and Computer Engineering, University of Manitoba, Winnipeg, MB, Canada.

([dagger]) The notation * is often used instead of [conjunction] for time-domain formulations [1]. For time-harmonic formulations, where * could easily be confused with complex conjugation, we have substituted [conjunction].

([double dagger]) Note that in [1] the notation {{Z}} is used to denote the average of [Z.sup.-] and [Z.sup.+] and an additional factor of 1/2 is included in the flux expressions.

([section][section]) While we do not rigorously prove this simplification the numerical results support the claim.

Caption: Figure 1. ERBC configuration. All sources are contained within [[GAMMA].sup.ERBC]. The zoomed view shows the interior elements used to interpolate to quadrature points on the ERBC surface (red). The Huygens' intergral is used to evaluate [[??].sup.sct,+] and [[??].sup.sct,+] on [[GAMMA].sup.ABC] (blue).

Caption: Figure 2. The interpolation-permutation framework used to compute flux contributions from neighbouring elements when the solution order varies from [V.sub.n] to [V.sub.n']. The first and second order elements also demonstrate the preservation of order when moving from an element to its boundary.

Caption: Figure 3. Meshes used for simulating planewave scattering from (a) a PEC cylinder and (b) penetrable cylinders. The meshes are the same for SM and ERBC boundary conditions. The ERBC Huygens' surface is shown in red, while the boundary of the penetrable cylinders is shown in black.

Caption: Figure 4. Comparison of (a) the real part of [E.sup.sct.sub.z] computed analytically and (b) using the DGM-TM- WH (magnetic wave equation) formulation at constant order p = 4.

Caption: Figure 5. Relative errors in Real ([E.sup.sct.sub.z]) for the ME, WE and WH formulations, with and without ERBCs, as a function of constant expansion order p for the cylinder targets: (a) PEC, (b) [epsilon]-target and (c) [mu]-target.

Caption: Figure 6. (a) Refined mesh for the e-target upon which a variable order from p = 3 to p = 5 is imposed based on element size. (b), the relative errors obtained from the variable-order DGM formulations are appended (as the last data point after order p = 6) to the [epsilon]-target error plot for the constant-order mesh, originally shown in Fig. 5.

Caption: Figure 7. Gradient dielectric target discretized as (a) a constant value based on the centroid position of each element in a fine mesh and (b) on a coarse mesh at each of the 10th- order nodes in each element.

Caption: Figure 8. The [mu]-target approximated using polynomials on (a) the coarse 96 element mesh using p = 10 and (b) a 2,938 element mesh using p = 4.

Caption: Figure 9. (a) A University of Wisconsin-Madison UWCEM breast model target, modelled using FEM and (b) DGM with a variable-order basis. (c) The fields produced by FEM and (d) DGM-TM-WE agree. White boxes denote the location of the parameters shown.

Table 1. Comparison of the DGM solutions for the e-target at 4th order, including the dimension of the matrix [dim.sub.K], the number of non-zeros (nnz) in the system, the fill-time [t.sub.fill] and the factor time [t.sub.factor]. Order p Element [dim.sub.K] [nnz.sub.K] Count DGM-TM-ME 4 848 38,160 1,973,384 DGM-TM-ME + ERBC 4 848 38,160 2,640,284 DGM-TM-WE 4 848 12,720 378,750 DGM-TM-WE + ERBC 4 848 12,720 452,850 DGM-TM-WH 4 848 25,440 1,514,615 DGM-TM-WH + ERBC 4 848 25,440 1,811,015 [t.sub.fill] (s) [t.sub.factor] (s) DGM-TM-ME 2.7 1.9 DGM-TM-ME + ERBC 4.4 3.4 DGM-TM-WE 1.8 1.0 DGM-TM-WE + ERBC 2.4 5.0 DGM-TM-WH 2.5 6.8 DGM-TM-WH + ERBC 3.7 11.0 Table 2. Summary of gradient-target results comparing the refined 2nd order reference solution to the 10th order variable permittivity approximation. Errors are computed relative to the reference solution. Order p Element [dim.sub.K] [nnz.sub.K] Count DGM-TM-ME 2 3,966 71,388 2,061,232 Reference DGM-TM-ME 10 96 19,008 3,292,750 DGM-TM-WE 10 96 6,336 615,648 DGM-TM-WH 10 96 12,672 2,461,782 [t.sub.fill] (s) [t.sub.factor] (s) Relative Error DGM-TM-ME 7.3 3.4 -- Reference DGM-TM-ME 3.3 4.2 0.0065 DGM-TM-WE 0.8 0.28 0.0064 DGM-TM-WH 2.2 0.95 0.0093 Table 3. Summary of simulations for the locally varying constitutive approximation of the [mu]-target for a coarse and fine mesh. Errors are computed relative to the analytic solution. Order p Element [dim.sub.K] [nnz.sub.K] Count DGM-TM-ME 4 2,938 132,210 6,839,268 DGM-TM-WE 4 2,938 44,070 1,315,650 DGM-TM-WH 4 2,938 88,140 5,262,600 DGM-TM-ME 10 96 19,008 3,294,592 DGM-TM-WE 10 96 6,336 615,648 DGM-TM-WH 10 96 12,672 2,461,782 [t.sub.fill] (s) [t.sub.factor] (s) Relative Error DGM-TM-ME 9.4 8.1 0.026 DGM-TM-WE 6.1 4.4 0.069 DGM-TM-WH 8.9 42.0 0.027 DGM-TM-ME 3.3 3.9 0.026 DGM-TM-WE 0.87 0.24 0.072 DGM-TM-WH 2.3 0.95 0.028

Printer friendly Cite/link Email Feedback | |

Author: | Jeffrey, Ian; Geddert, Nicholas; Brown, Kevin; LoVetri, Joe |
---|---|

Publication: | Progress In Electromagnetics Research |

Date: | Sep 1, 2015 |

Words: | 9147 |

Previous Article: | Broadband Nanoantennas for Plasmon Enhanced Fluorescence and Raman Spectroscopies. |

Next Article: | Plane-Wave Propagation in Electromagnetic PQ Medium. |

Topics: |