A High-Order Ultra-Weak Variational Formulation for Electromagnetic Waves Utilizing Curved Elements

,


Introduction
A Trefftz method [24] for approximating a linear partial differential equation is a numerical method using local solutions of the underlying partial differential equation as basis functions.The version we shall study here, the Ultra Weak Variational Formulation (UWVF) of Maxwell's equations, is a Trefftz type method for approximating the solution of Maxwell's equations on a bounded domain due to O. Cessenat and B. Després [5,6].The UWVF uses a finite element computational grid, classically composed of tetrahedral elements, and plane wave solutions of Maxwell's equations on each element.
A study of this method from the point of view of symmetric hyperbolic systems was presented in [13] where the inclusion of the Perfectly Matched Layer absorbing condition [3] into the UWVF was also described.In addition, several computational heuristics relevant to our study were also presented.As a result of this work, a first parallel implementation of the UWVF called ParMax was written.This was further developed at Kuava Inc, and the University of Eastern Finland and is the basic software used in this paper.ParMax uses MPI and domain decomposition to implement parallelism.
From the point of view of theoretical convergence analysis, it was shown in [4] that the UWVF for the related Helmholtz equation is a special discontinuous Galerkin (DG) method, and using duality theory convergence estimates could be proved.A more general DG approach, again for the Helmholtz equation, is taken in [7].Based on these studies, and following the proof of explicit stability bounds for the interior impedance boundary value problem for Maxwell's equations in [10], convergence of Trefftz DG methods (in particular the UWVF) for Maxwell's equations under rather strict geometric assumptions was proved in [11].A key point in this analysis is that it is not necessary to use tetrahedral elements, but a wide class of element shapes are theoretically covered and we shall return to this point in Section 3.
A restriction on the use of the UWVF is that the relative electric permittivity denote ϵ r and relative magnetic permeability µ r must be piecewise constant (constant on each tetrahedron in the mesh).This assumption could perhaps be weakened using generalised plane waves [14] or embedded Trefftz techniques [17] but these have not yet been extended to Maxwell's equations (the latter technique has been demonstrated in the open source package [23]).The restriction of piecewise constant media is relaxed if a Perfectly Matched Layer (PML) is used.There ϵ r and µ r are spatially varying tensors [13].
The main issue facing the UWVF is that the condition number of the global system rises rapidly as the number of plane wave directions increases.This in turn causes the iterative solution of the linear system to slow down.We adopt the approach of [13] choosing the number of plane wave directions on a given element on the basis of its geometric size (in wavelengths) to control ill-conditioning.
We use the same stabilised BiConjugate Gradient (BICGstab) scheme as in [13] for solving the global linear system resulting from the UWVF.Interesting new results on preconditioned iterative schemes can be found in [20].These results are based on the use of a structured hexahedral grid, whereas we focus on unstructured grids in this paper.
Within the UWVF scheme used here there is considerable scope for different choices of basis functions provided they form a complete family of solutions of Maxwell's equations.For example, plane wave basis functions (as used here), vector wave functions or the method of fundamental solutions (MFS) are discussed in [2].In the interesting case of MFS, [8] gives an application to wave guides.to the method.Our choice of plane waves in ParMax follows [5] and has the advantage that necessary integrals of products of basis functions can be performed in closed form on flat triangular faces in the mesh (the majority of faces).This greatly speeds up the assembly of matrices compared to quadrature that is needed for other choices of basis functions.
In this paper, i = √ −1 and we use the convention that the time variation of the fields and sources is proportional to exp(−iωt) where ω is the angular frequency of the radiation, and t is time.All results are then reported in the frequency domain.Bold face quantities are vector valued.The coefficients ϵ 0 and µ 0 are the permittivity and permeability of free space respectively.The wave number κ of the radiation is given by κ = ω √ ϵ 0 µ 0 .
To further fix notation and context, we now define the Maxwell system under consideration in this paper.Let Ω denote a Lipschitz bounded computational domain having unit outward normal ν and boundary Γ := ∂Ω.For a smooth enough vector field v, we define the tangential component v T on Γ by v T = (ν × v) × ν.Then, given the wave number κ > 0, a tangential boundary vector field g, piecewise constant functions ϵ r and µ r , and a parameter Q ∈ C with |Q| ≤ 1 we seek the complex valued vector electric field E that satisfies Here Z is the surface impedance (a positive real parameter).The boundary condition (1b) is of impedance type and is well suited to the UWVF.When Q = −1 it gives a rotated version of the Perfectly Electrically Conducting (PEC) boundary condition for the scatterered wave where we take g = −2iκE i T /Z and E i is the incident wave.When Q = 0 it corresponds to an outgoing condition that can be used as a low order radiation condition, while for Q = 1 we have a symmetry boundary condition.
This paper presents several novel extensions of the basic UWVF that are of considerable utility in practical applications.In particular: • The original UWVF uses a tetrahedral grid, however the error estimates in [11] hold for more general element types.Besides tetrahedral elements, we have implemented hexahedral and wedge elements.In this paper, hexahedral elements are only used in the PML region.
• Very often curved surfaces appear in applications, and we have implemented a mapping technique to approximate smooth curved surfaces.This allows us to use larger elements near a curved boundary.We shall show, using numerical experiments, that this improves the efficiency of the software by decreasing the overall time to compute a solution.Note that we only need to map faces in the mesh which simplifies the implementation, but we have to use quadrature to evaluate integrals on curved faces.
• We show how to implement resistive sheet transmission conditions across thin interfaces.
Related to this we have also implemented a combined total field and scattered field formulation to allow the solution of problems involving penetrable media.
• We point out that a low memory version can be used to solve very large problems by avoiding the storage of the most memory intensive matrix in the algorithm.
• We provide numerical results to justify the utility of the above innovations.
The paper proceeds as follows.In Section 2 we start with a brief derivation of the basic UWVF and describe the plane wave based UWVF.Then in Section 3 we describe the five contributions of this paper.We start with comments on new geometric element types and a brief discussion of implementing an algorithm for using scattered or total fields in different subdomains in the context of scattering from a penetrable object.Then we move to discuss resistive sheets, curved elements and quadrature, and finally a lower memory version of the method.Section 4 is devoted to numerical examples illustrating the UWVF and the previously mentioned modifications.In Section 4.1 we start with two examples of scattering from a PEC scatterer.The first PEC example is scattering from a sphere for which the Mie series gives an accurate solution for comparison (c.f.[19]).This example demonstrates the benefits of curved elements and different element types, as well as the use of very coarse meshes (compared to those used by finite element methods).The second example is X-band scattering from an aircraft model.Here again we use curved elements, but in addition use the low-memory version of the software.In Section 4.2 we consider two examples involving a resistive sheet.The first is a classical example having an exact analytical solution, and the second is a resistive screen surrounding a sphere.Next, in Section 4.3, we consider heterogeneous or penetrable scatterers.The first example is a dielectric sphere where we use the scattered/near field formulation and compare to the Mie series solution.The second example is scattering from a plasma.In Section 5 we present our conclusions.Finally in Appendix A we give an update to the basis selection rule used previously in [13] for the new element types.Then, in Appendix B, we present a comparison of ParMax with the edge finite element method for scattering from a dielectric sphere.

Derivation and properties of the basic UWVF
In this section we provide a sketch of the derivation of the UWVF sufficient to allow us to present the new features of this paper in the following section.For full details see [5,13].

A brief derivation of the UWVF
The version of the UWVF presented here is equivalent to that used in ParMax (from [5]) but with simplified notation.Consider a mesh of Ω of elements of maximum diameter h denoted by T h .An element K ∈ T h in this mesh is a curvilinear polyhedron (curvilinear tetrahedron, wedge, or hexahedron) with boundary denoted by ∂K and unit outward normal ν K .We now extend the parameter Z to a real piecewise positive constant defined on all faces in the grid.Following [5], we choose Z as follows.Let where | K denotes the restriction to K. The edge function μ is defined in the same way.Then Z = √ μ/ √ ε.Suppose ξ is a smooth solution of the adjoint Maxwell equation in K: Then taking the dot product of (1a) with ξ (including complex conjugation) and integrating by parts twice provides the following fundamental relation between the electric and magnetic fluxes on ∂K: Using the above fundamental identity, we can then prove equality (4) by expanding both sides of (4) and using (3).Equality (4) gives the conclusion of the "isometry lemma" (c.f.[5]).
To simplify the presentation, we define rescaled versions of the unknowns in [5] as follows: The next step is to rewrite (4) using the above quantities.In doing so we will use the function space of surface vectors Recalling that ξ| K satisfies the adjoint Maxwell system in K we can define T .Now suppose elements K and K ′ meet at a face in the mesh.On that face ν K = −ν K ′ .Also we have transmission condition requiring continuity of E T and (µ −1 r ∇ × E) T across the face so For a boundary face we can use the boundary condition (1b) to replace the corresponding term on that face in terms of χ K and g.Using the above results, we may rewrite (4) for every element K.We conclude that for χ K ∈ L 2 T (∂K) the following equation holds for any Y K ∈ L 2 T (∂K).This equation should hold for every K ∈ T h , and gives the UWVF for Maxwell's equations before discretization.Cessenat [5] proves the uniqueness of the UWVF solution, and existence follows because (1a)-(1b) is well posed.

The plane wave UWVF (PW-UWVF)
Now that we have the variational formulation (8) we can discretize it by using a subspace of L 2 T (∂K) on each element.It is important for efficiency that F K be easy to compute and this is where the plane wave basis is useful [5].On each element K we choose independent directions {d K,j } p K j=1 , ∥d K,j ∥ = 1, using the first p K Hammersley points [9] on the unit sphere.Then we choose ξ| K to be a linear combination of the p K plane wave solutions of the adjoint problem and K ∈ T h .Here x K,0 is the centroid of the element and the polarizations A K,ℓ,m are chosen to be unit vectors such that Then W h,p = Π K∈T h W K,h,p K where p denotes the vector of number of directions on each element.
The dimension of W h,p is N DoF = 2 K∈T h p K .In our work, p K is chosen according to the heuristic in [13] (see Appendix A for updates to this formula for larger numbers of directions and new element types).Then F K is easy to compute using the definition of the basis functions.The discrete PW-UWVF uses trial and test functions from W h,p in place of W in (8).
In our implementation, we use domain decomposition by subdividing the mesh according to Metis [16] where p K is used to estimate the work on each element.The elements are sorted using reverse Cuthill-McKee to minimise bandwidth.Then, enumerating the degrees of freedom element by element, the matrices and vectors corresponding to the terms in ( 8) can be computed.The left hand side of ( 8) gives rise to an N DoF × N DoF block diagonal Hermitian positive-definite matrix D, while the remaining sesquilinear forms on the right hand side of (8) give rise to a general sparse complex matrix C. The data term involving g gives rise to a corresponding vector ⃗ b.Denoting the vector of unknown degrees of freedom by ⃗ χ, we solve the global matrix equation using BICGstab where D −1 can be calculated rapidly element by element.Once ⃗ χ is known, the solution on each element can be reconstructed for post processing.For unbounded scattering problems, we use either the low order absorbing condition (Q = 0) on the outer boundary or a PML as in [13].For a scattering problem the far field pattern can the be calculated using an auxiliary surface containing all the scatterers in its interior in the usual way [5].

Towards industrial scale software
In this section we discuss the extensions to the basic UWVF in this paper.

New element types
In [11] error estimates are proved for general elements that have Lipschitz boundaries, are shape regular in the sense of that paper, and are star-shaped with respect to a ball centred at a point in the element.This allows a wide variety of elements.We have implemented curvilinear tetrahedral, wedge, and hexahedral elements.In ParMax we represent the boundary of each element as a union of possibly curvilinear triangles.The assembly phase is quickest if the triangles are planar since then quadrature can be avoided.
For this paper, we rely on COMSOL Multiphysics to generate the meshes.An example of a grid using tetrahedral, wedge, and hexahedral elements is shown in Fig. 3.Here we use wedge and hexahedral elements in the PML, and tetrahedral elements elsewhere.

Scattered/total field formulation
The numerical results we shall present are all of scattering type.The total field E is composed of an unknown scattered field E s and a given incident field E i so E = E s + E i .We will use plane wave incident fields (but point sources or other incident fields can be used) so where d ∈ R 3 is the direction of propagation and ∥d∥ = 1.The vector polarisation p ∈ C 3 is non-zero and satisfies d • p = 0.
To allow the use of a PML or other absorbing boundary condition, we need to compute using the scatterered field in the PML.But inside a penetrable scatterer we need to compute with the total field so that there are no current sources in the scatterer.This is standard for finite element methods, but not usual for the UWVF so we outline the process here.Suppose Ω is partitioned into two subdomains denoted Ω − and Ω + such that the PML (or a neighbourhood of the absorbing boundary if one is used) is contained in Ω + where the scattered field is used and where ϵ r = µ r = 1.The scatterer is contained in Ω − where possible ϵ r or µ r are no longer unity and the total field is used.Let Σ denote the boundary between Ω − and Ω + and assume that Σ is contained in the interior of Ω and is exactly covered by faces of the mesh.Then suppose that elements K − ⊂ Ω − and K + ⊂ Ω + meet at a face F ⊂ Σ.
Consider first K − .Reviewing the derivation of the UWVF outlined in Section 2.2 we see that we must rewrite the term in E on the right hand side of equation ( 4) in terms of χ K − and χ K + .In particular using the facts that ν K − = −ν K + and that on K + the scattered field E s is be approximated, so that in (5) Here we have also used the transmission condition that the tangential components of E and µ −1 r curl E are continuous across Σ.In this equation, the source function on F associated with Carrying out the same procedure on K + remembering this element supports the scattered field gives another equation and source function g K + ,F again relating E and E s .Thus a combined scattered and total field algorithm can be implemented by allowing for a source function g on the internal surfaces.

Resistive sheet
A similar procedure to that used to implement the scattered/total field UWVF can be used to derive appropriate modifications to include resistive [15] or conductive sheets [22].We only consider the resistive sheet.Suppose now that a surface Σ r ⊂ Ω is a resistive sheet, and that ν r denotes a continuous normal to the sheet.This surface may be open or closed and could intersect the boundary.We just assume that it is the union of a subset of the faces (possibly curvilinear) in the mesh.Suppose K + and K − are two elements that meet at a face F r ⊂ Σ r such that ν r points into K + .The geometry is shown in Fig. 1.

The resistive sheet approximation requires us to implement the following transmission conditions across Σ
where σ is the conductivity of the material in the layer and d is the thickness.Alternatively the resistivity of the sheet is given by with ϵ r = 1 + iσ/(ϵ 0 ω), Z 0 = µ 0 /ϵ 0 is the impedance of free space and κ 0 = ω √ µ 0 ϵ 0 is the wave-number.We now set η = σd.As for the scattered/total field version of UWVF we must write the term in E on the right hand side of (4) using the variables χ| K + and χ| K − on either side of F r .
Adding χ| K + to χ| K − , using the resistive sheet transmission conditions and noting that ν Now we rewrite the term in E on the right hand side of (4) using the UWVF functions.Using Equation (13) we have But, using (11), we have Then, using the above equality in (14) together with (13) again we have The new term introduces a new diagonal block into the matrix C (defined before (9)) and a perturbation to the off diagonal blocks coupling fields on K + and K − .In the current ParMax implementation, η is assumed constant on each mesh face of the resistive sheet and may be complex.

Curved elements and quadrature
We take a straightforward approach to approximating curved boundaries and quadrature.All the elements in ParMax have faces that are unions of possibly curvilinear triangles.Suppose Figure 2: Sketch of mapping from the reference face to the face of an element in the volume mesh.Here we sketch a quadratic map requiring that the midpoint of each edge in the curvilinear face be given.a curvilinear face F in the mesh is such that either an edge of F , or F itself are entirely contained in a smooth curvilinear subset of the boundary Γ.We approximate F by a mapping from a reference element F (with vertices a 1 = (0, 0), a 2 = (0, 1), and a 3 = (1, 0)) in the (s, t) plane to an approximation of F using a degree ℓ polynomial map F F : F → F .
In the important case of a quadratic map (ℓ = 2) we choose a i , i = 1, 2, 3 to be the vertices of K and take the remaining interpolation points by choosing a i,j to be a point on the smooth boundary approximately half way between a i and a j (mapped from a i,j = ( a i + a j )/2; see Fig. 2).We use the following nodal basis and define F F by We use F F ( F ) ≈ F for computing the UWVF matrices.Thus we need to evaluate integrals over each face F F ( F ). Using the reference element, for a smooth function g defined in a neighbourhood of It now suffices to define quadrature on the reference element via the Duffy transform.Suppose G(s, t) is a smooth function on F (in particular the integrand above) then setting t = (1 − s)ξ we can map to the unit square: Let (w g i , t g i ), 1 ≤ i ≤ N , denote the N -point Gauss-Legendre rule weights and nodes on (0, 1).Then for each s Next, using N point Jacobi quadrature weights and nodes on (0, 1) denoted (w J j , x J j ), we obtain finally Note that the quadrature has positive weights.

A low memory version
A simple low memory version of ParMax is easily available because we use BICGstab to solve the linear system.We compute D as usual (it is block diagonal), and the vector ⃗ b but do not compute the elements of C in (9).Then, as required by BICGstab, to compute D −1 C⃗ x for some vector ⃗ x we compute the blocks of C element by element and accumulate C⃗ x element by element.Then D −1 is computed element by element using precomputed LU decompositions.
Obviously computing the entries of C repeatedly greatly increases CPU time but this allows us to compute solutions to problem that would otherwise require very large memory to store C.For example, the solution of a scattering problem for a full aircraft at X-band frequencies is shown in Section 4.1.2.

Numerical examples
All results were generated using the computer clusters Puhti and Mahti at the CSC -IT Center for Science Ltd, Finland.Detailed descriptions of these supercomputers can be found from the CSC's website [1].Computational grids used in this work were prepared using COMSOL Multiphysics on a personal computer.In addition, the geometry model for aircraft used in Section 4.1.2is adapted from the COMSOL's application Simulating Antenna Crosstalk on an Airplane's Fuselage.For all numerical experiments, the incident electric field is a plane wave propagating in the direction of the positive x-axis polarised in the y-direction where the field is given by (10).

Scattering from PEC objects
In both PEC examples, we compute the scattered field and the incident field is used as a source via the PEC condition on the surface of the scatterer.

A sphere
This experiment is intended to demonstrate the advantages of multiple element types, and a cuvilinear approximation to a smooth curved boundary (the surface of the sphere).In particular, we study scattering from a PEC sphere placed in vacuum, ϵ r = µ r = 1, where the frequency of the incident field is f = 2 GHz (wavelength in air: λ 0 = 0.19946 m).The scatterer is a PEC sphere with radius of 1 meter or approximately 5λ 0 .
In order to demonstrate the use of several types of large elements and curvilinear grids, we use an unusually large computational domain.In particular, the sphere is placed with its origin at the center of a cube shaped computational domain [−1 − 15λ 0 , 1 + 15λ 0 ] 3 .
An absorbing boundary condition, (1b) with Q = 0, is used on the exterior surface geometry.In addition, a PML with thickness of 5λ 0 , and constant absorption parameter σ 0 = 1 is used within each side of the cube.
Two computational grids with different geometric approximations are used (see Fig. 3).For the first grid (mesh 1 ), we set the requested element size on the PEC sphere to be h s = λ 0 /5 which we will see provides a geometrically accurate surface representation using flat face elements such that the far field pattern predicted by UWVF is in good agreement with the far field computed via Mie scattering.For the second grid (mesh 2 ), the surface grid density is relaxed to h s = 3λ 0 .For both cases, the requested element size in the volume is set to h v = 10λ 0 .We request wedge and hexahedral elements in the PML region, and tetrahedra elsewhere.
We show results for three cases: 1) scattering calculated using mesh 1, 2) scattering computed using mesh 2 with flat facets approximating the surface of the sphere, and 3) the use of mesh 2 with a quadratic approximation to the boundary of the sphere (see Section 3.4).These results are computed on the Puhti system with 5 computing nodes and 10 cores per node.Table 1 summarises the notation used in reporting results and Table 2 gives details of the PEC computations including CPU time.
Figure 3 shows that up to 1,200 directions are used on some elements.As is done in the Mie series, the field scattered by a sphere can be approximated well by relatively few spherical vector wave functions, and these in turn can be approximated by special plane wave expansions involving many less directions [18].However our code is for general wave propagation problems, and the heuristic used in the appendix for giving the number of directions on a element always chooses the largest number consistent with a chosen condition number.This approach is intended to ensure good accuracy (within the conditioning constraint) for a general problem.4 shows the modulus of the y-component of the scattered electric field |E s y | on the z = 0 plane.There are clear differences between the results for mesh 1 and mesh 2 with flat facets.These are caused by the coarse surface grid in the second case.However mesh 1 with flat facets, and mesh 2 with a quadratic boundary approximation are in good agreement.As can be seen in Table 2, mesh 2 with quadratic boundary approximation is much cheaper in terms of CPU-time than mesh 1.
Very often the far field pattern of the scattered wave [19] is the quantity of interest for these calculations, and in particular the Radar Cross Section (RCS) derived from the far field pattern.In this paper, far field directions are defined in terms of the azimuth angle ϕ ( • ) as (cos(ϕπ/180), sin(ϕπ/180), 0). Figure 3: Cross-sections of the computational grids used to approximate scattering from a PEC sphere.The colorbar shows the number of plane wave directions on each element.The major difference is the grid density on the surface of the sphere.This figure shows how wedge and hexahedral elements can be usefully employed in the outer PML layer.In all figures showing grids, the colorbar shows the number of plane wave direction, see Eq. ( 16), for each element.
Figure 4: Snapshots of the scattered electric field component |E s y | for the meshes considered here.There is good agreement between mesh 1 and mesh 2 with curved faces.Mesh 2 with flat faces produces unacceptable error due to the coarse boundary approximation.
Figure 5 shows a comparison of the bistatic RCS predicted by the computational experiments with the UWVF and one computed by the Mie series.Clearly mesh 2 with flat facets produces an inaccurate far field pattern, whereas mesh 1 or mesh 2 with curved elements produces much more accurate predictions.

An aircraft at X-band frequency
The aircraft model used in this section is derived from a model available in COMSOL (application Simulating Antenna Crosstalk on an Airplane's Fuselage).We treat the aircraft as a curvilinear perfect conductor.The frequency of the incident field is f = 8 GHz so λ 0 = 0.03737 Figure 5: Bistatic RCS at 2 GHz for the PEC sphere.We show the RCS computed using meshes 1 and 2 compared to the Mie series.The bottom panel shows the difference between the numerical solution and the Mie series.Note that the 'Mesh 2, flat -Mie' curve is omitted from the bottom panel due to its significantly larger amplitude compared to the other two curves.For generating the computational grid, curved elements and a mesh size parameter of h s = 3λ 0 were employed on the aircraft's surface.Because we can use 10λ 0 sized elements away from the boundary, the entire grid can be created on a standard office computer.The grid consists of 697,783 tetrahedral elements with 142,731 vertices covering the computational domain.The surrounding PML layer is discretized using 413 hexahedral and 14,904 wedge elements.In addition, for this grid h min = 1.72 cm and h max = 0.55 m.The number of degrees of freedom N DoF is 685,245,422.
In Fig. 6, the computational grid on the aircraft's surface is depicted, along with the grid in two planes: z = −1 m and x = 0 m.
For this numerical experiment we used the supercomputer Mahti.When we tried this example storing the matrices C and D as usual (see Eq. ( 9)), we ran out of memory.An estimate for the memory requirement for solving the problem by keeping all necessary matrices in memory is 68 terabytes.The low memory version relaxes this memory requirement by a factor of 0.2.So we switched to the low memory version described in Section 3.5 to compute the results shown here.In Mahti, we used a total of 200 computing nodes and 100 CPU units per node.The total time for the calculation was 18 hours, including building the system matrix

Resistive sheets 4.2.1 Salisbury screen
We next model a "Salisbury screen" (W.W. Salisbury, U. S. Patent US2599944 A 1952).For simplicity, assume ϵ r = µ r = 1.Our standard incident field (10) propagates normally to a To the left of the resistive sheet the total electric field is where (0, R 2 , 0) T is the polarization of the reflected wave.Between the sheet and the PEC surface, where −H < x < 0, Here (0, q 0,2 , 0) T and (0, q 1,2 , 0) T are the polarizations of the left and right going waves respectively in the gap −H < x < 0.
Imposing the PEC boundary condition at x = 0 and the resistive sheet transmission conditions (11)-( 12) at x = −H shows that .
For given κ, H, and η we can compute total field in each region and compare to the analytic solution.As is well known, one choice of η gives R 2 = 0: A particularly interesting case occurs when cot(κH) = 0 or H = π/(2κ).Recalling that the wavelength of the radiation is λ 0 = 2π/κ, we see that zero reflection occurs when H = λ 0 /4.To test the UWVF for resistive sheets, we use a rectangular parallelepiped computational region with faces normal to the coordinate directions (see Fig. 9).The rightmost face x = 0 is PEC (Q = −1) the left-most face is an ABC with Q = 0 where we use a non-homogeneous absorbing boundary condition to excite the plane wave.There is no need for a PML since the solution is a wave propagating orthogonally to the absorbing boundary.
On the remaining faces Q = ±1 chosen so that the incident plane wave propagates along the box without distortion.We take the radiation to have frequency f = 2 GHz and place the resistive sheet one quarter wavelength (H = 3.75 cm) from the PEC surface.Figure 9 shows a surface of the grid used in the computations.In this case, the grid consists of 136 tetrahedral elements with 51 vertices.In addition, h min = 3.33 cm, h max = 16.18 cm, and N DoF = 9, 946.
Results are shown in Fig. 10.We plot the magnitude of the y-component of the total field as a function of x when y = z = 7.5 cm.Choosing η = 1 the magnitude of the field is flat to the left of the resistive sheet showing that this choice of η gives rise to no reflected wave.However when η = 0.5 the non-constant magnitude of the total field indicates that a reflected wave is present to the left of the sheet.
In this example just two plane waves per element would suffice to compute the solution, but our code is not tuned to this example and chooses the number of directions as if this is a general Maxwell problem via the condition number heuristic in Appendix A.

Sphere with resistive sheet
In this second experiment with resistive sheets, a PEC sphere with a radius of 1 meter is placed at the origin of a cube [−1 − 15λ 0 , 1 + 15λ 0 ] 3 .Surrounding this sphere is a spherical resistive sheet of radius 1 + λ 0 /4 and surrounding both is an artificial sphere of radius 1 + λ 0 .Outside this artificial sphere we compute the scattered field, and inside the total field (see Section 3.2).The incident field then gives rise to a source on the artificial boundary as detailed in Section 3.2.A PML with a thickness of 5λ 0 is applied to the inside each side of the cube, and the frequency of the incident field is set at f = 2 GHz.The incident field gives rise to a source located on the sphere.
To generate the computational grid, we utilized curved elements and a mesh size parameter of h s = 2λ 0 on the PEC and resistive sheet surfaces.The grid comprises of 428 wedge elements (representing the domain between the resistive sheet and PEC sphere) and 7,589 tetrahedral elements that cover the main domain of interest.Furthermore, the surrounding PML layer is discretized using a combination of 104 hexahedral and 936 wedge elements.The entire grid is composed of 2,550 vertices, with h min = 3.75 cm and h max = 1.48 m.A cross-section of the computational grid is shown in Fig. 11.In this case, the number of degrees of freedom N DoF = 4, 599, 234.We used the supercomputer Puhti with a total of 5 computing nodes and 10 CPU units per node to solve the three configurations for different resistive sheet parameters η.It took 18 (η = 0), 17 (η = 0.5), and 17 (η = 1.0) minutes CPU-time respectively to build the system matrices and then reach the solution with BICGstab.The solution was achieved after 194 iterations for η = 0, 169 iterations for η = 0.5, and 172 iterations for η = 1.0.
In Fig. 12, the total field component ℜ(E y ) on the plane z = 0 is shown for three choices of η.We can no longer expect invisibilty since the screen is curved, but it is evident that backscattering is decreasing as η increases.This is seen more clearly in Fig. 13 where we show the RCS in each case.

Heterogeneous models 4.3.1 A dielectric sphere
In this experiment, a penetrable sphere with a radius of 1 meter is centered at the origin inside the cube [−1 − 15λ 0 , 1 + 15λ 0 ] 3 , where λ 0 is the wavelength in vacuum.For the penetrable sphere, we assume ϵ r = 1.5 + 0.5i and µ r = 1, while we select vacuum parameters in other domains.The frequency of the incident field is f = 2 GHz and a PML with thickness of 5λ 0 is used on each side of the cube.An artificial spherical boundary with radius 1 + λ 0 is used to separate a scattered field region outside and a total field region inside this surface.The scattered-total field formulation in Section 3.2 is used to introduce a source on the artificial boundary.
We utilised curved elements and a mesh size parameter of h s = 3λ s , where λ s is a measure of the wavelength in the penetrable sphere computed using κ abs defined in Appendix A. Figure 14 shows cross-sections of the computational grids.To solve this problem, we used the supercomputer Puhti with a total of 7 computing nodes and 40 CPU units per node.More detailed information on the computations, including information on the grids and CPU time, is given in Table 3.
Figure 14: Cross-section of the computational grid for the penetrable sphere case.
In Fig. 15, we show the total field component ℜ(E y ) on the z = 0-plane for the two meshes.As expected (see also Table 3), mesh 2 with curved face elements produces a solution comparable to mesh 1, but faster.The accuracy of the far field pattern is compared to the Mie series solution in Fig. 16, and again demonstrates that mesh 1 and mesh 2 with curved faces give comparable results.

A plasma sphere
In this second experiment for penetrable objects, we assume a penetrable sphere with a radius of 0.25 meter is centred at the origin of a cube [−0.25 − 15λ 0 , 0.25 + 15λ 0 ] 3 .The material in Figure 16: Bistatic RCS for the penetrable sphere at 2 GHz, otherwise the layout is the same caption as for Fig. 5. the sphere is assumed to be a plasma modelled by setting ϵ r = −1.5 + 0.5i and µ r = 1.The frequency of the incident field is f = 2 GHz and a perfectly matched layer with thickness of 5λ 0 is used inside each side of the cube.The source is introduced on an artificial spherical surface of radius 0.25 + λ 0 .We utilise curved elements and a mesh size parameter h s = 3λ s , where λ s denotes the wavelength in the penetrable sphere, on the material discontinuity surface.The grid comprises of 2,959 tetrahedral elements that cover the main domain of interest.Furthermore, the surrounding PML layer is discretized using a combination of 92 hexahedral and 768 wedge elements.The entire grid has 1,304 vertices, with h min = 5.88 cm and h max = 1.30 m.Here, the degrees of freedom number N DoF is 2,234,416.
We used the supercomputer Puhti with a total of 5 computing nodes and 10 CPU units per node to solve the problem.It took 427 seconds from each CPU unit to build the system matrices C and D and then reach the solution with 103 bi-conjugate iterations.Snapshots of the total field are shown in Fig. 17.In Fig. 18 we show the computed RCS which shows remarkable agreement with the Mie series solution.

Conclusion
In this study, we explored and extended the Ultra-Weak Variational Formulation (UWVF) applied to the time-harmonic Maxwell's equations.Our research findings led to important contributions that enhance the efficiency and applicability of the UWVF method for solving electromagnetic wave problems.The paper shows a series of numerical examples validating the effectiveness of the new enhancements.Scattering problems from PEC objects were considered, highlighting the benefits of curved elements, different element types, and the low-memory version of the software.The applicability was further demonstrated through simulations of scattering from an full size aircraft, emphasising its potential for real-world industrial scenarios.This paper has not only extended the capabilities of the UWVF for electromagnetic wave problems but has also provided a comprehensive set of numerical results to underline the practical significance of these advancements.The integration of curved elements, various element shapes, and resistive sheets collectively contribute to the method's robustness and utility, making it a valuable tool for addressing complex electromagnetic problems.
The choice of mesh size for the surface triangulation needs further study in the case of large imaginary part of ϵ r or in regions of high curvature but these issues are beyond the scope of the paper.Related to this, an important direction for further work would be to refine the heuristics for choosing the number and direction of plane waves on each element.This is particularly needed for elements that might have large aspect ratio such as elongated wedges.

A Choice of basis
Because we have used new element types and larger numbers of directions in this paper compared to [13], we need new heuristics for choosing the number of plane wave directions on a particular element.We use the same technique as in [13].Computing on the reference element, in Fig. 19, the number of plane waves N ℓ is plotted as a function of (κ abs h av ) ℓ , when the maximum condition number of the matrix blocks of D ℓ is limited by the tolerances 10 5 , 10 7 , and 10 9 .Here κ abs = ω | √ ϵ r µ r | and the element size parameter h av is defined as the mean distance of the element's vertices from its centroid.This data is fitted by a quadratic polynomial function with the constraint that the polynomial gives at least 4 directions even on the finest grid: Results of this fitting are shown in Table 4.Although only computed for one element shape, these polynomials are used to set the number of directions for any given mesh.Generally a higher tolerance on the condition number results in more directions per element so greater accuracy, but too high a condition number slows BiCGstab unacceptably.

B Comparison to edge elements
There are many ways to solve the scattering example problems in this paper.To provide a comparison to a more familiar method, we now use ParMax and the edge finite element method to compute scattering from a dielectric sphere.The exact Mie series is used to assess accuracy at several frequencies.We have used the open source Netgen/NGSolve finite element package [21].This is a library of highly optimised C++ function supporting multithreading that can be used via a Python interface (called NGSpy).It includes the Netgen mesh generator.We use pth order Nédélec edge elements of full polynomial degree on a tetrahedral mesh with a rectilinear PML and an outer impedance boundary condition.A direct sparse matrix solver supplied with NGSpy was used to solve the linear system.This limited the overall size of the problem we could solve to approximately 1.2 million degrees of freedom.
In order to make the application of the FEM possible, the computational domain described in Section 4.3 is not used here.Instead, the computational domain is the cube  The computations were run on an Intel Xeon Gold 6138 CPU 2.00GHz with 20 cores, each having two threads per core, and 394GB of memory.Since NGSpy ran multithreaded computations, we ran ParMax with MPI and 20 cores.Elapsed time is reported for NGSpy and ParMax, including assembly, matrix factorization, solving, and computation of the farfield pattern.The error reported is the relative percentage error in the RCS compared to the Mie solution.The results are shown in Fig. 20.
Figure 20 (top) shows that the relative L2-error remains comparable between NGSpy and ParMax across the entire frequency spectrum used here.As illustrated in Fig. 20 (middle), the number of degrees of freedom required by NGSpy increases more rapidly than that for ParMax.Likewise, the comparison in Fig. 20 (bottom) reveals that the CPU-time required by NGSpy increases more significantly compared to ParMax.This could perhaps be ameliorated using an auxiliary space preconditioned iterative scheme [12] rather than a direct factorization, but this is not available to us in NGSpy.

Figure 1 :
Figure 1: Geometry and notation for the resistive sheet calculation.The normal ν is outward to Ω − .

Figure 6 :
Figure 6: Surface triangulation for the aircraft and tetrahedral, hexahedral, and wedge elements on two planes.The colorbar shows the number of plane waves per element.

Figure 7 :
Figure 7: Snapshots of the scattered electric field |E s y | for the aircraft model.In the top left panel we show results in the x − y plane at z = 0, top right is in the z − y plane at x = 0 and bottom left is in the x − y plane at z = 0. Clearly a strong shadow region and multiple reflections are evident.

Figure 8 :
Figure 8: Bistatic RCS for the aircraft at 8 GHz.

Figure 9 :
Figure 9: Computational grid for the Salisbury screen, showing smaller elements between the screen and PEC surface, expanding away to the left.

Figure 10 :
Figure 10: Magnitude of the y-component of the total electric field E y as a function x with y = z = 7.5 cm with η = 1 and η = 0.5, together with the analytic solution.The vertical dashed line marks the location of the resistive sheet.To the left of the sheet, no reflected wave is evident when η = 1 whereas a reflected wave is indicated when η = 0.5.The bottom panel shows the absolute value of the difference between the numerical and analytic solutions, with the difference multiplied by a factor of 1,000.

Figure 11 :
Figure 11: Left: 428 wedge elements forming the interior of the resistive sheet outside the PEC surface.Right: A cross-section of the overall computational grid.

Figure 12 :
Figure 12: Snapshots of the total electric field ℜ(E y ) for the three choices of resistive sheet parameter η (left: η = 0, middle: η = 0.5 and right: η = 1).The dashed red line marks the resistive sheet interface and the solid black line the artificial interface used to introduce the incident wave.Backscattering appears less for η = 1 compared to η = 0.5 or the pure PEC sphere.

Figure 13 :
Figure 13: Bistatic RCS for the PEC sphere and resistive sheet example at 2 GHz.The decreased backscattering due to the resitive sheet covering the sphere at η = 1 is clearly seen.

Figure 15 :
Figure 15: Snapshots of the electric field ℜ(E y ) on the z = 0 plane for the penetrable sphere.The solid white line shows the material interface and the solid black line marks the artificial interface used to introduce the incident wave.

For a comparison
of ParMax with edge finite elements in the case of the dielectric sphere, see Appendix B.

Figure 17 :
Figure 17: Snapshots of the total electric field.Top panel: a contour plot of |E y | in the plane z = 0.The solid white line shows the interface between vacuum and the plasma sphere.The solid black line marks the interface used to introduce the incident wave.Bottom panel: a plot of |E y | along the line y = z = 0.

Figure 18 :
Figure 18: Bistatic RCS for the plasma sphere at 2.0 GHz, comparing the computed RCS and Mie result.Top right corner shows the difference between the numerical and Mie series.

Figure 19 :
Figure 19: The number of basis functions N ℓ as a function of (κ abs h av ) ℓ when the basis dimension is chosen by constraining the maximum condition number of D ℓ .
[−1 − (1 + d pml )λ 0 , 1 + (1 + d pml )λ 0 ] 3 .A PML of relative width d pml = 0.3 with FEM and d pml = 1 with ParMax is added within this cube.For the FEM, we requested NGSpy to use the mesh size λ 0 (2p + 1)/(4π) outside the scatterer and increased p gradually from p = 3 to p = 8 as λ 0 decreased to maintain an accuracy of roughly 1% in the computed RCS.Inside the scatterer the mesh size is reduced by a factor of 0.4.

Table 1 :
Notation used in reporting computational results.

Table 2 :
Table giving details for the computational grids for the PEC sphere, and accuracy and timing.See Table 1 for definitions of the reported quantities.
m.The aircraft is 20.5 m or 547 wavelengths long and 17.8108 m or 475 wavelengths wide.A perfectly matched layer with a thickness of 5λ 0 is added to each side of the cuboid computational domain with side lengths(16.3334,16.3334,4.4826)m.

Table 3 :
Details of the computational results for the penetrable sphere.See Table1for definitions of the reported quantities.
mesh id surface N tetra elements N wedge elements N hexa elements N vertices h min (cm) h max (m) N DoF N BiCG iter L2-error (%) CPU-time (s)