Assessing Curl-Conforming Bases for Pyramid Cells

Successful three-dimensional finite element codes for Maxwell's equations must include and deal with all four types of geometrical shapes: tetrahedra, bricks, prisms, and quadrangular-based pyramids. However, pyramidal elements have so far been used very rarely because the basis functions associated with them have complicated expression, are complex in derivation, and have never been comprehensively validated. We recently published a simpler procedure for constructing higher-order vector bases for pyramid elements, so here we fill a gap by discussing a whole set of test case results that not only validate our new curl-conforming bases for pyramids, but which enable validation of other codes that use pyramidal elements for finite element method applications. The solutions of the various test cases are obtained using either higher order elements or multipyramidal meshes or both. Furthermore, the results are always compared with the solutions obtained with classical tetrahedral meshes using higher order bases. This allows us to verify that purely pyramidal meshes and elements give numerical results of comparable accuracy to those obtained with multitetrahedral meshes that use elements of the same order, essentially requiring the same number of degrees of freedom. The various results provided here also show that higher order vector bases always guarantee a superior convergence of the numerical results as the number of degrees of freedom increases.

Assessing Curl-Conforming Bases for Pyramid Cells Roberto D. Graglia , Life Fellow, IEEE, and Paolo Petrini , Member, IEEE Abstract-Successful three-dimensional finite element codes for Maxwell's equations must include and deal with all four types of geometrical shapes: tetrahedra, bricks, prisms, and quadrangularbased pyramids.However, pyramidal elements have so far been used very rarely because the basis functions associated with them have complicated expression, are complex in derivation, and have never been comprehensively validated.We recently published a simpler procedure for constructing higher-order vector bases for pyramid elements, so here we fill a gap by discussing a whole set of test case results that not only validate our new curl-conforming bases for pyramids, but which enable validation of other codes that use pyramidal elements for finite element method applications.The solutions of the various test cases are obtained using either higher order elements or multipyramidal meshes or both.Furthermore, the results are always compared with the solutions obtained with classical tetrahedral meshes using higher order bases.This allows us to verify that purely pyramidal meshes and elements give numerical results of comparable accuracy to those obtained with multitetrahedral meshes that use elements of the same order, essentially requiring the same number of degrees of freedom.The various results provided here also show that higher order vector bases always guarantee a superior convergence of the numerical results as the number of degrees of freedom increases.
Index Terms-Electromagnetic fields, finite-element methods, higher order vector elements, pyramidal elements, numerical analysis.

I. INTRODUCTION
A LMOST all practitioners of computational electromagnet- ics (CEM) avoid using pyramidal cells because the related vector bases are very complicated and never extensively tested; in fact very few authors have used pyramids so far, as can be seen from [1], [2], [3] and references therein.This has hindered the development of codes that use hybrid meshes smoothly; that is, meshes that employ higher-order pyramidal elements in addition to tetrahedra, bricks, and prisms, despite the fact that a reasonably extensive scientific literature on higher-order pyramidal elements has developed over the past twenty years [4]- [11].Things may now change, firstly because the method in [1], [2] to obtain conforming bases of arbitrarily high order for pyramids is, The authors are with the Dipartimento di Elettronica e Telecomunicazioni, Politecnico di Torino, 10129 Torino, Italy (e-mail: roberto.graglia@polito.it;paolo.petrini@polito.it).
Digital Object Identifier 10.1109/JMMCT.2023.3333563 in our opinion, decidedly simpler than those previously available in the literature.Secondly because this paper fills the lack of useful test cases for validating numerical codes that use pyramid cells, or at least their curl-conforming bases for Finite Element Method (FEM) applications.In fact, this paper assesses the modeling capabilities of the higher order bases for pyramid cells reported in [1], which differ from those in [3], presenting several results for pyramidal and rectangular cavity resonators useful for validating computer codes that use these cells and bases.This is done because cavity problems are the best benchmark for immediately checking whether a (new) basis avoids spurious modes or not.More specifically, this paper considers rectangular electromagnetic resonators for two very good reasons: first because they have a known and very simple solution, and second because their geometry can be meshed using only pyramids.The results obtained here add to those in [1] where we modeled a few pyramidal cavities with a single cell.
The fields inside the volume V of a closed metal cavity and the corresponding wavenumbers k are modeled by the vector Helmholtz equation whose discretization produces a generalized matrix eigenvalue equation Ae = k 2 Be with entries [12, Sections 6.5, 6.6] and B m and B n being vector basis functions.In the following we assume the material inside the cavity to be homogeneous with unitary relative permeability μ r and relative permittivity ε r .
The wavenumbers k are computed by means of a FORTRAN code that uses the LAPACK library and the routine DSPGV.In common with standard finite element implementations, the system of equations is usually constructed in a cell-by-cell manner where the integrals in (1) and ( 2) are evaluated throughout a single cell for all combinations of basis and testing functions, stored in a temporary "element matrix," and systematically transferred to the global system of equations.
This FEM technique can be applied to any conforming mesh made with a mixture of cells of different size and shape (that is tetrahedrons, prisms, bricks and pyramids), provided that curlconforming bases are used, as required by (1).Obviously, we refer to implementations that use also higher order bases.
Recall that all elements have shape functions and basis functions of polynomial form in the so-called "parent" space [12], with the exception of pyramid elements which must have shapefunctions and basis functions of fractional form for conformity with the other cells and to guarantee cell-to-cell tangential continuity, even in the case of curved cells [1].(In the parent space, the divergence-conforming basis functions of the pyramid element have fractional form to guarantee the continuity of the normal component of the field from cell to cell [2]).
For pyramids, the shape functions and basis functions take polynomial form in a different space, called "grandparent" space, where all the pyramid cells are mapped by a cube of unit side [1], [2].The polynomial order of the functions in the grandparent domain for the pyramid elements and in the parent domain for the other non-pyramidal elements is denoted by the integer p.For any polynomial order p, the curl-conforming functions are tangentially continuous on the boundaries in common to adjacent cells of different shape [1].
For greater clarity we recall that the use of pyramidal elements requires two mappings.The first from parent space to the observer's space (the so-called child space).The second from grandparent space to parent space.The Jacobian of the first mapping depends on the shape of the pyramidal cell in the child space, that of the second does not and is always given by the same polynomial term; see [1, eq. (11)].As happens for other elements, the unitary basis vectors of the pyramidal cell are obtained by differentiating the element position vector r with respect to three independent parent coordinates [1].From these vectors, as shown in [1 Table I], one gets the gradient vectors and the edge vectors that appear in the expression of the vector basis functions and their curls.Now, unlike what is usually done with the other elements, the FEM-matrix entries (1, 2) associated to pyramids are always computed by integrating polynomials in the grandparent domain [1], and not by integrating the "corresponding" fractional functions in the parent domain.However, the software used to produce the results presented in this paper is not optimized; that is, we trivially compute the volume FEM integrals by numerical integration on the unit cube of the grandparent domain along three directions, in cascade, without using sophisticated integration schemes that can reduce computation times [9], [13].
The rest of the paper is structured as follows.Section II provides several results related to a pyramid-shaped cavity resonator with unit sides, meshed with a single cell or multiple pyramidal cells of different aspect ratios.The numerical results are validated against those of a different code that uses higherorder tetrahedral elements.Section III reports several results for rectangular cavity resonators meshed by pyramidal cells and compares them with the well-known exact analytical solution.Also, in Section III, the accuracy of the results obtained using only pyramidal cells is compared with that obtained from the code using higher order tetrahedral elements.In summary, in this paper we use meshes formed only by quadrangular-based pyramidal cells; unlike what happens when pyramids are used (albeit rarely) as "fillers" of hybrid meshes; that is, as gluing elements inserted between other non-pyramidal cells.Readers may find it helpful to review [1], [2], for a detailed introduction to the notation and other background information.Preliminar results of this work were presented in [14].

II. RESULTS FOR THE EQUILATERAL PYRAMID
Table I reports the first six wavenumbers of an equilateral pyramid-shaped cavity obtained with bases of order p from 0 to 6, while Fig. 1 shows the trend of the condition number of the mass-matrix as the order used varies.The results obtained by meshing the cavity with a single pyramidal cell are compared in Table I with those obtained using four identical tetrahedral Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.I results obtained for p between 0 and 6, and are computed against the WN result reported in bold in the middle column of Table I.The second and third modes of the cavity have identical wavenumbers, as do the fourth and fifth modes.For p = 0, a single pyramid cell yields only four non-zero eigenvalues; for this p = 0 case we show with a blue square marker the errors of the first three modes only.cells, as done in [1 Table VII] for bases up to third order.The Table results show that r for p ≥ 3 we can sort the first six wavenumbers of the cavity in ascending order and with the correct multiplicity, regardless of the mesh we have used; r for p ≥ 5 the number of zero eigenvalues obtained with a single pyramid cell is higher than that obtained with 4 tetrahedral cells; r for p ≥ 6 the number of degrees of freedom (DoF) obtained by meshing the structure with a single pyramidal cell is higher than that required by meshing with four identical tetrahedrons.In essence, Table I makes it clear that to obtain good results for simply shaped cavities we must use at least third-order bases when using very few cells.Furthermore, in our case, the mesh composed of 4 identical tetrahedra does not destroy the symmetry of the modes and is therefore able to use the DoFs in an optimal way, provided that bases of order higher than the fifth are used.In principle we can expect to obtain results as good as those shown in Table I (for p = 6) by using many more cells than done here, and bases of order lower than the third.Unfortunately, this turns out to be completely inconvenient because tens or hundreds of thousands of DoFs would be needed precisely because the number of DoFs that guarantees good results does not depend linearly on the order of the base used.On the other hand, we cannot think to improve the results simply by increasing the order of the bases used beyond the sixth because, as shown in Fig. 1, the condition number of the mass matrix (CN) grows exponentially with the order of the bases in use (regardless of the shape of the cells), and is already of the order of 10 6 for sixth-order bases.To improve the accuracy of the results of Table I (for p = 6) it becomes necessary to refine the mesh rather than further increase the order of the base used.
In any case, we are not interested in determining whether the best results in Table I for sixth-order bases are those obtained with the single pyramidal cell model rather than those obtained with four tetrahedral cells.What really matters is that we have never detected spurious modes using a single pyramidal cell.Since an analytical result is missing, we assume that for p = 6 the two meshes lead to results of equal accuracy; that is, the reference wavenumber (WN) in the central column of Table I is simply the average of the two values obtained using the two different meshes with p = 6.As shown in Fig. 2, the results obtained using a single cell and those obtained with four identical tetrahedral cells converge similarly to the reference WN as p and the number of DoFs increase.
The absence of spurious modes using a single cell needs to be verified considering structures meshed with multiple cells.Therefore, let us now consider the equilateral pyramid of Fig. 3, meshed with four pyramids of different aspect ratio which varies by varying the parameter s v defined in the figure caption.The base of the yellow pyramid of Fig. 3 is always a square for 0 < s v ≤ 1/2, therefore its Jacobian is constant at all points inside it.Conversely, the Jacobian inside the other three pyramids is polynomial (not constant), unless s v = 1/2.Note that the yellow pyramid disappears for s v = 0.This limiting case is not considered in the sequel because it relates to a mesh formed by one pyramid (the one in blue in Fig. 3) and two tetrahedrons (in white).
Fig. 4 reports, versus s v , different results for the pyramidal cavity obtained with the meshes described in Fig. 3.The figure at top shows the percentage error on k for the first three modes, Fig. 3. On the left, an equilateral pyramid is meshed in various ways by four rectilinear pyramids of different aspect ratios that we define by assigning the position of their five vertices in the observer's space.The central column shows the pyramids seen from above, the right one the bases of the pyramids.In addition to the common apex, the four pyramids of the mesh also share a vertex located at the coordinate point s = s v along the diagonal going from southwest to northeast in the figures of the right column.The pyramid in yellow has a square base whatever the value of s v , for s v between 0 and 1/2.The meshes in the upper and lower row are obtained with s v = 1/2 and s v = 1/8, respectively.(The coordinate s v is normalized w.r.t. the length of the diagonal of the base of the equilateral pyramid).
obtained with the base of order p = 4.In the center the percentage error on k for the first mode only, obtained with bases of order 2, 3 and 4. Fig. 4 at bottom the condition number of the mass matrix as s v varies, obtained with bases of order 2, 3 and 4.
The results of Fig. 4 show that the wavenumbers of the second and third modes coincide for s v = 0.5 (symmetric mesh), and therefore have the same percentage error.This no longer happens if the mesh loses symmetry (s v < 0.5) while, for sloppier meshes (s v < 0.1) the errors increase, albeit slightly, and the CN increases (worsens) by a factor of the order of 10 3 .These results prove the robustness of our pyramid bases.

III. RESULTS FOR RECTANGULAR CAVITY RESONATORS
Let us consider a rectangular cross-sectional waveguide of length d closed by metal walls at both ends, that is the cavity resonator {x, y, z: 0 ≤ x ≤ a; 0 ≤ y ≤ b; 0 ≤ z ≤ d}.For this kind of cavities, the cartesian components of the eigenfields e and h are the product of three sinusoidal functions, one for each Cartesian variables x, y, z [15], [16].For example, the longitudinal components e z , h z are for the TE and the TM mode, respectively.Regardless of whether it is a TE or TM mode, the mode wavenumber depends on the integer mode-numbers m, n, p, although there are no TE modes with mode index p equal to zero, nor TM modes with mode index m or n equal to zero, while there are  .The results are obtained using the base of order p = 4 with a structured mesh made of only six pyramids, one for each face of the cavity, having a common apex at . This model produces 2,020 DoFs and 640 zero eigenvalues, with a mass matrix condition number of the order of 1.1 × 10 5 after diagonal preconditioning.The percentage error for the first 20 modes is shown in Fig. 5.The error is less than 0.1% for the first 18 modes and less than 0.5% for the first 37 modes.No spurious modes are observed.(We show the results for order p = 4 because this is a good compromise to obtain reduced computation times and achievable precision.) Since each face of the cavity is the quadrilateral surface on which only one pyramid rests, the quality of the results depends on how well the basis functions approximate, inside each pyramid, the sinusoidal trend of the eigenfields in the three Cartesian Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.II obtained with the pyramidal base of order four and a mesh of six pyramids.The figure reports the percentage error on the computed wavenumber with respect to the exact wavenumber (5), identified by the mode numbers m, n, p.These indices are reported only for some modes simply to clarify the trend of the error as their value increases.
and [0 ≤ z ≤ pπ].In fact, Fig. 5 clearly shows that every time the index m, n, or p increases by one unit, the accuracy of the results degrades.
To obtain results of quality close to that of Fig. 5 with zeroorder basis functions and a structured mesh made up of pyramids we would need at least 10 7 unknowns which we may obtain by subdividing the cavity into N × N × N bricks, with N ≈ 100 (see [12, §1.3]), and then define the usual 6 pyramids in each brick.(It goes without saying that if you are looking for quality results for the modes of an arbitrarily shaped cavity it is generally more convenient, in terms of CPU time, to immediately use brick or tetrahedral elements rather than pyramid cells everywhere.)The results that follow refer only to the cube-shaped cavity and are obtained with meshes of the type shown in Fig. 6, obtained by moving the apex in common to the 6 pyramids along the diagonal of the cavity of length √ 3.This model continues to have 2,020 DoF, but now the precision of the results is a function of the distance t between the apex and the corner, measured along the aforementioned diagonal by the normalized coordinate s = t/ √ 3. Fig. 7 at top shows the error on k for the three fundamental modes of the cubic cavity (i.e., modes 101, 011, and 110), and for the next two modes (TE and TM 111).At bottom we show the mass-matrix condition number versus s.
So far we have shown results for the cube-shaped cavity obtained with a few cells in the mesh.As mentioned, we can increase the number of cells by replicating the "basic cube" several times, after having divided it into 6 pyramids, or into 5 tetrahedra as in Fig. 8.If we use sub-cubes divided into 5 tetrahedra, we must rotate each of them appropriately to guarantee mesh conformity.The average percentage error of the first 5 modes of the cubic cavity vs the number of DoFs is shown in Fig. 9 for varying base order.The straight-lines in Fig. 9 indicate the speed of convergence as the number of DoFs Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and the polynomial order of the elements used increase.Fig. 9 shows that the numerical solutions with tetrahedral or pyramidal mesh are more or less equivalent.However, we point out that, given the same order and (almost the same) number of DoFs, the CPU times required by a pyramidal mesh are higher than those required by a mesh entirely made up of tetrahedra.
Finally, for completeness, Fig. 10 shows the mass-matrix condition numbers relative to the cases in Fig. 9 versus the number of degrees of freedom.Note that in this case the CN obtained with the pyramidal meshes essentially depend only on the order of the base used because they relate to structured meshes formed by perfectly identical pyramids arranged in a perfectly symmetrical way, as shown in Fig. 6 in the center.Conversely, the CNs obtained with tetrahedral meshes increase as the number of degrees of freedom increases since these meshes are formed by tetrahedra of different aspect ratios, as shown in Fig. 8.

IV. CONCLUSION
This paper presents a whole set of test case results validating the new higher-order curl-conforming bases for pyramidal cells we have recently published.These basis functions are built to guarantee the continuity of the tangential components between adjacent elements of the same order but different shape, so that it is possible to build more efficient computer codes that use hybrid meshes made with tetrahedra, bricks, prisms and quadrangular-based pyramids.The solutions of the various test cases discussed here are obtained using higher order elements or multipyramidal meshes or both, and the results are always compared with the solutions obtained with classical tetrahedral meshes using higher order bases.This allows us to verify that purely pyramidal meshes and elements give numerical results of comparable accuracy to those obtained with multitetrahedral meshes that use elements of the same order, essentially requiring the same number of degrees of freedom.The reported numerical examples show once again that higher order functions provide more accurate results than those obtainable with lower order elements.It is hoped that these results will facilitate the development of new electromagnetic solvers based on the use of hybrid meshes.

Manuscript received 20
September 2023; revised 9 November 2023; accepted 13 November 2023.Date of publication 16 November 2023; date of current version 24 November 2023.This work was supported in part by PRIN under Grant 2017NT5W7Z GREEN TAGS, and in part by the European Union under the Italian National Recovery and Resilience Plan of NextGenerationEU through Program "Telecommunications of the Future" (PE00000001-program "RESTART") and the Program PNRR M4C2 "Multiscale modeling and Engineering Applications" of the National Centre for HPC, Big Data and Quantum Computing under Grant HPC-CUP E13C22000990001.(Corresponding author: Roberto D. Graglia.)

Fig. 1 .
Fig. 1.Condition Number (CN) of the mass-matrix versus the number of degrees of freedom obtained using bases of order p from 0 to 6 by meshing an equilateral pyramid with a single pyramidal element (in blue), or with four identical tetrahedral elements (in red).

Fig. 2 .
Fig. 2. Percentage errors for the first 6 wavenumbers of an equilateral pyramid versus the number of degrees of freedom.Errors refer to TableIresults obtained for p between 0 and 6, and are computed against the WN result reported in bold in the middle column of TableI.The second and third modes of the cavity have identical wavenumbers, as do the fourth and fifth modes.For p = 0, a single pyramid cell yields only four non-zero eigenvalues; for this p = 0 case we show with a blue square marker the errors of the first three modes only.

Fig. 4 .
Fig. 4. Results obtained with the meshes in Fig. 3 by varying the s v parameter.

Fig. 5 .
Fig. 5. Results for the rectangular cavities of TableIIobtained with the pyramidal base of order four and a mesh of six pyramids.The figure reports the percentage error on the computed wavenumber with respect to the exact wavenumber (5), identified by the mode numbers m, n, p.These indices are reported only for some modes simply to clarify the trend of the error as their value increases.

Fig. 6 .
Fig. 6.A cube is meshed in different ways by six square-based rectilinear pyramids by varying the position of their common apex.In the figure the pyramids resting on the east and west faces of the cube are in yellow, the others are transparent.The common apex is located at the normalized coordinate point s on the diagonal joining the southwest vertex of the cube, where s = 0, to the northeast vertex, where s = 1.The six pyramids are identical if s = 1/2, as shown in the center.The left and right figures show the cases of s = 1/4 and s = 3/4, respectively.

Fig. 7 .
Fig. 7. Figure considers the cube cavity meshed with six pyramids having a common apex that varies along the diagonal of the cube of length √ 3, as described in Fig. 6.The results obtained with the fourth order base are shown with respect to the normalized apex-corner distance s.At top we show the percentage error on k for the first five modes; at bottom the CN mass matrix after diagonal preconditioning.

Fig. 8 .
Fig.8.A unit cube is divided into five tetrahedra by cutting off every other vertex; that is, four out of eight vertices, as shown in the figure on the left.On the left, the blue colored tetrahedron in the center of the cube is regular; that is, its edges have equal length.In the center and on the right we show the remaining tetrahedra two at a time, in different colors.

Fig. 9 .
Fig. 9. Average percentage error for the first five modes of a cubic cavity obtained with elements of different order p. Meshes of different densities are obtained by first dividing the cavity into cubes, and then each cube into five tetrahedra (top figure) or six pyramids (bottom figure).The straight lines indicate the speed of convergence as the number of degrees of freedom and the polynomial order of the elements increase.The results at top were first presented in [17].

Fig. 10 .
Fig. 10.Mass-matrix condition number related to the cases in Fig. 9. On the left the results obtained with tetrahedral cells, on the right those with pyramidal cells.

TABLE I FIRST
SIX WAVENUMBERS OF THE EQUILATERAL PYRAMID-SHAPED CAVITY