Analytical ﬁnite element matrix elements and global matrix assembly for hierarchical 3-D vector basis functions within the hybrid ﬁnite element boundary integral method

. A hybrid higher-order ﬁnite element boundary integral (FE-BI) technique is discussed where the higher-order FE matrix elements are computed by a fully analytical procedure and where the gobal matrix assembly is organized by a self-identifying procedure of the local to global transformation. This assembly procedure applys to both, the FE part as well as the BI part of the algorithm. The geometry is meshed into three-dimensional tetrahedra as ﬁnite elements and nearly orthogonal hierarchical basis functions are employed. The boundary conditions are implemented in a strong sense such that the boundary values of the volume basis functions are directly utilized within the BI, either for the tangential electric and magnetic ﬁelds or for the asssociated equivalent surface current densities by applying a cross product with the unit surface normals. The self-identiﬁed method for the global matrix assembly automatically discerns the global order of the basis functions for generating the matrix elements. Higher order basis functions do need more unknowns for each single FE, however, fewer FEs are needed to achieve the same satisﬁable accuracy. This improvement provides a lot


Introduction
The Finite Element Boundary Integral (FE-BI) method (Jin, 2002;Tzoulis and Eibert, 2005;Eibert and Hansen, 1997) is an efficient numerical technique for solving electromag-netic field problems.Traditional finite element methods rely on utilizing the local information of the FEs.The fixed local node order forces the local matrix elements to be transformed into global ones.Facing low order (LO) basis functions, the local-global transformation is easy as edge related elements only follow the edge directions.When it comes to higher order (HO) basis functions (Jin, 2002;Sun et al., 2001;Ismatullah and Eibert, 2009;Jorgensen et al., 2004;Graglia et al., 1997;Jorgensen et al., 2005), the basis functions are also related to faces, volumes or even more complicated structures.The local-global transformation procedure introduces then considerably more difficulties.In this paper, a self-identified hierarchical basis function method is illustrated.This method effectively overcomes the problem mentioned above and provides more feasibility within FE-BI.Without fixing the node order or the sequence order of the basis functions for the local FEs, the self-identified hierarchical basis function organization allows a simple assembly of the global equation system.Simultaneously, this method guarantees the compatibility between FE and BI (Ismatullah and Eibert, 2009) fluently.All arbitrarily shaped components are meshed into tetrahedra (Volakis et al., 1998) apart from perfect electric conductors (PEC) or perfect magnetic conductors (PMC) where E and H are forced to vanish inside the volume.As FE-BI solves for the field distribution inside the volume together with the corresponding equivalent surface currents (Ismatullah and Eibert, 2009;Rao et al., 1982), the self-identified hierarchical basis functions describe the distribution of fields within the tetrahedra.When observation points tend to the enclosed boundary surface, the boundary condition determines the continuity of E and H fields (Harrington, 1961;Bladel, 1964;Mautz and Harrington, 1978).So it is useful to guarantee that the basis functions for FE are the same as the basis functions for the BI.
Derived from the equivalence theorem, basis functions for BI are used to compose equivalent 2-D surface currents (Rao et al., 1982;Chew et al., 2001;Ylä-Oijala and Taskinen, 2003).The currents are relevant to the surface unit normal vector and the polarization of the fields, thus the tangential component of FE basis functions on the surface is perpendicular to the current basis functions.With respect to the corresponding sources (electric current J s or magnetic current M s ) and the surface normal unit vector, the corresponding subspaces ensure the compatibility between FE and BI.
The LO basis functions have shortcomings when the simulation accuracy and the large number of unknowns are considered.Precision and efficiency of LO are difficult to improve with increasing numbers of unknowns.The solution with LO basis functions demands the mesh size to be around λ/20 to λ/8.With coarser meshes, the elements may introduce inaccurate waveforms for field reconstruction.
The well-known mixed order basis functions are successful for electromagnetic field distributions and surface current reconstructions.Rao-Wilton-Glisson (RWG) (Rao et al., 1982) basis functions are inherited as LO.As RWG is very effective for BI, it is implemented as the first order of Rotational Subspace (0th order) in FE tetrahedra.The Nedelec HO basis functions also form the first order of Gradient Subspace (1st order), the second order Rotational Subspace (2nd order), an so on, where this paper is restricted to 2nd order.Apart from BI, LO and HO basis functions are also utilized within FE and they improve the accuracy for field computations.In FE-BI, 0th and 1st order basis functions for FE and BI are easy to match.Both of them are edge related and follow the same direction of the edge vector.The situations for 2nd order basis functions are more complicated.2nd order basis functions are face related, so that to achieve compatibility between FE and BI, the basis functions of them have to maintain the same global node order.The tetrahedral FE basis functions defined in Table 1 have a format represented by subscripts k = (i, j ) and k = (r, s, t) as illustrated in Fig. 1.As elements of matrices need row and column positions, the subscripts (m i , m j ) and (m r , m s , m t ) are introduced for row basis functions and subscripts (n i , n j ) and (n r , n s , n t ) are assigned to column basis functions.In this work, (m i , m j ) and (n i , n j ) contain the global order of local node numbers for edge related basis functions, whereas (m r , m s , m t ) and (n r , n s , n t ) represent the global order of nodes for face related basis functions, where k = (i, j ) and k = (r, s, t) store the local node number of finite elements.
In FE analysis, the self-identified hierarchical basis functions are derived from the geometrical information of the tetrahedra.In the mesh file, the data structure for each tetrahedra contains six edge identities and four face identities.Each edge is constructed with two node numbers and the order of the two nodes determines the edge global direction.Every face has the identity of the corresponding outside boundary triangles, inside volume triangles or inside boundary triangles.The outside boundary triangles are described by three edges in certain directions.The first edge gives out the first two nodes of the outside boundary triangle, the last node can be found through another two edges.The order of these three nodes are inherited as global node order.The inside volume triangles and inside boundary triangles are constructed by three nodes directly and the node order is viewed as the global node order.Through the index of edge and face identities, the tetrahedral FEs can easily consult the corresponding edges and triangles.Thus the determined global node order can be set for the LO and HO basis functions.As shown in Fig. 1, (m i , n i ) always represent the starting point of the edge, (m j , n j ) represent the ending point and (m r , n r ), (m s , n s ), (m t , n t ) represent the node order of the triangle.Practically, the local node numbers are arrayed in the unique global order and assigned to the corresponding subscripts.When generating the system matrices, the assigned subscripts are set to the corresponding positions into the list of basis functions.Thus elements of system matrices are automatically assigned to global edges and faces.If HO basis functions are implemented into FE, the matrix elements can be calculated analytically and precisely.Since these results are commonly not available in FE literature, it is a major contribution of this paper to present these analytical matrix elements up to 2nd order.As the order of basis functions is enlarged, the accuracy of the boundary integral (BI) should also be improved.As it turns out, the integration order for the testing surface integrals should be increased, the adaptive numbers of quadrature points in the singularity cancellation technique have to grow and larger maximum numbers for spherical harmonic expansion terms are needed within the Multilevel Fast Multipole Method (Eibert, 2005).
HO achieves satisfactory accuracy with larger mesh size and it provides a better solution for non uniform finite elements.Good radar cross section (RCS) results of PEC structures coated by dielectric materials are acquired.Since LO is inherited by the hierarchical bases, orthogonality or nearorthogonality of basis functions is usefull for HO.Based on the structure of tetrahedra, system matrices built from nearly orthogonal basis functions converge faster and the solution is more accurate, so that the flexibility of mesh size provided by HO gives a more feasible solution.Meanwhile, HO improves the accuracy and also reduces the number of finite elements.
FE-BI solutions for coated spheres and the Flamme aircraft referring to the self-identified hierarchical nearly orthogonal basis functions are explicitly illustrated.The material of the layered sphere is homogeneous, isotropic and lossy.A variety of FE-BI simulations up to 3 million unknowns based on self-identified hierarchical basis functions are presented.The accuracy of HO testing cases is good, and the simulation results based on HO basis functions are also compared with LO situations.

Finite element variational formulation
Consider the configuration of an arbitrary volume as shown in Fig. 2. A typical FE model consists of a finite volume V a with possibly anisotropic and inhomogeneous materials inside.The materials are characterized by relative permittivity ↔ r (r) and relative permeability ↔ µ r (r).A d is the assembled enclosed envelope and n(r) is the normal vector pointing out of V a .
Assuming a field expansion E = u n α n and H = i n α n as well as a suppressed time dependence e j ωt , a linear system of equations (Jin, 2002;Tzoulis and Eibert, 2005) is achieved as (1) The system matrices [R mn ], [S mn ], [T mn ] and the right hand side vector [w m ] are defined as where a m and a n are field basis functions, k 0 = ω √ 0 µ 0 is the free space wave number, Z 0 = √ µ 0 / 0 is the free space intrinsic impedance, and J d is a volume current density source.
The hierarchical basis functions are based on the normalized barycentric (simplex) coordinates in FE tetrahedra (λ 0 , λ 1 , λ 2 , λ 3 ), where λ 0 + λ 1 + λ 2 + λ 3 = 1.The functions are also related to the gradients ∇λ i (i = 0, 1, 2, 3) and the volume of the tetrahedron V T .Equation ( 2) is related to the curl of the basis functions.Equations ( 2) and (3) form the system matrices related to electric field unknowns [u n ].Equation ( 4) generates the matrices related to surface magnetic field unknowns [i n ] where the integral region is the envelope of V a .Equation ( 5) is related to the current excitation inside the volume.As V a = V e and A d = A e , integration intervals of Eqs. ( 2), ( 3), ( 4) and ( 5) are the combination of finite elements.With simplex coordinates, the properties of basis functions allow for an analytical solution of the matrices in Eq. ( 1).

3-D HO self-identified basis functions
The general format of FE-BI basis functions for different orders is displayed in Table 1.A i (i = 0, 1, 2, 3) is the local face area vector related to surfaces of the tetrahedron pointing into the volume, V T is the tetrahedral volume.Face vectors are independent in the tetrahedron and the surface tangential part of field basis functions compared with surface current basis functions are perpendicular.The field basis functions distribute inside the finite element volume.When tending to the tetrahedral boundary edges and faces, the corresponding tangential parts of edge basis functions for 0th and 1st order exist on the adjacent connected surfaces, they are zero on other edges or faces.The tangential components of face based functions for 2nd order only exist on the related face, but vanish on the other unrelated faces and on all edges.
The compatibility between FE and BI is easily achieved by utilizing self-identified hierarchical bases defined by global node order.The identified hierarchical basis functions α i are composed of rotational first order (0th order) represented by a 1 , gradient first order (1st order) represented by b 1 and rotational second order (2nd order) related to c 2 and d With self-identified hierarchical basis functions, analytical solutions for system matrices can be achieved through the properties of simplex coordinates in tetrahedra.The inner products of simplex coordinates (Lapidus and Pinder, 1982) are given as The dot and cross products between face vectors are defined as Combining identified basis functions and their properties (Table 1, Eqs. 6-10), the system matrix [R mn ] is calculated in Table 2 and [S mn ] is given in Table 3. From Eqs. ( 2) and (3), matrices [R mn ] and [S mn ] are symmetric, so that the elements at symmetric positions in Tables 2 and 3 are identical and it is explicit that the functions are the same along columns and rows.
System matrices [T mn ] depend on the surface boundaries of the volume, so that the dimension of the hierarchical space is reduced.In fact Eq. ( 4) can be written as and the hierarchical space turns out to be With properties of surface triangle finite elements (Lapidus and Pinder, 1982), it is shown that Based on Eq. ( 11), [T mn ] is computed as From Eq. ( 5), [w m ] can be written as where the current source J d can be used directly for the product with self-identified 3-D basis functions and vector [w m ] can be easily solved.Thus, with analytical solutions for the hierarchical basis functions, the system matrices of FE integrals can be generated analytically and this avoids any numerical error accumulations and generates accurate matrices.(Sun et al., 2001).

Surface integral formulation
With Huygen's principle (Ismatullah and Eibert, 2009;Harrington, 1961;Chew et al., 2001), the radiation sources can be replaced by equivalent electric surface currents J s and magnetic surface currents M s on the volume enclosed boundary A d for evaluation of free space radiation.The numerical solution for the surface currents are based on the well known EFIE, MFIE and CFIE (Ismatullah and Eibert, 2009;Rao et al., 1982;Ylä-Oijala and Taskinen, 2003), given as MF I E : where the vector operators and K are calculated as  are 309 644.The running time was 4 763.4 s.For 2nd order, the mesh size is set 0.1 m, and the mean edge length is 10.93 cm, with minimum edge length 6.48 cm and maximum edge length 18.08 cm.The total unknowns are 70 448.The running time was 419.6 s.
The mesh size of LO is around λ/8 and the mesh size of HO is enlarged up to around λ/3. HO with coarser mesh and fewer unknowns achieves also accurate result as LO with finer mesh.

Coated Sphere II
For the second coated sphere testing case, numerical RCS from different orders are also compared with MIE scattering.

415
The second coated sphere contains a PEC sphere core with radius 0.5 m, and the PEC core is enclosed with a dielectric layer with thickness 0.0025 m.The properties of the dielectric layer are presented with ǫ r = 2.5 − 0.5j and µ r = 1.0.The incident wave is 3 GHz and propagating toward +z di-  ( E inc and H inc are incident electric and magnetic fields, G 0 is the scalar Green's function in free space and α is a number between zero and one.J s and M s are given as Thus the electric and magnetic currents can be written as With Ylä-Oijala and Taskinen, 2003), J s and M s are expanded as where f n is the 2-D surface hierarchical basis function for BI, i n and u n in Eq. ( 27) are coefficients with respect to the surface elements and the total number of BI unknowns is The discretized solution for MoM can be achieved through Galerkin's process and fast solutions as given by (Ismatullah andEibert, 2009, 2008;Notaros, 2002;Eibert, 2005;Notaros, 2008) can be utilized.The definition of HO f n is in 2-D derived from 3-D HO a n , thus the system matrices from MoM and FE will be compatible based on the same geometrical structure information of the object.The system matrices from BI are described in Ismatullah and Eibert (2009).
operated on a Server with Intel(R) Xeon(R) CPU E5630 @ der of self-identified basis functions for FE-BI are compared with MIE scattering as shown in Fig. 3.The coated sphere consists of a PEC core with radius 0.9 m and a layer of dielectrics with thickness 0.1 m.The dielectric layer properties are given by ǫ r = 3 − 0.1j and µ r = 1.0.The incident wave is 550 MHz and propagating towards +z direction and the electric field is 100 V/m along x direction (E x = 100 V/m).
For 0th order, the mesh size was set to 0.04 m in Hypermesh software (HyperWorks ( 2012)), and the mean edge length is 4.65 cm, with minimum edge length 2.24 cm and 400 maximum edge length 8.11 cm.The total unknowns are 154 822.The running time was 1 314.8 s.For 1st order, the same mesh as for 0th order was utilized.The total unknowns tric layer are presented with ǫ r = 2.5 − 0.5j and µ r = 1.0.The incident wave is 3 GHz and propagating toward +z di-420 rection and the electric field is 100 V/m along x direction (E x = 100 V/m).The results for the RCS are shown in Fig. 4. 4. Bistatic RCS of coated sphere II @ 3 GHz on xz cut half plane (ϕ = 0 • ).

Linear algebraic equation system
To solve the electric field [u] and the magnetic field [i], the subsystem from FE in form of Eq. ( 1) and the subsystem generated by BI (Ismatullah and Eibert, 2009) must be combined as a complete system.The BI subsystem based on EFIE may introduce resonances into the final system, so it is necessary to utilize CFIE with similarly satisfiable accuracy.As a result, the subsystems can be regarded as where is the sub-matrix derived from FE-BI for corresponding u and i. M FE u is square, symmetrical and sparse, M FE i is rectangular and sparse, M BI u is rectangular and fully staffed, M BI i is square, symmetrical and fully occupied, V FE,BI are excitation vectors for FE-BI.Thus, the complete combined system is written as The complete system solves the electric and magnetic fields simultaneously, thus equivalent surface electric and magnetic currents can be achieved.

Numerical results
To testify the accuracy of the analytical matrix elements and the global matrix assembly in FE-BI, several numerical simulation results are shown in this section.A cogent demonstration is to utilize a coated sphere testing case, where a PEC sphere is enclosed by a layer of dielectric material.The analytical RCS is well known as MIE Scattering (Balanis, 1989).Good matching of RCS between analytical solution  and numerical method verifies the efficacy of FE-BI.As a more complicated testing case, a second sphere is displayed.With higher frequency and finer mesh, more unknowns are handled.Moreover, an example of FE-BI application in very large scale simulations is shown through the RCS of the Flamme aircraft.As 0th order of FE-BI has been verified in many published articles (Tzoulis and Eibert, 2005;Eibert and Hansen, 1997;Eibert, 2007), it can be utilized as a reference for HO FE-BI.Efficiency of FE-BI based on different orders of self-identified basis functions are presented.
The sphere simulations were performed on a PC with Intel(R) Core(TM)2 Quad CPU Q9550 @ 2.83 GHz processor, installed memory (RAM) 16.0 GB and 64-bit operating system.The simulation of the Flamme aircraft was operated on a Server with Intel(R) Xeon(R) CPU E5630 @ 2.53 GHz (2 processors), installed memory (RAM) 96.0 GB and 64bit operating system.All simulations were computed on one core.

Coated sphere I
For the coated sphere, the RCS from 0th, 1st and 2nd order of self-identified basis functions for FE-BI are compared with MIE scattering as shown in Fig. 3.The coated sphere of a PEC core with radius 0.9 m and a layer of dielectrics with thickness 0.1 m.The dielectric layer properties are given by r = 3 − 0.1j and µ r = 1.0.The incident wave is 550 MHz and propagating towards +z direction and the electric field is 100 V/m along x direction (E x = 100 V/m).For 0th order, the mesh size was set to 0.04 m in Hypermesh software (HyperWorks, 2012), and the mean edge length is 4.65 cm, with minimum edge length 2.24 cm and maximum edge length 8.11 cm.The total unknowns are 154 822.The running time was 1314.8 s.For 1st order, the same mesh as for 0th order was utilized.The total unknowns are 309 644.The running time was 4763.4 s.For 2nd order, the mesh size is set 0.1 m, and the mean edge length is 10.93 cm, with minimum edge length 6.48 cm and maximum The mesh size of LO is around λ/8 and the mesh size of HO is enlarged up to around λ/3. HO with coarser mesh and fewer unknowns achieves also accurate result as LO with finer mesh.

Coated sphere II
For the second coated sphere testing case, numerical RCS from different orders are also compared with MIE scattering.The second coated sphere contains a PEC sphere core with radius 0.5 m, and the PEC core is enclosed with a dielectric layer with thickness 0.0025 m.The properties of the dielectric layer are presented with r = 2.5 − 0.5j and µ r = 1.0.The incident wave is 3 GHz and propagating toward +z direction and the electric field is 100 V/m along x direction (E x = 100 V/m).The results for the RCS are shown in Fig. 4.
For 0th order, the mesh size was set to 0.01 m, and the mean edge length is 0.858 cm, with minimum edge length 0.250 cm and maximum edge length 1.631 cm.The total unknowns are 411 339.The running time was 3525.6 s.For 1st order, the same mesh as 0th order is utilized.The total unknowns are 822 678.The running time was 4271.3 s.For 2nd order, there are two testing cases.One uses the same mesh as 0th and 1st order.The total unknowns are 1 973 604.The running time was 7725.5 s.Another simulation utilizes mesh size 0.03 in Hypermesh, and the mean edge length is 2.554 cm.The minimum edge is 0.250 cm and maximum edge length 4.316 cm.The total unknowns are 202 522.The running time was 1357.6 s.
The numerical RCS results are compared with MIE.The mesh size of LO is around λ/8, while HO with finer mesh size is also working well though a lot more unknowns are  settled.When the mesh size of HO increases up to around λ/4, HO with coarser mesh for FE-BI maintains good precision as the results from LO with finer mesh as shown in Fig. 4.

Flamme
The Flamme case is an application of FE-BI for very large scale simulation.The Flamme is located in the xy plane, with nose heading along the +x axis, as shown in Fig. 5.The Flamme is enclosed by a layer of lossy dielectric material with thickness of approximately 1 cm.The permittivity of the dielectric material is r = 1.21 − 10j and the permeability is µ r = 1.The simulation frequency is 2.5 GHz.The incident plane wave propagates towards −x direction, with electric  field (E z = 100 V/m).To visualize absorbing effects of the lossy dielectric material, a PEC Flamme simulated with BI with 0th order of basis functions is utilized for comparison.The RCS of PEC and layered Flamme in different cut planes are shown in Figs.6-10.The PEC Flamme is simulated through BI with 0th order of self-identified basis functions, the layered Flamme is simulated through FE-BI with 0th, 1st and 2nd order of self-identified basis functions.As the efficacy of 0th order with finer mesh has been verified, here it is used as a reference.The RCS comparison shows that most of input power goes over the Flamme.In scattered directions, the scattered power is evidently absorbed by the dielectric material.For PEC, the mesh size of the PEC Flamme was set to 0.01 m, and the mean edge length is 1.011 cm, with minimum edge length 0.100 cm and maximum edge length 2.559 cm.The total unknowns are 692 952.The number of BI electric currents is 692 952 and the number of BI magnetic currents is 0. The number of levels for MLFMM is 8 and the peak memory consumption was 4083.746MBytes.The running time was 75 529.3s.For 0th order, the mesh size of the layered Flamme was set to 0.01 m, and the mean edge length is 1.083 cm, with minimum edge length 0.044 cm and maximum edge length 2.671 cm.The total unknowns are 2 081 547.The number of BI electric currents is 690 960 and the number of BI magnetic currents is 602 284.The number of levels for MLFMM is 8 and the peak memory consumption was 11 384.43MBytes.The run time was 23 511.6 s.For 1st order, the mesh size of layered Flamme was set to 0.02 m, the mean edge length is 1.770 cm, with minimum edge length 0.065 cm and maximum edge length 3.933 cm.The total unknowns are 1 252 430.The number of BI electric currents is 427 144 and the number of BI magnetic currents is 360 976.The number of levels for MLFMM is 7 and the peak memory consumption was 8825.734MBytes.The running time was 12 141.3 s.For 2nd order, the mesh size was set to 0.02 m, and the mean edge length is 1.770 cm, with minimum edge length 0.065 cm and maximum edge length 3.933 cm.The total unknowns are 2 941 242.The number of BI electric currents is 712 284 and the number of BI magnetic currents is 602 004.The number of levels for MLFMM is 8 and the peak memory consumption was 18 536.68MByte.The running time was 47 871.3 s.

Conclusion
Self-identified hierarchical 3-D vector basis functions were generated for the hybrid finite element (FE) -boundary integal (BI) technique, where analytical solutions for the FE matrix elements have been presented up to 2nd order.Self-identified basis functions provide feasibility for FE and effectively maintain compatibility with BI.Going from 1st to 2nd order, FE-BI allows for a mesh size increase from λ/8 up to λ/3.From coated sphere testing cases, good accuracy was found and the Flamme simulations displayed that FE-BI based on self-identified basis functions can be applied for very large scale simulations.

Figure 1 .
Figure 1.The definition of subscripts based on a single tetrahedron.Every vertex index can be used as a row index m or as a column index n.

Figure 2 .
Figure 2. The general geometrical configuration for finite elementboundary integral method.

Figure 5 .
Figure 5.The geometry of a Flamme airplane.

Figure 11 .
Figure 11.Electric surface current density |J | in A/m distribution of a covered Flamme airplane with plane wave incidence @ 2.5 GHz.
2 .Properties of identified hierarchical basis functions are displayed in Table 1.(i, j ) and (r, s, t) in Table 1 contain the local node numbers in the series defined by global order.Also, the curl of identified hierarchical basis functions is given.

Table 1 .
Hierarchical basis functions and properties within the tetrahedron

Table 2 .
Analytical solution for [R mn ] matrix elements.

Table 3 .
Analytical solution for [S mn ] matrix elements.