EM modelling of arbitrary shaped anisotropic dielectric objects using an efficient 3D leapfrog scheme on unstructured meshes

The standard Yee algorithm is widely used in computational electromagnetics because of its simplicity and divergence free nature. A generalization of the classical Yee scheme to 3D unstructured meshes is adopted, based on the use of a Delaunay primal mesh and its high quality Voronoi dual. This allows the problem of accuracy losses, which are normally associated with the use of the standard Yee scheme and a staircased representation of curved material interfaces, to be circumvented. The 3D dual mesh leapfrog-scheme which is presented has the ability to model both electric and magnetic anisotropic lossy materials. This approach enables the modelling of problems, of current practical interest, involving structured composites and metamaterials.


Introduction
In anisotropic materials, the electromagnetic material parameters, such as permittivity, permeability and conductivity, may vary in the different crystal directions, so that they must be treated as tensors. It is assumed that already 1000 years ago, before the invention of magnetic compasses, Vikings used crystals, a naturally occuring anisotropic material, in Norse sagas referred to as sunstones to navigate on open water on cloudy days. In accordance to researchers these sunstones could have been calcite crystals where anisotropy leads to the phenomenom of birefrigence (crystalline materials with different indices of refraction with different crystallographic directions). Their sunstone came within 1 % of the true location of the sun [1]. Nowadays anisotropic materials offer many new and interesting perspectives in engineering. A thin anisotropic coating may, for example, significantly change the radar cross section of an aircraft. Composites, anisotropic materials with applications initially limited to stealth bombers, satellites and space shuttles become part of our everyday life. Due to their advantages with respect to mechanical strength and weight compared to metals they are now used in civil aircrafts, trains, automobiles, trucks, sports equipment and so on. Especially in plane and cars electromagnetic compatibility is an issue which can be dealt with using numerical simulations. Other applications are the design of patch antennas where anisotropy can be used as design parameter [2]. Furthermore anisotropy is the basis of metamaterials, which are new materials with electromagnetic properties which cannot be found in nature, e.g. a material may have a negative index of refraction, which could be employed in the design of invisible cloaking devices [3].
Analytical solutions to wave propagation problems in electromagnetics are mainly restricted to problems involving simple geometrical shapes and diagonal, uniaxial or biaxial, tensors [4,5]. Numerical techniques are required for the solution of the majority of problems, which involve arbitrary shaped objects. The second order accurate standard Yee algorithm implemented on a pair of staggered orthogonal Cartesian meshes is often the favored computational solution technique because of it's low operation count and its low storage requirements. Unfortunately in the case of a curved boundary where the physical boundary doesn't conform to the orthogonal mesh a very fine mesh is required due to errors induced by the stair stepping edges. In earlier work [6], we demonstrated the capability of a generalized Yee algorithm adapted to unstructured meshes to accurately model the radar cross section (RCS) of arbitrarily shaped lossy dielectric objects. For isotropic cases our method shows significant savings with respect to memory and time with respect to the standard FDTD scheme due to the unstructured mesh we employ. Here, we describe the extension of the method to deal with anisotropic materials, such as composites.
There are generally two methods to deal with anisotropic materials. Firstly you use the constitutive equation to replace the displacement field in Maxwells equations by the electric field [7]. Another possibility is to obtain the displacement field and afterwards use it in the constitutive equations. We adopt the latter method which has been proposed by [2]. This approach was originally presented within the context of a total field formulation, but the unstructured mesh extension adopted here employs a scattered field formulation.

Problem formulation
The integral form of Maxwell's equations is employed [8]. For a three dimensional lossy dielectric medium, Ampère's and Faraday's Laws are expressed, in a scattered field form, as and Here, ∂ A denotes a closed curve bounding a surface A, dA is an element of surface area, directed normal to the surface, dl is an element of contour length, in the direction of the tangent to the curve, t denotes time andĪ is the unit matrix. In addition,ε is the electric permittivity tensor,μ is the magnetic permeability tensor,σ andσ m are the electric and magnetic conductivity tensors respectively and ε 0 and μ 0 denote the electric permittivity and magnetic permeability of free space respectively.The subscripts (.) inc and (.) scat refer to incident and scattered field components, with the total fields regarded as being formed as the sum of the corresponding incident and scattered fields. The vectors D, E, B and H represent the electric flux density, or displacement field, the electric field, the magnetic flux and the magnetic field respectively. The constitutive equations, relating the electric field intensity E to the electric flux density D and the magnetic field intensity H to the magnetic flux density B, may be expressed as The incident field is assumed to be a monochromatic plane wave, generated by a source located in the far field, which has the form E inc = E 0 cos(ωt − k · r), where E 0 is the electric field vector, k is the wave vector, r is the position vector and ω is the angular frequency. From the known incident electric field, the incident magnetic field may be determined, using Faraday's Law, as wherek is the unit wave vector and η = √ μ o /ε o is the impedance of free space.

Discrete equations
As the electric and magnetic field vectors are orthogonal, we employ a primal polyhedral mesh, to store the electric field and displacement field projections. The Voronoi dual mesh is used to store the magnetic field and magnetic flux projections. For illustration, we assume here that the primal mesh consists of tetrahedra, generated by a Delaunay method [9]. By construction, each Voronoi face is a perpendicular bisector of the corresponding Delaunay edge and each Delaunay face is perpendicular to the corresponding Voronoi edge, fulfilling our orthogonality requirements. This means that in an ideal mesh, each Delaunay edge is perpendicular to a surface bounded by a closed loop of Voronoi edges. Similarly, each Voronoi edge is surrounded by a closed loop of Delaunay edges. The Delaunay mesh is assumed to have N D e edges and the Voronoi mesh N V e edges. When the leapfrog scheme is used for time discretization, it will be second order accurate if the unknowns are located at the midpoints of these edges. The unknown at the centre of the ith Delaunay edge corresponds to the projection, (D scat,i , E scat,i ), of the scattered electric field onto the direction of the edge, as illustrated in Fig. 1a. The unknown at the centre of the jth Voronoi edge corresponds to the projection, (B scat, j , H scat, j ), of the scattered magnetic field onto the direction of the edge, as illustrated in Fig. 1b. In the scattered field formulation, the incident field is a known function, while the scattered field is unknown. At the interface boundaries, the material parameters in the Eqs. (1) and (2) are not constant. In this case, we average the values ofε,μ,σ andσ m at a dielectric interface, leading to the valuesε av ,μ av ,σ av andσ m av . The method for determining these averaged values is detailed in Sect. 5. Direct discretization of Ampère's Law and Faraday's Law then leads to the equations and where t denotes the time step, the superscript n denotes an evaluation at time level t = n t, l D i represents the length of the ith Delaunay edge and A V i corresponds to the area of the Voronoi face spanned by the Voronoi edges surrounding Delaunay edge i. Similarly, l V j represents the length of the jth Voronoi edge and A D j corresponds to the area of the Delaunay face spanned by the Delaunay edges surrounding Voronoi edge j. The numbers j i,k , k = 1, . . . , M V i refer to the M V i edges of the Voronoi face corresponding to the ith Delaunay edge, while the numbers i j,k , k = 1, . . . , M D j refer to the M D j edges of the Delaunay face corresponding to the jth Voronoi edge. In addition, F,ê i denotes the dot product (projection) of any field vector F along the i th edge. We will use both F,ê i and F ·ê i to represent the scalar product between a field and a unit edge vector. The projection of the scattered electric field vector onto Delaunay edge i is denoted by E scat,i , while E scat | i denotes the scattered electric field vector at the location of the ith Delaunay edge. Defining the quantitiesā enables us to write Eqs. (5) and (6) as Here, B n+0.5 scat, j and D n+1 scat,i are projections onto the Delaunay and Voronoi edges respectively, whereas the quantities D n scat i , E n+0.5 inc i , and H n−0.5 scat j , H n inc j represent field vectors computed at the centre of the ith Delaunay edge and the jth Voronoi edge respectively. These field values have to be determined from their corresponding stored projections and this calculation, which is not direct, is described in detail in Sect. 6. In contrast to the isotropic case, these equations cannot be updated in one step, as vector-matrix multiplications are involved. These staggered equations are used to advance the solution in a leapfrog manner. The magnetic field is updated over the dual mesh at the half time step, using Eq. (11), and the electric field is updated over the primal mesh at the full time step, using Eq. (12).
The scheme is based upon the projections of the field unknowns at the edge centres. This allows us to use unstructured meshes and to model electromagnetic scattering, even for objects of arbitrary shape. Full details of the mesh generating process that is employed, to ensure Delaunay and Voronoi meshes with the correct properties, may be found elsewhere [8][9][10]. Here, the additional challenge that is faced, is that we have no direct access to the full field vectors B, H D and E at a given location. Nevertheless, we are able to get accurate values for the approximated field vectors, that are needed for the matrix-vector multiplication of the updating equations, and the results are then projected to the corresponding Delaunay edgeê i or Voronoi edgeê j . This updating process is explained in detail in Sect. 8.

Mesh generation
Structured meshes that are employed for wave propagation problems are generally constructed to have a uniform element size, δ, that is related to the characteristic wavelength, λ, of the problem. A value of δ in the range λ/30 to λ/10 is typical for practical applications of the Yee scheme on regular cartesian grids [11]. Consider now the problem of generating a body fitted mesh of this form for use with the co-volume algorithm outlined above. The computational domain surrounding a general scattering obstacle is discretised employing a hybrid mesh, which is generated in four stages. In the first stage, an unstructured triangulation of the surface of the scatterer is produced [12] and the triangulation is then placed inside a hexahedral box. The region inside this box is discretised using a regular cartesian mesh of cubes, of a prescribed edge length δ. Cubes within a prescribed distance of the scatterer, or lying internal to the scatterer, are removed in a second stage, to create a staircase shaped surface that completely encloses the scatterer. In the third stage, a point distribution is specified to completely cover the unmeshed region, with those points from the distribution that lie either inside the scatterer or outside the staircase surface being removed. The fourth stage consists of using the points that remain to generate an unstructured tetrahedral mesh in the region between the surface triangulation of the scatterer and the staircase surface. The problem of fitting the hexahedral and tetrahedral meshes is overcome by placing a pyramidal element on each exposed square face, leading naturally to a consistent mesh. The unstructured mesh is optimised to ensure that both the primal and the dual mesh are of the highest possible quality. The approach adopted is to relax the requirement that a dual edge must be a bisector of the corresponding Delaunay edge. At the same time, the corresponding dual mesh vertex is moved to a point which still ensures orthogonality between the two grids and which lies inside the corresponding primal element. Primal elements with a common circumcentre, and hence a corresponding Voronoï edge length of zero, will automatically be merged during the solution process, creating general polyhedral cells.

Boundary conditions
In scattering problems, the incident wave is assumed to be generated by a source located in the far field and the physical solution domain is infinite in extent. The numerical simulation of the scattering problem is undertaken on a finite computational domain. For example, the computational domain employed for the problem of simulating scattering by an anisotropic dielectric sphere, located in free space, is illustrated in Fig. 2. The infinite real physical domain has been truncated and, at the truncated outer computational boundary, the scattered field should only consist of outgoing waves only. The modelling of this requirement is achieved by adding a wave damping perfectly matched layer (PML) [13] to the truncated exterior far-field boundary. In earlier work [14,15], we have shown that it is convenient to represent the PML region with an assembly of regular hexahedral computational cells.

PEC boundary conditions
In a scattered field formulation, the condition is applied at the surface of a perfect electric conductor (PEC).
Here, n is the unit outward normal vector to the PEC surface. Without changing Eqs. (11) and (12), we can strongly impose the electric field unknowns, at the set of edges forming the PEC interface, to satisfy the condition of Eq. (13). Within this leapfrog scheme, we can also model thin resistive or PEC sheets, by assigning the sheet conductivity only to the Delaunay edges forming the interface.

Material interface boundary conditions
When the boundary is an interface between two different media, the update Eqs. (1) and (2) require integration across the interface. These integrals are evaluated by assigning a weighted average value to the material parameters, based upon the mesh structure.
In a previous publication [6], we obtained material parameters at the interface by a weighted arithmetic mean average and compared the results produced with those obtained by using the arithmetic mean average, the harmonic mean average and the geometric mean average. We demonstrated that the weighted arithmetic mean average resulted in improved accuracy on unstructured meshes. Here, we adopt the same form of averaging, but applied now to every component of the material parameter tensors. In the isotropic case, the scalar material properties, ε and σ , associated to the electric field, are stored on the Delaunay edges and the scalars μ and σ m associated to the magnetic field are stored on the Voronoi edges. For an anisotropic material, the scalars ε, σ , μ and σ m become second order tensorsε,μ,σ andσ m and, for each component of these tensors, we employ the averaging used in the isotropic case. In Fig. 3a, b, the colour code employed in the diagram of Fig. 2 is adopted, with blue indicating Delaunay edges in medium (1), green indicating Delaunay edges in medium (2) and red the Delaunay edges forming the interface. The Voronoi edges, surrounding a given Delaunay edge at the interface, are indicated in black. In Eqs. (2) and (1), the quantitiesε,μ,σ andσ m lie inside the integrals. In the discretized Eqs. (12) and (11), each component is averaged over the closed surface loop and the averaged values can be taken outside the integrals. For these Delaunay and Voronoi edges, each component (q, l) of the material property tensors is evaluated, using weighted average formulae, as Here q and l can take the values 1, 2 or 3, corresponding to the x, y or z directions respectively, while M V i refers to the number of Voronoi edges surrounding a given Delaunay Edge i. As there are two sub-volumes associated to each Voronoi edge, we have to sum over 2M V i Voronoi edges. These account for the contribution of the material parameter assigned to each of the cells surrounding Delaunay edge (.) Del,i , weighted by a coefficient w k that corresponds to the volume spanned by the two endpoints of the Delaunay edge, the intersection point of the Voronoi edge with the Delaunay face and the position of the circumcentre of the cell. For example, to obtain ε av q,l Del,i , in Fig. 3a, w 1 would be the volume spanned by the points N 1, N 2, C1 and Fi1, while ε q,l Cell,1 would correspond to the permittivity component of the element Cell1. The points labelled N belong to the Delaunay mesh, the points labelled C are the circumcentres of the corresponding cells and Fi refers to the intersection point of the Voronoi edge with a face spanned by Delaunay edges.
Each component of the tensors of magnetic permeabilitȳ μ and magnetic conductivityσ m , which are linked to the Voronoi edges, is obtained, by averaging, as μ av q,l V or, j = g 1 μ q,l Cell1 + g 2 μ q,l Cell2 g 1 + g 2 σ m av q,l V or, j = g 1 σ m q,l Cell1 + g 2 σ m q,l Cell2 The lengths of the Voronoi edges, g1 and g2, inside cell1 and cell2 are the distances between the intersection point Fi1 of the Voronoi edge (.) V or, j with the Delaunay face and the circumcentre of the cell. For example, in Fig. 3b, the distance between the points C1 and Fi1 is g 1 .

Obtaining approximated field vectors from edge projections
In Sect. 3 it was noted that the main difficulty with Eqs. (11) and (12) are the matrix-vector multiplications. This is because we only store incomplete field vectors, due to our projection based nature of the updating scheme. The challenge is now to obtain the corresponding field vectors D scat and B scat and associate them with Delaunay edge i and Voronoi edge j respectively. Unfortunately, we cannot get exact full field vector components from field to edge projections. However, we can approximate the full field components at any location in the mesh. To achieve this, we assume that, in R 3 , with a set of three orthogonal vectors v 1 , v 2 , v 3 , a general vector x can be reconstructed as in terms of the projections As the mesh is assumed to be unstructured, we cannot use this method directly to obtain an exact field vector. The trivial solution would be to consider all the edges connected to one node, as depicted in Fig. 4b, and solve a system of equations. This will give us a field vector at node n 1 . As one edge is always formed by two points, we can do the same for point n 2 . Finally, we average these field vectors and project to the corresponding edge. For example, to obtain the displacement field vector D scat on the edge i, denoted by the red line in Fig. 4a, we have to construct the set of equations andê i is the normalized Delaunay edge vector corresponding to Delaunay edge i. The matrixP, based upon the x, y, z components of the vectorê i and the projections D scat ·ê i are known. The displacement field vector at a node belonging to Delaunay edge i is approximated by considering the sum of the system in Eq. (18) for each Delaunay edge connected to that node, as depicted in Fig. 4b. In this case, we solve the system Here, D scat is the unknown vector,P and D scat,q ·ê q are known, N is the number of Delaunay edges connected to the node and e i 1 = e i x , e i 2 = e i y and e i 3 = e i z . We solve this system of equations locally, node by node, until we have an approximated field vector at all nodes of the Delaunay mesh. Finally, we have to link the computed field vector to the corresponding edges. A Delaunay edge is assumed to connect the two nodes n 1 and n 2 and the displacement field vector, associated to the Delaunay edge i, is obtained by averaging the electric field vectors at nodes n 1 and n 2 , i.e. D scat,i = D scat,n 1 + D scat,n 2 /2. The same procedure is applied for approximating the magnetic flux vectors B scat , but using now the Voronoi edges. However, unlike the Cartesian Yee scheme, where we have one update equation for each component of the field vectors, the approximations of the full vector fields obtained using this method are not good enough for time iterating the field components. It is found that error accumulation causes the algorithm to become unstable. In the following sections, we suggest how this difficulty may be circumvented.

Local orthogonal unit vectors
When using the integral formulation of Maxwell's equations, it is not the displacement field vector D scat that is updated, but rather its projection, D scat,i = D scat · e i , onto a Delaunay edge e i . In the case of an isotropic material, the electric permittivity and the magnetic permeability are scalars, so that updating the fields only involves scalar multiplication between the field projections and the scalar material properties. For an anisotropic material, the integrals in Ampère's and Faraday's laws contain matrix-vector multiplications between material tensors (ε,μ,σ,σ m ) and the fields (D scat , B scat ). To deal with these matrix-vector multiplications, we first create two linearly independent vectors for each Delaunay and Voronoi edge. Using the stabilized Gram-Schmidt orthonormalization procedure, we generate three orthonormal vectors, as illustrated in Fig. 5. The first vector e 1 , to which the two linearly independent vectors (e 2 , (e 3 ) are added, remains unchanged during the whole process. Each set of three orthogonal vectors represents one local coordinate system, leading to as many local systems as Delaunay and Voronoi edges. Due to the discretisation, we can reconstruct approximated field vector components for each local coordinate system, using the field projections of the surrounding Delaunay or Voronoi edges, as explained in Sect. 6. The field vectors, obtained in this way, are projected onto the two orthogonal vectors forming the local frame. The first vector of each subset, which remains unchanged during the orthonormalization procedure, can be immediately updated using the projection equation employed in the isotropic case. For this projection, no error is induced by the field averaging.

Coordinate transformation
The material tensors are expressed in the global reference frame formed by the orthonormal vectorsx,ŷ andẑ. Three orthonormalized vectorsx ,ŷ ,ẑ form the basis of each local coordinate system. The (.) always refers to unit vectors and (. ) to vectors or vector components of the local coordinate system. The Jacobian matrixJ , defined by can be used to transform a vector, or a matrix, from the global to a local frame. Here (x , y , z ) refers to the coordinate in the local coordinate system and (x, y, z) to the global coordinate system. Each component of the Jacobian can be interpreted as an amplification factor, describing how one coordinate in a given reference frame stretches, shrinks or rotates with respect to another coordinate in another reference frame. In our case, the Jacobian is pure a rotation matrixJ R which can be directly calculated as Maxwell's equations are form invariant [16,17], which means that, in the local coordinate system, they may be expressed as The electric and magnetic fields, in local and global coordinates, are related and a material parameter tensor,M say, is transformed into the local coordinate system using the operator transformation Due to the form invariance, it follows that we can absorb the coordinate transformation completely into the material properties. Note that the determinant of the rotation matrix is unity, i.e. det(J R ) = 1.

Time updating scheme
As mentioned previously, Eqs. (11) and (12) cannot be updated simultaneously, as the scattered field vectors D scat and B scat are not immediately available. In this section, the updating process will be explained in detail. For simplicity, we will restrict consideration to updating the electric field projection E n scat , as the magnetic field H n scat will be similarly updated. We first apply the coordinate transformation to the termsā ε+ ,ā ε− defined in (7) andε −1 av from Eq. (12), so that These are stored at the corresponding Delaunay edges before entering the time loop. Within the time iteration loop, we first calculate and store the right hand side of the Eq. (10) for each Delaunay edgeê i . This is readily accomplished, as the magnetic components in the circulation term are available, from the previous iteration, and the full vector E n+0.5 inc,i is also a known function. Before updating the projection D n scat to D n+1 scat , Eq. (20) is employed, to obtain the vectors D n scat and Z D as defined above . These vectors are projected in each of the three orthogonal directions of every local frame, e.g.ê 1 ≡ê 1 ,ê 2 ,ê 3 from Fig. 5, and the matrix vector multiplication is performed. Equation (12) in the local frame takes the vector form We can use the constitutive equation to obtain the electric field projection. In practice, we only have to consider the first line ofε

Validation
A series of examples, involving scattering of an incident plane wave by an anisotropic sphere, is included here to demonstrate the numerical performance of the algorithm that has been described. The algorithm is validated by comparing the results produced with those obtained from the open source program Discrete Dipole Scattering (DDSCAT) [18], which is an implementation of the frequency domain discrete dipole approximation [19]. The incident wave has free space wavelength λ 0 = 1m and it propagates in the x direction.
In each case, the electrical length of the sphere is 2λ 0 . The mesh employed is illustrated in Fig. 6 and has an average edge length of λ 0 /20. For each example, the PML region is located at a minimum distance of λ 0 from the scatterer and the PML is discretised using 10 layers of hexahedra. The minimum distance between the inner boundary of the PML and the surface of the scatterer is represented by 8 cells. The complete mesh consists of 876, 116 cells, 1, 673, 527 Delaunay edges and 2, 076, 019 Voronoi edges. The distribution of the radar cross section of the cross-and co-polarized scattered waves is computed as where the phasor amplitudesȆ andH , expressed in spherical coordinate system (r, θ, φ), are calculated by taking the Fourier transform of the time domain solution. The quantity displayed in each case is the radar cross section RC S(θ, φ) = 10 log 10 (χ ) Due to the spherical symmetry of these examples, we only display the RCS values in the range 0 • to 180 • . Further details of the RCS computation and the nature of the mesh can be found elsewhere [6].
It is worth noting that our unstructured mesh scheme only needs 8 points per wavelength to obtain accurate results in free space [20]. However, inside a dielectric, the wavelength λ Diel is less than the wavelength in free space λ 0 . For the anisotropic case, it is estimated that a mesh spacing of λ 0 /20 mesh should be sufficient. As we include a test case with full anisotropic tensors for both electric and magnetic properties, we decided to use the same mesh for all the test cases.

Magnetically uniaxial non-lossy anisotropic sphere
The first test case involves an uniaxial permeability tensor and is devised to check the updating of the magnetic field projections. For the sphere, the material parameters are used. In this case the only term in Eq. (5) that involves matrix multiplication is (μ av − μ 0Ī ) ∂ ∂t H n inc j . The matricesā μ+ andā μ− reduce to the unit matrix and the update of of magnetic field from the constitutive Eq. (31) only requires multiplication by the coefficient (μ −1 av ) 11 . Figure 7 shows the computed distributions of both the cross-polarised (σ θφ ) and the co-polarised (σ θθ ) RCS compared with the distributions obtained from using the discrete dipole approximation (DDA). It can be seen that the RCS distributions are in excellent agreement, apart from the differences in the troughs, which are typical of comparisons between timedomain and frequency-domain approximations. In this case the only term that require matrix multiplication in Eq.
, as the other matrices a ε+ andā ε− reduce to the unit matrix. The update of the electric field from the constitutive Eq. (31) only involves multiplication by the coefficient (ε −1 av ) 11 . Figure 8 shows good agreement between the computed RCS distributions and those produced with the DDA method. There seems to be a big difference between our solution and the solution from DDA for the co-polarized RCS distribution in Fig. 8a but this is mainly due to the logarithmic scaling we are using highlighting even small differences between the results.

Electrically lossy anisotropic sphere
The next example involves an anisotropic permittivity tensor, together with electric conductivity. The material parameters are employed. In this case, the matrix multiplications are required in Eq. (6) and the time updating scheme described in Sect. 8 was used. The RCS distributions computed with the 3D-leapfrog and with the DDA scheme are compared in Fig. 9. Steady state conditions were attained in the time domain approach after ten cycles of the incident wave.
In this example, matrix multiplication is required in Eq. (5) and the scheme described in Sect. 8 was used for the magnetic field update. The time domain solver reached steady state after ten cycles of the incident wave. The comparison between the RCS distributions computed with the current time domain approach and frequency domain method is given in Fig. 10.

Computational cost of the anisotropic model
The next example includes the use of full anisotropic tensors for both electric and magnetic properties shown in Fig. 11.
Comparison with DDA results is not possible in this case, as the DDA code only allows for the modelling of non-magnetic materials. To investigate the time penalty that results from the use of the anisotropic model, two simulations were performed on the same mesh. In the first example, the sphere was taken to be an isotropic lossy dielectric material, with μ > 1, ε > 1, σ > 0, σ m > 0. In the second example, the sphere was modelled as an anisotropic lossy dielectric material, so that μ = μ rĪ ,ε = ε rĪ ,σ > [0] 3×3 ,σ m > [0] 3×3 . It was found Fig. 11 Scattering of a plane wave by a fully anisotropic dielectric sphere of electrical length 2λ: a contours of H z shown on a cross-section through the mesh; b contours of E y shown on a cross section through the mesh that the computational cost for the example involving the anisotropic sphere was ten times the cost of the solution for the isotropic sphere. This extra cost mainly arises from equation (20), which implies a requirement to solve a system of three equations for each Voronoi and Delaunay node. Considering these additional costs, is this scheme competitive when compared to the standard FDTD method? In our first paper [6], we showed that the accuracy of this scheme, with a λ/15 unstructured mesh, compares favourable with that of the standard FDTD method, with a λ/90, λ/120 structured mesh, for objects of curved shape. This means that we get comparable results when using a mesh that is 6 to 8 times coarser. A standard FTDT Yee's cell has 12 edges and 6 faces where the electric and magnetic field components are stored respectively. This corresponds to 18 degrees of freedom. Discretising one Yee's cell requires 6 tetrahedra, in the worst case scenario, leading to 19 edges and 14 faces. As we store the electric field projections at the Delaunay edges and the magnetic field projections at Voronoi edges intersecting with the faces, we then have 33 degrees of freedom in total. Although the number of degrees of freedom have nearly doubled, this is more than compensated by the fact that we use a mesh that is, at least, 6 times coarser. Furthermore, because we are working in three dimensions, this should also be taken into account. This implies that, for the same volume, an unstructured mesh with 33 degrees of freedom produces the same accuracy as a structured mesh with 6 3 × 6 = 1296 degrees of freedom. The extra operations, needed for averaging and reconstructing the vectors described in Sect. 6, imply a cost penalty. For each node, we require 18 averaging operations for each field, in addition to 3 × 3 operations of vector reconstruction for each edge. For the 6 tetrahedra lying inside one cube, the number of operations will be 2×8×18+9×(19+14) = 585 for all degrees of freedom, taking into account the 19 electric field vectors and 14 magnetic field vectors. For one cell of a Cartesian mesh, 8 operations are required for averaging the two offset components of a field vector, leading to 48 extra operations on one cell node or 144 per Yee's cell. However, in our case dealing with curved boundaries, to achieve the same level of accuracy, the standard scheme will actually cost 6 3 × 48 = 10368 operations per same computational volume whereas we only have 585 operations.

Transmission of a narrow band pulse
As final example, we evaluate the transmission efficiency of a pulse through a radome made of anisotropic dielectric containing conducting fibres. Radomes of this type are typically used to conceal aircraft communications or radar systems. The radome consists of half an ellipsoid, with a lateral radius of 0.5m, a length of 1m and a thickness of 0.05 m. Figures 12  and 14 display a radome mesh with 20 points per wavelength. This mesh is only used for visualization, as it is too coarse for a 1GHz pulse. For the actual simulation, a finer mesh with 40 points per wavelength is used. In [6], we showed that the number of element layers employed to represent a given thickness of material has no impact on the results, provided that we meet the minimum number of points per wavelength requirement.
The incident plane wave used to illuminate the Radome is a narrow band pulse where τ denotes the pulse width, k is the the wave vector and r is the general position vector. Composites are more complicated to model than standard anisotropic materials, due to the inner structure, e.g orientation of the fibres. In Fig. 13a we illustrate a cut through a composite slab with the fibres oriented in theŷ direction. Due to the specific orientation of the slab, the material frame (x ,ŷ ,ẑ ) and the global frame (x,ŷ,ẑ) are identical. In the case of the curved shell of a radome, the situation is more complicated, as the orientation of the fibres changes in space. In this case, for each location on the shell, we have a material frame which, in general, differs from the global frame as illustrated in Fig. 13b.
In the previous sections, we have only used the coordinate transformation to pass from a global frame (x,ŷ,ẑ) to a local frame linked to each Voronoi or Delaunay edge (x ,ŷ ,ẑ ).
For a composite, we first have to make a coordinate transformation from the material frame, linked to the orientation of the fibres, (x ,ŷ ,ẑ ) to the global frame (x,ŷ,ẑ). Then we pass from the global frame to a local frame linked to each Voronoi or Delaunay edge, according to (x ,ŷ ,ẑ ).
Here,J Ri , i = 1, 2, are the two transformation matrices, whereJ R2 is identical toJ R from Eq. (23). First, we have to create an orthonormal material frame. As in the preceding sections, an orthonormal system simplifies the coordinate transformation because the transformation matrix becomes a simple rotation matrix. The procedure here is illustrated in Fig. 14a. Initially, we only consider Voronoi edges at the dielectric interface. Each Voronoi edge (Vor1,blue) crossing the face of a tetrahedron is surrounded by three Delaunay edges. We select two of these edges, (Del 1 ,Del 2 ). By construction, they are parallel to the surface. Using the cross product, we create a vector x = Del 1 × Del 2 perpendicular to the dielectric surface. Finally, we take the cross product to obtain the vector y = z × Del 1 and assign z = Del 1 . After normalizing, we have one orthonormal material frame (x ,ŷ ,ẑ ) for each face of the dielectric interface. The next step consists in linking the material frame to the cells inside the composite. To achieve this, we compare the distance between the intersection points of a Voronoi edge with the interface and the circumcentre of all cells C 1 inside the dielectric. After finding the smallest distance, we link the local material frame (x ,ŷ ,ẑ ) from the surface to the corresponding cell. This process is illustrated in Fig. 14b.
x · x x · y x · z y · x y · y y · z z · x z · y z · z ⎤ ⎦ J R2 = ⎡ ⎣ x · x x · y x · z y · x y · y y · z z · x z · y z · z ⎤ ⎦ (40) Fig. 13 Composite material: a material frame b Distance between cell centres and material frame located at surface The second coordinate transformation converts the parameters from the global to the local frame and this process has already been described in the preceding sections. If we consider, for example, Eqs. (7,8), the material parameters ε,σ,μ,σ m have to be replaced byε ,σ ,μ ,σ m . This is true for each equations in which the material parameters appear. Our radome is characterized by the material parameters Due to the orientation of the fibres, the properties of our composite are the same in them and thel directions, but differ along then axis. For our material parameters, this means thatM(1, 1) =M(2, 2) =M (3,3) , and the other components are 0, forε ,σ andμ =Ī andσ m = 0. Typically, a material that minimally attenuates the electromagnetic signal transmitted or received by the antenna is used. The transmission is evaluated as which corresponds to the ratio of the amplitude of the total electric field divided by the amplitude of the incident electric field at a point r 0 inside the radome. The transmission of 1 GHz narrowband pulse through a composite radome is represented in Fig. 15.

Conclusion
A Yee type algorithm has been implemented, on an appropriately generated unstructured mesh, to model electromagnetic wave scattering by bodies consisting of electrically and magnetically anisotropic and conducting dielectric materials. The implementation has been successfully validated by comparison with the results obtained using the discrete dipole approximation. The generalization from isotropic to anisotropic materials allows us now to accurately model anisotropic objects of complex geometry.
In future work, we expect to extend the method to allow for the modelling of anisotropic dispersive materials. Bi-anisotropic dispersive materials, such as metamaterials, which involve coupling electric and magnetic fields in the constitutive equations could also be included. This can be achieved using the Z-Transform [21][22][23]. It is also proposed to incorporate a multi-scaling procedure, to allow for the modelling of more complex materials, such as composites.