Hierarchical Curl-Conforming Vector Bases for Pyramid Cells

Advanced applications of the finite element method use hybrid meshes of differently shaped elements that need transition cells between quadrilateral and triangular faced elements. The greatest ease of construction is obtained when, in addition to triangular prisms, one uses also pyramids with a quadrilateral base, as these are the transition elements with the fewest possible faces and edges. A distinctive geometric feature of the pyramid is that its vertex is the point in common with four of its faces, while the other canonical elements have vertices in common with three edges and three faces, and that is why pyramids' vector bases have hitherto been obtained with complex procedures. Here we present a much simpler and more straightforward procedure by shifting to a new paradigm that requires mapping the pyramidal cell into a cube and then directly enforcing the conformity of the vector bases with those used on adjacent differently shaped cells (tetrahedra, hexahedra and triangular prisms). The hierarchical curl-conforming vector bases derived here have simple and easy to implement mathematical expressions, including those of their curls. Bases completeness is demonstrated for the first time, and results confirming avoidance of spurious modes and faster convergence are also reported.


I. INTRODUCTION
Very accurate three-dimensional (3D) models that balance computational efficiency with geometric flexibility are obtained using hybrid meshes of tetrahedra, hexahedra (bricks) and prismatic cells together with interpolatory [1], [2] or hierarchical vector bases [3], [4], [5] of high order. All these divergence-and curl-conforming bases are now reported in a single book [6]. Curl-conforming bases with continuous tangential components across adjacent cells in the mesh are used in Finite Element Method (FEM) applications to discretize the vector Helmholtz operator. Hierarchical bases are preferred in applications that implement p-adaptive refinement techniques because they easily allow for selective expansion using a different order in different regions of the computational domain, which often leads to reduced computation time and more accurate results [7], [8], [9]. It is difficult to create conforming hybrid-meshes without pyramid elements. (The definition of conforming mesh is given, for example, in [6,Chap. 2].) Pyramids with a quadrilateral and four triangular faces are in fact the most obvious natural fillers for discretizations mostly formed by other differently-shaped cells. For example, pyramids are required when one has to link a coarse to a dense mesh of bricks (see [10,Fig. 2] and [11, Fig. 1]). This paper presents the pyramid's bases that conform to and complete the hierarchical curl-conforming families reported in [6]. Our new basis functions are obtained directly, simply by multiplying the lowest-order vector functions given in [11] by a suitable set of generating scalar polynomials thereby using, in essence, the same technique already used in [3], [4], [5], [6] to obtain the bases for the other 3D elements. More importantly, this paper provides ready-to-use expressions of the vector basis functions, normalized to obtain relatively low mass-matrix condition numbers (CN) in FEM applications.
As with the other hierarchical volumetric elements in [6], our new pyramid bases have four distinctive features: (a) the vector basis functions are subdivided from the outset into three different groups of edge, face, and volume-based functions; (b) each basis function is obtained by using one generating edge, face or volume-based scalar polynomial whose analytical expression involves all the dependent variables that describe the cell; (c) the generating polynomials are either symmetric or antisymmetric with respect to the variables that describe each edge and face of the cell, and are organized hierarchically; (d) in each group, the generating polynomials are mutually orthogonal for inner product defined by the integral on the volume, the face or the edge of the cell.
To put our work in perspective, note that in the late 1990s practitioners were satisfied with using pyramidal bases at most of the first order [11], [12], [13], as it was not clear the number of volume-based functions needed for higher orders. Research in this area has then received a new impetus from 2010 onwards, with important results published mainly in applied mathematics journals (see [14]- [18] and references therein). Unfortunately, the specialized literature has privileged theoretical aspects so that it is difficult to extract from this literature ready-to-use recipes for computational applications. In fact, until now, the use of higher order pyramid bases has been hampered by the complex procedures used to generate them. Also, as discussed in the present paper, in the volume, inside the pyramid, the vector basis functions have a fractional form that is more complicated than that proposed in [11], as it became apparent in other papers [13], [15], [16], [17], [18] published after [11]. Indeed, as noted in [19], for orders higher than the first, the bases in [11] lead to spurious modes because they are constructed with multiplicative polynomials of the parent variables that vary linearly in the so called parent pyramidal cell. The correct multiplicative functions are obtained here by working in a different space of scaled parent variables, which we call the grand-parent space, where the pyramidal cell is mapped by a cube on which we compute the FEM volume integrals; in the grandparent cube the vector basis functions and their curl, as well as the multiplicative functions and their gradient, all have polynomial form.
Although the bases in [11] can yield spurious modes they are complete in their space and, on the pyramid boundary, they can take the same "polynomial" expression as those of differently shaped elements. This has pronted us to reverse the derivation technique adopted elsewhere [17] as the edge and face based functions reported here have been derived starting from the known expression of the functions of adjacent, differently shaped elements (so to speak, by imposing a tangential continuity at the boundary of the cell), while the volumebased functions are obtained thanks to simple analytical and geometric considerations. Our volume-based functions differ from those given in [16], [17]; their order does not exceed that of the edge-and face-based vector functions whose number and order is dictated by continuity with adjacent elements. In our opinion, this amounts to saying that the main problem that researchers have had up to now to build the pyramid's bases is to find the simplest and most direct way to build the volumebased vector functions, and then prove the completeness of the whole family which, to the best of our knowledge, is given here for the first time in the Appendix. The proof in the Appendix clarifies the correct order the multiplicative polynomials must reach in the grand-parent space to generate the vector basis functions of the desired order.
The remainder of the paper is structured as follows. In Section II we explain why the basis functions and their curls shall belong to a common (polynomial) subspace to be identified from the outset. Section III describes the geometric representation of the pyramidal element while the variables used to define the polynomial bases are discussed at length in Section IV, which also reports the lowest order base. The higher order bases are given in Section V, while numerical results are reported in Section VI. Readers may find it useful to review [11] for a detailed introduction to the notation and other background information, as well as the fundamental paper [17] that inspired ours.

II. BASIS FUNCTIONS' SUBSPACES
Polynomial vector bases accompany tetrahedral, brick and prismatic cells in a completely natural way because only three edges and three faces branch off from the vertices of a tetrahedron, a brick and a triangular prism. Conversely, since the pyramid has four edges and faces converging at one vertex, it is necessary to explain why to look for polynomial bases for the pyramid as well. In this regard, we observe that FEM problems are usually formulated in terms of a single main unknown: either the magnetic field H or the electric field E. The other field is determined by the curl of the main field, as imposed by one of the Maxwell's curl equation. For example, to clarify, the electromagnetic fields that can exist within a 3D cavity formed by perfectly conducting walls, containing homogeneous or inhomogeneous dielectric Out[1642]= Fig. 1. The shape functions map the parent pyramid on the left to the child pyramid in the center. In the grandparent space (η, ξ 5 ) the shape functions and the basis functions take polynomial form, while the pyramid is described by the cubic cell shown on the right. and magnetic materials, may be found by solving one of the following vector Helmholtz equations [6] where ϵ r (x, y, z) and µ r (x, y, z) denote the relative permittivity and relative permeability functions, respectively. (Section VI reports some results for pyramid-shaped cavities obtained by solving Helmholtz's equation numerically.) FEM formulations in terms of H or E do not show significant differences except for the different way used to impose the boundary conditions. In fact, curl-conforming bases are derived with no reference to the main field we want to represent since E and H swap in the dual problem. In all 3D FEM applications the curl-conforming basis functions Ω are expressed in the most convenient way, i.e. using three coordinate gradient vectors ∇ξ a , ∇ξ b , ∇ξ c (which are curl-free); the curl of the basis functions can then be expressed in terms of unitary basis vectors ℓ a , ℓ b , ℓ c (see [1], [6], and the following Section III). Thus, for the pyramid, we may write where A i and B i turn out to be non-polynomial functions of the parent variables. Without specifying the subspace to which they belong, here we can only require that (3) and (4) are finite energy functions. Of course this is too little since H and E are interchangeable and constrained by Maxwell's curl equations, and especially because we may wish to rewrite Ω in terms of the unitary basis vectors. It is therefore reasonable to expect and require that vector basis functions and their curls belong to a much better specified common subspace, perhaps a polynomial one as it happens for all the other elements.
Since we know the lowest order pyramidal basis functions from [11], in Section IV we easily identify the polynomial subspace in common to the basis functions and their curls. There we see that the coefficients A i and B i in (3,4) are polynomials of variables other than the parent ones, although related to them. The higher order volume-based vector functions are obtained by inspection once we understand which variables we must use to get vector bases of polynomial form.

III. ELEMENT GEOMETRY REPRESENTATION IN CHILD
AND PARENT SPACE With reference to Fig. 1, any child pyramid in the observer's space r = (x, y, z) is a mapping of the same parent pyramid obtained through shape functions of five parent variables ξ = (ξ 1 , ξ 3 ; ξ 2 , ξ 4 ; ξ 5 ) that vary linearly across the element. The faces are numbered to correspond to the indexing of the associated parametric coordinates; that is, the j-th face of the pyramid lies on the parametric surface ξ j = 0. The quadrilateral face lies on the coordinate surface ξ 5 = 0, while the four triangular faces have coordinate ξ j = 0, with j=1 to 4. We choose as independent coordinates ξ 1 , ξ 2 and ξ 5 , so that ∇ξ 5 · (∇ξ 1 × ∇ξ 2 ) is strictly positive (see Fig. 1), while ξ 3 and ξ 4 are dependent coordinates. Gradient (∇ξ a ) and edge (ℓ ab ) vectors are defined in the child-space and calculated as summarized in Table I. Notice that the triangular faces 1 and 3 are opposite to each other and have in common only the vertex of the pyramid, as happens for faces 2 and 4. Henceforth, the subscripts that label the triangular faces are always computed modulo 4 meaning, for example, ξ j−1 = ξ 4 for j = 1, and ξ γ+1 = ξ 1 for γ = 4. Similarly, the edges are given a twoindex label deriving from the two faces sharing the edge. The pyramid's edges shared by two triangular faces are called triangle edges to distinguish them from the mixed edges lying on the quadrilateral face and in common to only one triangular face; this naming is the same used in [17]. The edges indicated by two dummy indices γ and γ + 1 are always triangle edges and, to lighten the notation, we will often set The ξ variables are linked by the dependence relations [11] ξ 1 + ξ 3 + ξ 5 = 1 ξ 2 + ξ 4 + ξ 5 = 1 that hold also on each face of the cell. In fact, the parametric coordinates ξ 1 , ξ 3 and ξ 5 used to describe the triangular face ξ 2 = 0 (as well as the ξ 4 = 0 face) are dependent on each other as in the first of (6). The coordinates ξ 2 , ξ 4 and ξ 5 which describe the faces ξ 1 = 0 and ξ 3 = 0 are dependent on each other as in the second of (6). Finally, for the quadrilateral face ξ 5 = 0, one has the two dependency relations obtained from (6) by setting ξ 5 = 0. On the pyramid's tip one has ξ 5 = 1 while eqs. (6) force the remaining four coordinates to vanish.

IV. LOWEST ORDER BASIS FUNCTIONS AND NEW PARADIGM TO DERIVE HIGH ORDER BASES
One of the main objectives of this Section is to clarify which variables to use to get pyramidal bases of polynomial type, because certainly these bases are not of polynomial type if we specify them using the parent variables introduced in the previous Section. Thus, in addition to the parent coordinates, we introduce four scaled coordinates, for j = 1, 2, 3, 4 with dependence relations (for ξ 5 ̸ = 1) stemming from (6) These coordinates, used also in [17], transform surface integrals on the triangular face ξ j = 0 of the parent pyramid (i.e. a triangular simplex) into integrals on a unit square, while volume integrals on the parent pyramid become integrals on the unit cube of the (η, ξ 5 )-space, as shown below in (10), The kernel on the right side of (11) could be used to cancel singular terms that the function G(η, ξ 5 ) could have at the vertex of the pyramid, in ξ 5 = 1. For example, tip singularities could occur if the tip of the pyramid is an isolated vertex of the structure under consideration. To model physically permitted singularities one should develop and use singular basis functions, which are beyond the scope of this paper.
Most important of all is that the vector basis functions and their curl take a simple polynomial form in terms of the grandparent variables (η, ξ 5 ), as summarized in Table  I, whereas they have fractional form when using ξ parent variables [11]. The same happens to the shape functions; for example, the interpolatory shape functions used in [11] take a polynomial form by using (7) to replace the parent variables with the new grandparent ones.
In this regard, we immediately claim the fundamental result used to derive our bases which stems from (8), namely that the curl of any linear combination of terms such as η α 1 η β 2 ξ z 5 (1 − ξ 5 )∇ξ a (where the subscript a in ∇ξ a is 1, 2, or 5), and the gradient of any linear combination of terms such as η α 1 η β 2 ξ z 5 (1 − ξ 5 ) take a polynomial form in the space (η, ξ 5 ). (Of course, each term of these linear combinations can have different values of the exponents α, β, and z.) With reference to Table I, the factor (1 − ξ 5 ) 2 guarantees polynomial form to Ω γ5 and its curl, while the proof that the polynomials Ω γδ agree with the above fundamental result is obtained by highlighting their non-zero curl component which gives (for γ = 1 to 4, and δ = γ + 1) The curl-free polynomial vector on the right of (12) guarantees tangential continuity between adjacent elements and cannot be omitted. In fact, recall that the basis functions of Table I   The element Jacobian (J ) and the gradient vectors (∇ξ a ) are expressed in terms of the unitary basis vectors ℓ 1 , ℓ 2 , ℓ 5 which are the derivatives of the element position vector r with respect to the independent coordinates ξ 1 , ξ 2 and ξ 5 . The element edges are formed by intersection of pairs of zero-coordinates surfaces, and the edge vectors are directed along the cross product of the associated coordinate gradients. The edges are given a two-index label deriving from the two coordinate indices appearing in this cross product. The unitary basis vectors determine the following eight edge vectors For example one gets Dependencies among higher order functions arise from linear combinations of the bases which contain one of the following identities as a factor (with Ω ab = −Ω ba ) The completeness relations are: which matches with that of the curl-conforming zeroth-order functions of the adjacent element having brick, tetrahedral or prismatic shape [11]. Once again, (12) is the sum of two polynomial functions: • The curl-free component is polynomial because it is the divergence of a polynomial that contains the factor (1 − ξ 5 ). • The non-zero curl component, clearly polynomial, has also polynomial curl since it contains the factor (1 − ξ 5 ). It is fascinating how our fractional functions of space ξ unfold into polynomials simply by mapping the pyramid to a cubic cell whose vertices are intersection points of 3 edges. As mentioned in Section II, it is reasonable to insist that this also applies to all other higher order functions, if only for uniformity and for energy considerations based on (11). We therefore state the following POSTULATE: In terms of the scaled variables the curl-conforming basis functions of the pyramid are polynomials with polynomial curl.
As postulated, the higher-order functions are obtained simply by multiplying the lowest order functions of Table I by polynomials of the space (η, ξ 5 ), and the order of a function is that of its multiplicative polynomial. The lowest order functions of Table I have zero order (in fact, on the boundary, they match the zero-order functions of adjacent elements). Expressions (7)-(13) wipe old habits [11] and allow to derive the pyramid's basis functions by shifting to a new paradigm: 1) The vector components of the basis functions and of the curl of the basis functions are polynomials of the variables η, ξ 5 . Unisolvency and base completeness must be proved in the space described by the grandparent variables. 2) Each higher order vector function is obtained by multiplying one vector function of zero order by a scalar generating polynomial which, in turn, is the product of normalized orthogonal polynomials (the same was done in [3], [4], [5], [6]).
3) The multiplicative polynomials are defined in the grandparent cubic cell of Fig. 1 (whose vertices are points of intersection of only three edges and faces). 4) On the pyramid border, the multiplicative polynomials that generate the edge and face based functions coincide with those for the adjacent elements, no matter what shape they have. According to this new paradigm, once the importance of the (1 − ξ 5 ) factor previously discussed is understood, we immediately derive the volume-based bubble functions and obtain all the others by simply imposing, on cell's boundary, the tangential continuity of the pyramid basis functions with the other known expressions valid for differently shaped cells. So keep in mind that it is only for the sake of brevity that conformity with differently shaped elements is "demonstrated" in Section V "after" the axiomatic definition of the basis functions.
Also note that the dependency relations (6,9) ensure that homogeneous polynomials of order p of the ξ-domain are polynomials of order p in ξ 5 of the (η, while, the other way around, the product of (1 − ξ 5 ) p times a polynomial of order p in η j is a polynomial of global order p in ξ j , ξ j+2 . In particular we get (16) where the coefficients b m depend on the coefficients a k , as well as on the value of p. Eq. (16) is useful to asses the global order of the multiplicative polynomials we introduce in Section V. Notice that (14) neglects the condition that ξ j and ξ j+1 are zero at the tip of the pyramid, while (15) takes automatically into account the fact that (14) must vanish at ξ 5 = 1. The same is true for (16) which vanishes at ξ 5 = 1 for p ≥ 1, a condition that we would have to impose if we used the expression on the right of (16).

V. HIGH-ORDER BASIS FUNCTIONS
The edge and face-based functions are derived by considering two groups: one consisting of functions associated with the mixed edges, and another associated with the triangle edges. Each function is the product of a zero-order vector function times a multiplicative scalar polynomial that we construct by using the normalized polynomials Q n (x, y) of Table II and the normalized Jacobi polynomials of Table III that we readily evaluate for any order using their appropriate recurrence relations [6], [20], [21]. We already used these polynomials to construct the hierarchical bases reported in [3], [5]; for example, the edge-based polynomials shown in the upper part of [3, Table II] coincide with Q n (ξ a , ξ b ) by where P n is the Legendre polynomial of order n. Properties: Recurrence relation: Referring back to (16), note that (6), (9) yield and that the polynomials of Table III are functions of (2z −1), with Eq. (7) and the first of (18) shows that, apart from the normalization constant, the polynomials (17) coincide with the shifted scaled Legendre polynomials used also by other authors [17], [22]. On the pyramid border, the edge and face-based functions match those given in [6] for the tetrahedron, the brick and the prism and this guarantees the mutual orthogonality of the edge-based functions on the edges and faces of the pyramid, as well as the mutual orthogonality of the face-based functions on the pyramid's faces. The edge and face-based functions are not mutually orthogonal for integrals on the volume of the pyramid. They can be rendered mutually orthogonal on the volume by adding to each of these functions an appropriate linear combination of the volume-based functions.

A. Volume-Based Functions
Volume-based vector functions have zero tangential component on all cell faces and are obtained by inspection observing that the gradients ∇η j = −∇η j+2 are orthogonal to both The polynomials P n (z) in the left column are obtained by rescaling the shifted Jacobi polynomials P faces η j = 0 and η j+2 = 0. The three fundamental bubblefunctions of the lowest order possible have a non-physical singular curl in ξ 5 = 1 and therefore do not belong to the polynomial base. In fact, the lowest order polynomial functions with polynomial curl are obtained simply by multiplying (20) by (1 − ξ 5 ), which immediately allows us to write all the higher-order polynomial bubble functions in the form In terms of the zeroth order vector functions of Table I, the bubbles (21) read as follows being ∇ξ 5 the sum of the functions Ω 12 , Ω 23 , Ω 34 , Ω 41 associated with the triangle edges, see Table I at bottom left. Before proceeding further, let us see the differences with respect to the bubble functions proposed by other authors, bearing in mind that our bubbles (21) contain the factor (1 − ξ 5 ) n with n = 1, which is the minimum value of n that guarantees polynomial curls. A factor (1 − ξ 5 ) n appears in the expression of the bubble functions of the fourth family of [17, eqs. (9.26), (B.35)] and in all those of [16], but in those expressions n depends on the maximum value between two of the indices i, j, k and can be greater than unity. Since the other basis functions have no factors that depend on the values of a pair of indexes, it is not clear to us why some bubbles should contain a factor (1 − ξ 5 ) n instead of the factor (1 − ξ 5 ) as in the case of our bubbles (21). By substituting the functions of the fourth family of [17] for our Ω B3 ijk , assuming they are equivalent, we have found the occurrence of spurious modes. We also observe that the expression of the higher order functions in [17] involve integrated Legendre polynomials, while we always use orthogonal Jacobi polynomials. It is certainly possible to use integrated Legendre polynomials instead of orthogonal Jacobi polynomials (normalized as in Table III); we have tried it and we have seen that by doing so the mass matrix CN increases by at least two orders of magnitude.
Our p-order complete curl-conforming set has a total of 3p 2 (p + 1) bubble functions since each Ω B1 , Ω B2 and Ω B3 generates p 2 (p + 1) functions hierarchically organized as follows where Ω B ijk stands for Ω B1 ijk , Ω B2 ijk , or Ω B3 ijk . Table IV reports the curls of the bubble functions in terms of the unitary basis vectors of Table I. For the maximum value p assumed by the indexes i, j, k, the component in ∇ξ 5 of the bubbles (21) is of order (p + 1) in η 1 , η 2 and ξ 5 . The ∇ξ 1 component of the bubbles (21) is of order p in η 1 , and (p + 1) in η 2 ; similarly, the ∇ξ 2 component is of order p in η 2 , and (p+1) in η 1 . The same happens for the vector functions associated with the triangle and mixed edges that we consider in the following two sub-sections. All this does not happen by accident, and it is used to demonstrate the completeness of the bases in the Appendix.

B. Functions Associated with the Mixed Edges
Higher order vector functions associated with the mixed edge γ5 (for γ = 1, 2, 3, 4) are obtained simply by multiplying the zeroth-order function by (p + 1)(3p + 2)/2 hierarchical polynomials obtained for k + i ≤ p. Set (25) has (p + 1) edge-based polynomials H γ5 00k that naturally form a hierarchical subset for k ranging from 0 to p, in the sense that for p = 0 the subset contains only the polynomial H γ5 000 , while for p ≥ 1 one has to increment the set of order (p − 1) with the polynomial H γ5 00p .
We then have p(p + 1) face-based polynomials H γ5 0jk that vanish on face γ (where η γ = 0) and are based on the ξ 5 = 0 face; they are hierarchically organized as follows • for p = 1 the set consists of the polynomials H γ5 010 , H γ5 011 ; • for p > 1 the set of order (p − 1) is augmented by adding p polynomials H γ5 0pk , obtained for k = 0, 1, . . . , p − 1, plus p polynomials H γ5 0jp obtained for j = 1, . . . , p. Finally we have p(p + 1)/2 face-based polynomials H γ5 i0k that vanish for ξ 5 = 0 and are different from zero on the ξ γ = 0 triangular face; that is, they are based on face γ and hierarchically organized as follows • for p = 1 the set consists of the polynomial H γ5 100 ; • for p > 1, one has to increment the set of order (p − 1) with p polynomials H γ5 i0k obtained for k = 0, 1, . . . , p − 1 with i = p − k (Note that the number of polynomials indicated here holds for any mixed edge and does not depend on the shape of the cell, be it a pyramid or a triangular prism [5], [6].) To show that the edge-based functions H γ5 00k Ω γ5 conform to the basis function of adjacent cells we observe that on the quadrilateral face ξ 5 = 0 the simplified expression matches that of the edge-based polynomials given in [5] for the prism, and in [4] for the brick while, to prove conformity on the triangular face η γ = 0, recall that and then observe that for (ξ γ+1 , ξ γ−1 ) = (ξ a , ξ b ), the polynomials (27) coincide with the edge-based polynomials reported in the upper part of [3, Table II] obtained for χ cd = ξ 5 = (1−ξ a −ξ b ). Equally simply it is proved that on the triangular face the multiplicative polynomials H γ5 00k coincide with those of the triangular prism given in [5, Table III].
As for the face-based polynomials H γ5 0jk we note that on the quadrilateral face ξ 5 = 0 we have η γ = ξ γ , so that we get which is the same expression holding for the face-based polynomials of the brick cell reported in [4, eq. (12)] (The Na and Nb values are obtained from Table II.) Lastly, and in a similar way, to demonstrate the conformity of the polynomials H γ5 iok based on the triangular face ξ γ = 0 with those of the adjacent cell (possibly of tetrahedral or prismatic shape), it is sufficient to compare their expression with those reported in [5, Table III] and [3, Table II].
As already done in [3], [4], [5], the sign of the vector functions (23) must be adjusted if the orientation of the mixed edge is opposite to the "reference" one, or different from that of the adjacent cells, as a consequence of (27) and of the symmetry relations of the Legendre polynomials
3) The gradient of (30) takes the polynomial expression That said, the higher order vector functions are simply obtained by multiplying the lowest-order function Ω γδ by (p + 1) 2 hierarchical polynomials which all have a gradient of polynomial form in the space (η, ξ 5 ). This implies that the curl of (32) takes a polynomial form in the space (η, ξ 5 ) since H γδ ijk , Ω γδ and its curl have polynomial form in this space. (The polynomial form of ∇Q k (ξ 5 , ξ γδ ) stems from (31).) Set (33) has (p + 1) edge-based polynomials H γδ 00k ; that is, for p = 0 this subset is simply formed by H γδ 000 while, for p ≥ 1, one has to increment the subset of order (p − 1) with the polynomial H γδ 00p . Set (33) then has p(p + 1)/2 polynomials H γδ 0jk based on the ξ δ = η δ = 0 triangular face, organized as follows • for p = 1 the set consists of the functions H γδ 010 ; • for p > 1, one has to increment the set of order (p − 1) with p functions H γδ 0jk obtained for k = 0, 1, . . . , p − 1 with j = p − k Finally, the p(p + 1)/2 polynomials H γδ i0k based on the ξ γ = η γ = 0 triangular face are hierarchically organized as follows • for p = 1 the set consists of the polynomial H γδ 100 ; • for p > 1, one has to increment the set of order (p − 1) with p polynomials H γδ i0k obtained for k = 0, 1, . . . , p − 1 with i = p − k The curls of the functions (32) are given in Table VI. The vector functions (32) naturally conform to the basis functions of adjacent cells. For example, along the γδ-edge one gets ξ γδ = 1 − ξ 5 and which proves that H γδ 00k conforms to the edge-based polynomials given in [3], [5], apart from a possible sign adjustment. Conformity on the triangular faces γ and δ descends from the fact that, on these faces, H γδ 00k takes the same form as (27) H γδ

FACE-BASED FUNCTIONS' CURL -With reference to (37) or (38) we first compute the gradients
Note that H γδ i0k and ∇H γδ i0k are obtained by exchanging γ with δ in the expressions that hold for H γδ 0jk and ∇H γδ 0jk . The edge vectors components of the curl of the functions (37) is then computed by using the following formulas for γ = 1, 2, 3, 4; δ = γ + 1 computed modulo 4, and with ∇ξ a × ∇ξ b = ℓ ab /J We recommend quadrature formulas with samples inside the element to avoid numerical problems due to the expressions of S a k , S b k provided above. Formulas to compute S a k and S b k for x and/or y equal to zero are however given below.
For k ≥ 1 use the following limits: For k ≥ 1 use the following limits: Similarly, the conformity of the face-based functions with those in [3], [5] is demonstrated by comparing their expression with those reported in [5, Table III] and [3, Table II].
In applications, the functions can be replaced by where, as indicated by the subscripts used are the zero-order curl-conforming functions of the twodimensional simplexes defined by the triangular faces ξ δ = 0 and ξ γ = 0, respectively; with At last, we clarify the order of the multiplicative polynomials (33) so to prove the completeness of our bases as in Appendix I. Q k (ξ 5 , ξ γδ ) is a polynomial of degree 3k of the (η γ , η δ , ξ 5 ) space whose highest degree term in ξ 5 is proportional to (ξ 5 η γ η δ ) k . Its higher degree term in η γ is proportional to [η γ (1 − η δ )(1 − ξ 5 )] k , therefore proportional to ξ k γ on face ξ δ = 0. The highest degree term in η δ is proportional to [η δ (1 − η γ )(1 − ξ 5 )] k ; that is, it is proportional to ξ k δ on the ξ γ = 0 face. Similarly, one can show that the highest degree term of all the polynomials (33) obtained using the hierarchical indexing rule reported above is proportional to (η γ η δ ξ 5 ) p , and that we never obtain terms proportional to η α γ , η α δ or ξ α 5 with α > p. Considering the expression (12) this is equivalent to saying that, for the maximum value assumed by the indexes i, j, k, the component in ∇ξ 5 of the functions (32) is of order (p+1) in η 1 , η 2 and ξ 5 . The ∇ξ 1 component is of order p in η 1 , and (p+1) in η 2 ; similarly, the ∇ξ 2 component is of order p in η 2 , and (p + 1) in η 1 . The same happens for the vector functions associated with the mixed edges considered in the previous sub-sections.

D. Degrees of Freedom
Not all the higher order face-based vector functions introduced so far are independent of each other. To form a pth order base, we must discard the dependent functions and count the total number of Degrees of Freedom (DoF). There are in all 8(p + 1) edge-based functions since the pyramid has 8 edges and all the edge-based functions are independent. Then, regarding the face-based functions, we observe that each lowest-order vector functions associated with the three edges that bound the triangular face ξ γ = 0 generates a set of p(p + 1)/2 face-based functions but one of these sets must be discarded because on a face we have only two independent tangent vectors, as per the dependency relationship in Table I. So, for p ≥ 1, on the ξ γ = 0 face we have a total of p(p + 1) face-based functions. Similarly, since the quadrilateral face is bounded by four mixed-edges, we have four different subsets of vector functions based on face ξ 5 = 0. Two of these are dependent and must be discarded as per the dependency relationship in Table I   • one component × (p + 1) DOF's × 8 edges = 8(p + 1) edge degrees of freedom; • two components × p(p + 1)/2 DOF's × four triangular faces = 4p(p + 1) face degrees of freedom; • two components × p(p + 1) DOF's × one quadrilateral face = 2p(p + 1) face degrees of freedom; • three components × p 2 (p + 1) DOF's × one element = 3p 2 (p + 1) degrees of freedom interior to a pyramid for a grand total of degrees of freedom per pyramid equal to DoF# = (p + 1)(8 + 6p + 3p 2 ) = 3 p 3 + 5 p as obtained in [17] for p = p + 1. The number of DoF of the pyramid is always lower than that of the brick, while it remains higher than that of the triangular prism for p ≥ 2 (see Fig. 2). Note that the pyramid and the brick have the same number of interior DoF [6].

VI. NUMERICAL RESULTS
The principal concern of hierarchical bases tends to be the matrix conditioning arising from their use (see for example [6], paragraph 5.3.4, from page 242 onwards). To establish that the rate of growth in condition number (CN) for our new pyramidal bases is not substantially worse than that of the hierarchical bases for the other differently shaped cells (bricks, triangular-prisms, and tetrahedrons) we show below a couple of illustrative results obtained with rectilinear elements. We do this for the sake of brevity as our bases can also be used on curved pyramids appearing in hybrid meshes. To better asses (and then predict) the CN behavior one should consider and study many hybrid meshes that use differently shaped elements. Fig. 3 reports results for the individual element mass-matrix condition numbers for hierarchical vector bases of different order obtained by considering rectilinear cells with equal edges and of unitary length, computed with 3D-integrals over the child cells. (Note that the CNs shown in Fig. 3 do not depend on the cells edge-length in the child space since the unitary basis vectors and the Jacobian of the transformation from parent to child space are constant for the cells considered by Fig. 3.) The CNs of the pyramid shown in Fig. 3 are obtained by using for the triangular faces the basis functions associated with the triangle edges. Fig. 4 compares the individual element mass-matrix condition numbers of the equilateral pyramid of Fig. 3 with those for pyramids obtained by moving one vertex of the base of the equilateral pyramid along its diagonal, doubling and tripling the length of this diagonal as depicted in Fig. 5. These pyramids have equal height and equal length for one of the diagonals of their base, but flat quadrilateral base of different shape. The ratio between the longest and the shortest side of each cell, commonly known as Aspect Ratio (AR), is given in the captions of Fig. 4 and 5. Note that unlike the equilateral pyramid, the distorted pyramids considered in Fig. 4 do not have a constant Jacobian. In view of the (expected) results of Fig. 4, we recommend using cells with AR near unity and less than 3 when using bases of order higher than the first.  13 on the right. The Jacobian J of the transformation from parent-to-child space is constant for an equilateral pyramid, while for distorted pyramids it can vary within the cell. For example, we have J = K for the pyramid shown on the left, while J = K(1 + η 1 + η 2 ) for the pyramid shown in the center, and J = K(1 + 2η 1 + 2η 2 ) for the pyramid shown on the right. equation. Results obtained using a single pyramidal cell are compared with those obtained using four identical tetrahedral cells. These results not only demonstrate that our pyramid bases avoid spurious modes, but are also able of providing better results with fewer unknowns than that required by a tetrahedral cell-based model. To reduce the number of degrees of freedom of the tetrahedral cell-based model we also tried with a mesh formed by only two identical tetrahedral elements instead of four (these results are not reported here), and we have seen that this heavily impacts on the symmetry of the cavity modes and leads to errors much higher than those reported in the lower part of Table VII. VII. CONCLUSIONS This paper presents a general procedure to obtain higher order hierarchical curl-conforming vector basis functions for pyramidal elements. The functions ensure the continuity of the proper vector components across adjacent elements of equal order but different shape. Properties of the vector basis functions are discussed in detail. The reported numerical examples show that higher order functions provide more accurate results than those obtainable with low-order elements.