Automated generation and symbolic manipulation of tensor product finite elements

We have extended the domain-specific language UFL, the Unified Form Language, to support symbolic product operations on finite elements. This allows representation of scalar- and vector-valued finite element spaces on quadrilateral, triangular prism and hexahedral cells. In particular, we concentrate on spaces that correspond to scalar and vector identifications of tensor product finite element differential forms. This includes fully continuous, curl-conforming, div-conforming, and discontinuous spaces in both two and three dimensions. We have made corresponding extensions to FIAT, the FInite element Automatic Tabulator, to enable numerical tabulation of these spaces. This facilitates the automatic generation of low-level code that carries out local assembly operations, within the wider context of solving finite element problems posed over such function spaces. We have done this work within the code-generation pipeline of the software package Firedrake; we make use of the full Firedrake package in this paper to present numerical examples.


Introduction
Many different areas of science benefit from the ability to generate approximate numerical solutions to partial differential equations. Traditionally, the implementation of suitable numerical schemes has been a labour-intensive task, often requiring a bespoke program written in Fortran, C++, or some other low-level language. In stark contrast, the FEniCS Project (Logg et al., 2012a) allows a user to express discretisations, based on the finite element method, in UFL (Alnaes et al., 2014;Alnaes, 2012): a concise, high-level language embedded in Python. The corresponding "kernels", containing low-level code, are automatically generated by FFC, the FEniCS Form Compiler (Kirby and Logg, 2006;Logg et al., 2012b), making use of FIAT (Kirby, 2004(Kirby, , 2012. These local assembly kernels are executed on each cell 1 in the mesh, and the resulting linear systems can be solved using a number of third-party libraries. This has the overall effect of reducing development time and increasing productivity: the automated code generation process is extensively tested, and significantly fewer lines of code need to be written by the end user, reducing the likelihood of serious bugs.
Firedrake is an alternative software package which presents a similar -in many cases, identical -interface to that of FEniCS. Firedrake reuses components of the FEniCS Project in order to automatically generate low-level kernels, but the execution of these kernels over the mesh is performed in a fundamentally different way. This leads to significant performance increases relative to FEniCS 1.4, across a range of different problems . Firedrake and FEniCS have other distinguishing features: for example, they support a wide range of arbitrary-order finite element families. These may be utilised by numerical analysts proposing novel discretisations of existing PDEs. However, a limitation of Firedrake and FEniCS has been the lack of support for anything other than unstructured meshes consisting of simplicial cells: intervals, triangles or tetrahedra.
There are several reasons why a user may wish to use a mesh of non-simplicial cells. Our main motivation is geophysical simulations. These are often governed by highly anisotropic equations in which gravity plays an important role. In addition, they often require high aspectratio domains: the vertical height of the domain may be several orders of magnitude smaller than the horizontal width. These domains allow a decomposition which has an unstructured horizontal 'base mesh' but has regular vertical layers -we will refer to this as an extruded mesh. From a mathematical viewpoint, the vertical alignment of cells minimises difficulties associated with the anisotropy of the governing equations. From a computational viewpoint, the vertical structure can be exploited to improve the performance relative to a generic unstructured mesh.
There are a few ways to generate an orderly hierarchy of cells in an increasing number of spatial dimensions. One method is to add a single new point to an existing lower-dimensional cell; this leads to the family of simplices: points (0D), intervals (1D), triangles (2D), tetrahedra (3D) and so on. Another method is to take a geometric product of intervals, say [0, 1] n , giving points (0D), intervals (1D), quadrilaterals (2D), hexahedra (3D) and so on. This can be generalised by allowing products of any set of simplices, not just intervals. There are other related cells such as the square-based pyramid, which is formed by adding a single new point to a product of intervals. However, we will only consider product cells in this paper. In two and three dimensions, this gives rise to three non-simplicial cells: quadrilaterals, triangular prisms and hexahedra.
Since our cells have a product structure, we will focus on producing finite element families that can be expressed as (sums of) products of existing finite elements. This covers many, though certainly not all, of the important finite element spaces on product cells. We pay special attention to element families relevant to finite element exterior calculus, a mathematical framework that leads to stable mixed finite element discretisations of partial differential equations (Arnold et al., 2006(Arnold et al., , 2010. This paper hence describes the extension of the Firedrake code-generation pipeline to enable support for finite element discretisations on cells which are products of simplices. This enables the automated generation of low-level kernels representing finite element operations on such cells. We remark that, due to our geophysical motivations, Firedrake has complete support for extruded meshes whose base mesh is built from simplices. At the time of writing, however, it does not support generic unstructured quadrilateral, prismatic or hexahedral meshes. This paper is structured as follows: in section 2, we provide the mathematical background for product finite elements with minimal reference to implementation details. In section 3, we describe the software extensions in detail. In section 4, we present several numerical experiments using the extruded mesh functionality of Firedrake to demonstrate the correctness of our implementation. Finally, in section 5 and section 6, we give some limitations and other closing remarks.

Mathematical preliminaries
This section is structured as follows: in subsection 2.1, we give the definition of a finite element. In subsection 2.2, we briefly define the sum of finite elements. In subsection 2.3, we describe the various types of inter-cell continuity of the finite element spaces that we use. In subsection 2.4, we give examples of finite element spaces on simplicial cells compatible with finite element exterior calculus. In subsection 2.5 and subsection 2.6, which form the main part of this section, we define the product of finite elements and state how these products can be manipulated and combined to produce elements compatible with finite element exterior calculus. In subsection 2.7, we make some remarks on numerical quadrature rules in product cells. Finally, in subsection 2.8, we re-state subsection 2.5 and subsection 2.6 in the language of differential forms, which provides a far more natural setting for the underlying operations.

Definition of a finite element
We will follow Ciarlet (1978) in defining a finite element to be a triple (K, P, N) where • K is a bounded domain in R n , to be interpreted as a generic reference cell on which all calculations are performed, • P is a finite-dimensional space of continuous functions on K, • N = {n 1 , . . . , n dim P } is a basis for the dual space P -the space of linear functionals on P -where the elements of the set N are called nodes.
Let Ω be a compact domain which is decomposed into a finite number of non-overlapping cells. Assume that we wish to find an approximate solution to some partial differential equation, posed in Ω, using the finite element method. A finite element together with a given decomposition of Ω produce a finite element space.
A finite element space is a finite-dimensional function space on Ω that is used in calculations; for example, the approximate solution will live in such a space. Essentially, there are two things that need to be specified to characterise a finite element space: the form that a function may take within a single cell, and the amount of continuity required between neighbouring cells.
The former is given by, in the simplest cases, the composition of P with (the inverse of) some map from the reference cell to a given physical cell. Almost all finite elements take P to be some subspace of polynomials, often simply all polynomials up to a given degree. A basis for P is very useful in implementations of the finite element method. A common choice is a nodal basis in which each of the basis functions Φ 1 , . . . , Φ dim P vanish at all but one node: (2.1) in practice, this condition may be slightly relaxed. Basis functions from different cells can be combined into basis functions for the finite element space on Ω. The inter-cell continuity of these basis functions is related to the choice of nodes, which is the core topic of subsection 2.3.

Sum of finite elements
Suppose we have two finite elements U = (K, P A , N A ) and V = (K, P B , N B ) defined over the same reference cell K. If the intersection of P A and P B is trivial, we can define the direct sum U ⊕V to be the finite element (K, P, N), where Even though the intersection of P A and P B is trivial, it is not generally true that concatenating the existing nodal bases for P A and P B gives a nodal basis for P. However it is not always necessary to use a strictly nodal basis and simply concatenating the existing bases may be sufficient within an implementation.

Sobolev spaces, inter-cell continuity, and Piola transforms 2.3.1 Sobolev spaces
Sobolev spaces are used to characterise the smoothness of solutions to partial differential equations. One of the strengths of the finite element method is that the approximate solution lies in some finite-dimensional subspace of a larger Sobolev space containing the exact solution. Sobolev space notation is, therefore, often invoked to characterise the smoothness of functions in finite element spaces -in particular, the inter-cell continuity.
We start by defining L 2 (Ω) to be the space of square-integrable functions on Ω. We say that f is in the Sobolev space H n (Ω) if all weak derivatives up to order n exist and are in L 2 (Ω). In particular, the Sobolev space H 1 (Ω) requires the function and all weak first derivatives to exist and be square-integrable; this is sometimes written as where ∇ f should be interpreted in a weak sense.
We remark that square-integrability is guaranteed, in the context of finite element spaces, since the domain is compact and the functions are polynomials on each cell of Ω. On the other hand, the existence of weak derivatives acts as a genuine constraint. If f is in a finite element space, the statement f ∈ H 1 (Ω) is equivalent to f being continuous between cells. Alternatively, one can say that f is single-valued at any vertex, edge or facet 2 between neighbouring cells.
Given a function u in some vector-valued finite element space, we can also say u ∈ H 1 (Ω) (resp. L 2 (Ω)) if each scalar component is in H 1 (Ω) (resp. L 2 (Ω)), as defined above. However, vector fields admit a more fine-grained treatment of inter-cell continuity than scalar fields. For example, we say that u ∈ H(div; Ω) if both u and its weak divergence exists and are square integrable: (2.5) As before, we note that the important condition is the existence of the weak divergence of u.
If u is in a (vector-valued) finite element space, the statement u ∈ H(div; Ω) is equivalent to the normal component of u being single-valued at any facet between a pair of neighbouring cells.
We can define H(curl; Ω) similarly. In three dimensions, we demand that both u and its weak curl exist and are square-integrable: (2.6) In two dimensions, we interpret curl as the scalar-valued operator, known sometimes as rot. For u in a vector-valued finite element space, the statement u ∈ H(curl; Ω) is equivalent to the tangential component of u being single-valued at any edge or facet between cells. Further descriptions of these Sobolev spaces can be found in Boffi et al. (2013).
We make two brief remarks before finishing this subsection. Firstly, we will generally suppress the domain Ω from our notation: we will simply write H 1 , H(curl), H(div), and L 2 . Secondly, it is clear that the Sobolev spaces mentioned above have some trivial inclusion relations: H 1 , H(div), and H(curl) are all subspaces of L 2 , and vector-valued functions in H 1 are in both H(div) and H(curl). However, when we make casual statements such as V ⊂ H(div), it is implied that V ⊂ H 1 , i.e., we have made the strongest statement possible. In particular, we will use L 2 to denote a total absence of continuity between cells.

Inter-cell continuity and the geometric decomposition of nodes
The set of nodes N, from the definition in subsection 2.1, are used to enforce the continuity requirements on the 'global' finite element space. This is done by associating nodes with topological entities of K -vertices, facets, and so on. When multiple cells in Ω share a topological entity, the cells must agree on the value of any degree of freedom associated with that entity. For example, a degree of freedom on a facet leads to coupling between any cells that share a facet, while a degree of freedom on a vertex leads to coupling between all cells adjacent to a single vertex. In the context of implementation, there may be both inter-and intra-entity symmetry constraints on the placement of nodes, particularly on an unstructured mesh.
In many communities, a node is merely synonymous with function evaluation at a given point in K. However, the more general definition, given in subsection 2.1, admits a much wider range of nodes. One common example is function moments against some space of polynomials. Another example is function evaluation influenced by geometric properties of K, such as the component of a vector-valued function normal to a facet. The exact choice of nodes is often not so important in determining the finite element space. However, the association of nodes with topological entities is crucial -this is sometimes called the geometric decomposition of nodes.
For H 1 elements, functions must be single-valued on vertices, edges and facets: 1. Nodes are firstly associated with vertices.
2. If the function is not yet uniquely specified in the entire cell (i.e., if the dimension of the polynomial space, P, exceeds the number of nodes on vertices), additional nodes are associated with edges until the function is specified uniquely on edges.
3. If the function is not yet uniquely specified in the entire cell (i.e., if the dimension of P exceeds the number of nodes on vertices and edges), additional nodes are associated with facets until the function is specified uniquely on facets.
4. If the function is not yet uniquely specified in the entire cell (i.e., if the dimension of P exceeds the number of nodes on vertices, edges and facets), any remaining nodes are associated with the interior of the cell.
For H(curl) elements, which are intrinsically vector-valued, the component(s) of the function tangential to edges and facets must be single-valued: 1. Nodes are firstly associated with edges until the tangential component of the function is specified uniquely on edges.
2. If the function is not yet uniquely specified in the entire cell (i.e., if the dimension of P exceeds the number of nodes on edges), additional nodes are associated with facets until the tangential components of the function are specified uniquely on facets.
3. If the function is not yet uniquely specified in the entire cell (i.e., if the dimension of P exceeds the number of nodes on edges and facets), any remaining nodes are associated with the interior of the cell.
For H(div) elements, which are also intrinsically vector-valued, the component of the function normal to facets must be single-valued: 1. Nodes are firstly associated with facets until the normal component of the function is specified uniquely on facets.
2. If the function is not yet uniquely specified in the entire cell (i.e., if the dimension of P exceeds the number of nodes on facets), any remaining nodes are associated with the interior of the cell.
L 2 elements have no continuity requirements: 1. All nodes are associated with the interior of the cell. Since the interior of the cell is not shared between multiple cells, this does not lead to any continuity constraints on the global mesh.
Note that some element families do not strictly follow the above -for example, the lowestorder Crouzeix-Raviart element is only in L 2 , but is continuous at the midpoint of edges. It therefore has nodes associated with edges, but not enough to specify the function uniquely on edges. However, as such elements do not fit into the framework of finite element exterior calculus, we shall not consider them further.

Piola transforms
For functions to have the desired amount of continuity on the global mesh, it is important that they undergo an appropriate mapping from the reference cell to the physical cell. Let X represent coordinates on the reference cell, and x represent coordinates on the physical cell. Each physical cell can be represented by some (invertible) map from the reference cell, say x = φ ( X).
To map an H 1 or L 2 function from the reference cell to the physical cell, the simplest mapping suffices. Letf ( X) be a function defined over the reference cell. The corresponding function f ( x) defined over the physical cell is then (2.7) We will refer to this as the identity mapping.
However, we may wish to have continuity of the normal or tangential component of the vector field in physical space; (2.7) clearly does not preserve this. H(div) and H(curl) elements therefore use Piola transforms when mapping functions from reference space to physical space. We will use J to denote Dφ ( X), the Jacobian of the coordinate transformation. H(div) functions are mapped using the contravariant Piola transform, which preserves normal components: while H(curl) functions are mapped using the covariant Piola transform, which preserves tangential components: Further details on the Piola transform can be found in, for example, Boffi et al. (2013).
Another property of the Piola transforms is that they map between the tangent spaces of the reference cell and the physical cell. This is particularly important when solving an equation on, for example, a triangulation of the surface of the sphere. In such cases, the Jacobian matrix J is not square, and pseudodeterminants or pseudoinverses must be used in (2.8) and (2.9). More details can be found in, for example, Rognes et al. (2013).

Simplicial finite element spaces within finite element exterior calculus
The work of Arnold et al. (2006Arnold et al. ( , 2010 on finite element exterior calculus provides principles for obtaining stable mixed finite element discretisations on a domain consisting of simplicial cells: intervals, triangles, tetrahedra, and higher-dimensional analogues. In full generality, this involves de Rham complexes of polynomial-valued finite element differential forms linked by the exterior derivative operator. In 1, 2 and 3 dimensions, differential forms can be naturally identified with scalar and vector fields, while the exterior derivative can be interpreted as a standard differential operator such as grad, curl, or div. The vector-valued element spaces are slightly unusual since they only have partial continuity between cells: they are in H(curl) or H(div), which have been discussed already. The element spaces themselves were, however, already well-known in the existing finite element literature for their use in solving mixed formulations of the Poisson equation and problems of a similar ilk.
We give some explicit examples of de Rham complexes of finite element spaces, compatible with finite element exterior calculus, on intervals, triangles, and tetrahedra. When specifying finite element spaces, we have used the notation from the Periodic Table of Finite Elements (Arnold and Logg, 2014), which is compatible with the notation used in UFL.
On an interval (one spatial dimension), the standard example is (2.10) where P n represents the space of functions which are piecewise polynomials of degree n, continuous between cells, and DP n represents the space of functions which are piecewise polynomials of degree n, possibly discontinuous between cells -this is often denoted by P DG n in other literature. The P spaces, being continuous between cells, are subspaces of H 1 , while the DP spaces are only in L 2 . In one dimension, the exterior derivative corresponds to the usual derivative operator.
On triangles (two spatial dimensions), the examples given in Arnold et al. (2010) where RTF n represents the Raviart-Thomas element (Raviart and Thomas, 1977) and BDMF n represents the Brezzi-Douglas-Marini element (Brezzi et al., 1985). Both of these vector-valued element families are partially continuous between cells: across a facet (edge), the component normal to the facet must be continuous, but the tangential component can be discontinuous. These families are therefore in H(div). Here, the exterior derivative corresponds to ∇ ⊥ for the scalar P n element, a pointwise-rotated version of the usual 2D gradient operator. For the vector element, the exterior derivative corresponds to the usual 2D divergence operator here. Crucially, the identity ∇ · ∇ ⊥ = 0 mimics the exterior derivative identity d 2 = 0.
Another example, highlighted in Cotter and Shipton (2012) as being relevant for numerical weather prediction, is where B 3 represents a cubic 'bubble' function, vanishing at cell edges, and BDFM 2 represents a member of the vector-valued Brezzi-Douglas-Fortin-Marini element family (Brezzi and Fortin, 1991), also in H(div).
Before continuing, we remark that there are two main conventions in the literature for the numbering of certain element families. We have used the convention in which the lowestorder RTF element, containing some -but not all -linear polynomials, is labelled RTF 1 , rather than RTF 0 . Similarly, the BDFM 2 element contains some -but not all -quadratic polynomials.
In two dimensions, there is a second class of examples arising from an alternative identification of differential 1-forms with vector fields. For the examples given in Arnold et al. (2010), taking the alternative identification leads to the following: (2.14) and P n ∇ −→ BDME n−1 (2.15) the RTE and BDME families are essentially pointwise-rotated versions of the RTF and BDMF families. They therefore have continuous tangential component across the edges between cells, while the component normal to the edge may be discontinuous. These families are in H(curl) rather than H(div). When using this alternative identification of 1-forms with vector fields, the exterior derivative corresponds to the usual 2D gradient operator for the scalar P n element, and ∇ ⊥ · for the vector element, essentiallyẑ · ∇×. As before, we have the identity ∇ ⊥ · ∇ = 0.
On tetrahedra (three spatial dimensions), the examples given by Arnold et al. correspond to where N1E and N2E denote the Nédélec edge elements of the 1st and 2nd kind, and N1F and N2F denote the Nédélec face elements of the 1st and 2nd kind. These element families were introduced in Nédélec (1980Nédélec ( , 1986. The edge elements provide subspaces of H(curl), while the face elements provide subspaces of H(div). These spaces on tetrahedra will be less important to us since we will mostly consider taking the product of spaces on intervals and triangles. In three dimensions, the exterior derivative d now corresponds to each of the standard differential operators in turn: grad, curl, and div. The identity d 2 = 0 is encoded in the identities ∇ × ∇ = 0 and ∇ · ∇× = 0.

Product finite elements
In this section, we discuss how to take the product of a pair of finite elements and how this product element may be manipulated to give different types of inter-cell continuity. We will label our constituent elements U and V , where U := (K A , P A , N A ) and V := (K B , P B , N B ) following the notation of subsection 2.1. We begin with the definition of the product reference cell, which is straightforward. However, the spaces of functions and the associated nodes are intimately related, hence the discussion of these is interleaved.

Product cells
Given reference cells K A ⊂ R n and K B ⊂ R m , the reference product cell K A × K B can be defined straightforwardly as follows: (2.20) The topological entities of K A × K B correspond to products of topological entities of K A and K B . If we label the entities of a reference cell (in R n , say) by their dimension, so that 0 corresponds to vertices, 1 to edges, . . . , n − 1 to facets and n to the cell, the entities of K A × K B can be labelled as follows:  It is important to distinguish between different types of entities, even those with the same dimension. For example, if K A is a triangle and K B an interval, the (2, 0) facets of the prism K A × K B are triangles while the (1, 1) facets are quadrilaterals.

Product spaces of functions -simple elements
Given spaces of functions P A and P B , the product space P A ⊗ P B can be defined as the span of products of functions in P A and P B : where the product function f · g is defined so that In the cases we consider explicitly, at least one of f or g will be scalar-valued, so the product on the right-hand side of (2.22) is unambiguous. A basis for P A ⊗ P B can be constructed from bases for P A and P B . If P A and P B have nodal bases respectively, a nodal basis for P A ⊗ P B is given by (2.25) the right-hand side uses the same product as (2.22).
While this already gives plenty of flexibility, there are cases in which a different, more natural, space can be built by further manipulation of P A ⊗ P B . We will return to this after a brief description of product nodes.

Product nodes -geometric decomposition
Recall that the nodes are a basis for the dual space (P A ⊗ P B ) , and that the inter-cell continuity of the finite element space is related to the association of nodes with topological entities of the reference cell.
Assuming that we know bases for P A and P B , there is a natural basis for (P A ⊗ P B ) which is essentially an outer product of the bases for P A and P B . Let n i, j denote a "product" of n j , the j'th node in N B -though we will refrain from actually defining the product of nodes until subsubsection 2.5.5. If n (A) i is associated with an entity of K A of dimension p and n (A) j is associated with an entity of K B of dimension q then n i, j is associated with an entity of K A × K B with label (p, q).
This geometric decomposition of nodes in the product element is used to motivate further manipulation of P A ⊗ P B to produce a more natural space of functions, particularly in the case of vector-valued elements.

2D
We take the reference cells K A and K B to be intervals, so the product cell K A × K B is twodimensional. The elements on intervals are scalar-valued and are either in H 1 or L 2 . There are, therefore, three cases that we will consider: the product of two H 1 elements, the product of two L 2 elements, and the product of an H 1 element with an L 2 element (or vice-versa). A summary of the below is given in Table 2.1.
The vertices of the constituent intervals are guaranteed to have nodes associated with them. Taking the product guarantees that there are nodes associated with the vertices of the product cell -with label (0, 0). This suggests that the product element is in H 1 , i.e., fully continuous. It can be further verified that, due to the product structure of both the nodes and the basis functions, the nodes on the facets with labels (0, 1) and (1, 0) are sufficient to uniquely determine the function on the facets. The functions are, therefore, single-valued on facets, as required. The constituent elements are scalar-valued, and so the product is scalar-valued; no further manipulation is required.
The H 1 element is guaranteed to have nodes associated with the vertices of its reference cell, while the L 2 element only has nodes associated with the interior. There are therefore no nodes associated with the vertices of the product cell, but there are nodes associated with the (0, 1) facets. It is clear that the product element does not have enough continuity to be in H 1 , but yet it seems to have 'too much' continuity to only be in L 2 . In particular, it is continuous across some facets but not across others. Recall that, in two dimensions, the vector-valued H(div) and H(curl) elements have nodes associated with facets, but not with vertices. It therefore seems plausible that we can interpret the product as an H(div) or H(curl) element. The constituent elements are both scalar-valued, and so their product is scalar-valued. To make a vector-valued element, the scalar needs to be padded out into a two-dimensional vector, as we have a two-dimensional reference cell.

Product (1D × 1D) Components Modifier Result
Mapping Table 2.1: Summary of 2D product elements If we interpret the scalar as the first component of the vector, this gives us an H(div) element. If we take the scalar to be the second component of the vector, we get an H(curl) element. If we had instead taken the product of an L 2 element with an H 1 element, the reverse would be true. In all these cases, we would now need to use an appropriate Piola transform when mapping functions from reference space into physical space. It can be easily verified that, due to the product structure of the basis functions and nodes, there are indeed enough nodes on the facet to guarantee continuity of the normal or tangential component. Note that the scalar-valued product element is a perfectly legitimate finite element, and it is not compulsory to manipulate it into a vector-valued element. However, the vector-valued elements are more useful and fit more naturally within Finite Element Exterior Calculus, as we will see in subsection 2.6.
L 2 ×L 2 : The constituent intervals only have nodes associated with their interior, and so the product element only has nodes associated with the interior of the product cell, with label (1, 1). This does not give rise to any continuity constraints, and hence the product element leads to a fully discontinuous finite space over the global mesh. The product is again naturally scalar-valued, and no further manipulation is required.

3D
We now consider the case where K A ⊂ R 2 and K B is an interval, so the product cell K A × K B is three-dimensional. The elements on a two-dimensional reference cell may be either scalar-or vector-valued, and may be in H 1 , H(curl), H(div) or L 2 . The elements on a one-dimensional reference cell, as before, are scalar-valued and are either in H 1 or L 2 . This leads to eight distinct cases. However, in two dimensions, H(curl) and H(div) elements are closely related: either can be transformed into the other by a 90 degree pointwise rotation of values. This reduces the number of cases to six, which we will cover below. A summary is given in Table 2.2.
Note: In the pictures below, we have taken the two-dimensional cell to be a triangle. However, the discussion is equally valid for quadrilaterals.
H 1 × H 1 : As in the two-dimensional case, the vertices of the constituent reference cells are guaranteed to have nodes associated with them, and so the product cell also has nodes associated with its vertices. Similarly, the edges and facets of the product cell have a sufficient number of nodes in order to uniquely determine the function on the interfaces between cells. The product element is again in H 1 , i.e. fully continuous. The constituent elements are scalar-valued, and so the product is naturally scalar-valued; no further manipulation is required.
H 1 × L 2 : As in the two-dimensional case, the H 1 element is guaranteed to have nodes associated with the vertices of its reference cell, while the L 2 element on an interval only has nodes associated with its interior. There are, therefore, no nodes associated with the vertices of the product cell, but there are nodes associated with the (0, 1) edges. The product is naturally scalar-valued, but the element appears to have some degree of continuity without being in H 1 .
In three dimensions, it is H(curl) elements that have nodes associated with edges but not with vertices. If we pad the scalar-valued product into a three-dimensional vector whose first two components are zero, this indeed produces an H(curl) element which requires an appropriate Piola transform. It can be verified that there are a sufficient number of nodes on (0, 1) edges and (1, 1) facets to specify the tangential component of the function uniquely.
H(div)/H(curl) × H 1 : Recall that both H(div) and H(curl) elements in two dimensions have nodes associated with edges (facets) of the cell. H 1 elements on an interval have nodes associated with vertices. The product therefore has nodes associated with the (1, 0) edges of the product cell, but no nodes are associated with the vertices. This again suggests interpreting the product as an H(curl) element. The product naturally takes values in R 2 , since the two-dimensional element is vector-valued and the one-dimensional element is scalarvalued. However, an H(curl) element in three dimensions must take values in R 3 , so the product needs padding, at minimum. If the two-dimensional element is in H(curl), it is sufficient to pad the product into a three-dimensional vector whose third component is zero. If the two-dimensional element is in H(div), padding alone does not give an H(curl) product element in three dimensions; this would give an element whose normal and tangential components are continuous across the (2, 0) facets but whose normal component is continuous across the (1, 1) facets.
Instead, the two-dimensional product must be rotated by 90 degrees before being padded into a three-dimensional vector. In both cases, an appropriate Piola transform is required, and it can be verified that there are a sufficient number of nodes on (1, 0) edges, (1, 1) facets and (2, 0) facets to specify the tangential component of the function uniquely.
L 2 × H 1 : In this case, nodes are only associated with the interior of the two-dimensional element. The one-dimensional element, on the other hand, has nodes associated with the vertices of the interval. The product therefore has nodes associated with the (2, 0) facets, and possibly with the interior. Recall that, in three dimensions, H(div) elements do not have nodes associated with vertices or edges, but do have nodes associated with facets. This element can therefore be interpreted as an H(div) element.
The product is naturally scalar-valued. As before, this is a perfectly legitimate finite element in its own right, but it is more useful to pad the scalar-valued product into a vector-valued element whose first two components vanish. This then gives a three-dimensional H(div) element, which requires an appropriate Piola transform. It can be verified that there are a sufficient number of nodes on the (2, 0) facets to specify the normal component of the function uniquely.
H(div)/H(curl) × L 2 : The two-dimensional element has nodes associated with the edges/facets of the two-dimensional cell. The one-dimensional element, on the other hand, only has nodes associated with the interior of the interval. The product therefore has nodes associated with the (1, 1) facets, and possibly with the interior, but not with any vertices or edges. This suggests interpreting the product as a three-dimensional H(div) element. Again, the product naturally takes values in R 2 , and needs to be padded. If the two-dimensional element is in H(div), it is sufficient to pad the product into a three-dimensional vector-valued element whose third component vanishes. If the two-dimensional element is in H(curl), padding alone would give a product element whose tangential component is continuous across (1, 1) facets, but whose normal component across (2, 0) facets vanishes (and is hence continuous). Again, the product must be rotated by 90 degrees before padding; this gives an H(div) element in three-dimensions. In both cases, an appropriate Piola transform is required, and it can be verified that there are a sufficient number of nodes on the (1, 1) facets to specify the normal component of the function uniquely.
L 2 × L 2 : As before, all nodes are associated with the interior of the constituent reference cells, and so all the nodes of the product element are associated with the interior of the product cell. The product element is naturally discontinuous and scalar-valued, and no further manipulation is required.

Product nodes -2D and 3D
This section covers the definition of the nodes of the product element. As mentioned previously, historically there have been a wide range of possible nodes. Here we will focus on just a few, but enough to construct all elements mentioned in the paper.
The simplest nodes are evaluation of a scalar-valued function at a given point, and evaluation of a component of a vector-valued function at a given point. More complicated nodes are evaluation of the normal or tangential component of a vector-valued function on an edge or facetused in H(div) and H(curl) elements. As in the previous subsection, the most effective way to describe the process is through an exhaustive list.

2D
The reference cells K A and K B are intervals; the elements on intervals are scalar-valued and are Product (2D × 1D) Components Modifier

Result
Mapping Table 2.2: Summary of 3D product elements. The elements marked with † are of little practical use; they are 2-vector valued, but the tangent space to the cell is 3-dimensional. No mapping has been given for these elements; the Piola transformations from a 3D cell require all three components to be defined.
either in H 1 or L 2 . All nodes are simply point evaluation.

3D
The reference cell K A ⊂ R 2 , while K B is an interval. If the element on K A is in H 1 or L 2 , the nodes are point evaluation. If the element is in H(curl) (resp. H(div)), the nodes on the facets are evaluation of the tangential (resp. normal) component of the vector-valued function, while the interior nodes are evaluation of the two components of the function at given points. The elements on K B are scalar-valued, and so all nodes are point evaluation.
H 1 × H 1 : The nodes of the product element are simply point evaluation at the set of points which is the Cartesian product of the sets of points of the 2D and 1D element nodes. L 2 × L 2 : The nodes of the product element are simply point evaluation at the set of points which is the Cartesian product of the sets of points of the 2D and 1D element nodes.

Consequences for implementation
It is clear that we need an operator that takes the product of two existing elements; this will be called OuterProductElement. This will generate a new element whose reference cell is the product of the reference cells of the constituent elements, as described in subsubsection 2.5.1. It will also construct the product space of functions P A ⊗ P B , as described in subsubsection 2.5.2, but with no manipulation (e.g., padding with zeros). The basis for P A ⊗ P B is as defined in (2.24) and (2.25). The nodes are topologically associated with topological entities of the reference cell as described in subsubsection 2.5.3.
To construct the more complicated vector-valued finite elements, we introduce additional operators HCurl and HDiv which form a vector-valued H(curl) or H(div) element from an existing OuterProductElement. This will modify the product space as described in subsubsection 2.5.4 by padding the existing product with zeros (after rotation, if applicable). These operators will also modify the nodes, as described in subsubsection 2.5.5, and set an appropriate Piola transform to be used when mapping functions. The basis is still of the form given in (2.24) and (2.25), but the values must be manipulated in the same way as the product space.  generalises the results of finite element exterior calculus on simplices to cells which can be expressed as geometric products of simplices. It also describes the construction of a specific complex of finite element spaces on hexahedrons (and, implicitly, quadrilaterals). When the differential forms are identified with scalar-and vector-valued functions, these correspond to well-known spaces such as the scalar-valued Q r and Q DG r , and various well-known vector-valued spaces as introduced in Brezzi et al. (1985); Nédélec (1980Nédélec ( , 1986. Within finite element exterior calculus, there are element spaces of interest which cannot be expressed as a tensor product of spaces on simplices -see, for example, Arnold and Awanou (2014) -but we shall not consider such spaces in this paper.

Product finite elements within finite element exterior calculus
Recall from subsection 2.4 that finite element exterior calculus makes use of de Rham complexes of finite element spaces. In one dimension, the complex takes the form where U 0 ⊂ H 1 and U 1 ⊂ L 2 .
In two dimensions, there were two types of complex: where U 0 ⊂ H 1 , U 1 ⊂ H(div) and U 2 ⊂ L 2 , and where U 0 ⊂ H 1 , U 1 ⊂ H(curl) and U 2 ⊂ L 2 . As mentioned previously, this arises due to two possible identifications of differential 1-forms with vector fields.
Note that the vector-valued space is a direct sum of two 'product' spaces that have had the HDiv or HCurl operator applied to them. It can be easily verified that the intersection of the two spaces is trivial, and so the usual sum is a direct sum. In addition, the nodes fit together naturally: for example, HDiv(U 0 ⊗V 1 ) has nodes on (0, 1) facets, while HDiv(U 1 ⊗V 0 ) has nodes on (1, 0) facets. Any interior nodes of HDiv(U 0 ⊗V 1 ) are evaluation of the first component of the vectorvalued function at certain points, while any interior nodes of HDiv(U 1 ⊗ V 0 ) are evaluation of the second component. A similar statement can be made for the H(curl) element.
Again, the vector-valued spaces are direct sums of 'product' spaces that have had the HDiv or HCurl operator applied to them. Again, the nodes fit together naturally: for example, in the H(curl) element, HCurl(U 0 ⊗V 1 ) has nodes on (0, 1) edges, while HCurl(U 1 ⊗V 0 ) has nodes on (1, 0) edges. HCurl(U 0 ⊗V 1 ) may have nodes on (1, 1) facets; if so, these are evaluation of one of the two tangential components. HCurl(U 1 ⊗V 0 ) may have nodes on (1, 1) facets; if so, these are evaluation of the other tangential component. HCurl(U 1 ⊗ V 0 ) may also have nodes on (2, 0) facets; if so, these are evaluation of both tangential components. HCurl(U 0 ⊗ V 1 ) may have interior nodes; if so, these are evaluation of one of the three vector components. HCurl(U 1 ⊗ V 0 ) may have interior nodes; if so, these are evaluation of the other two vector components. A similar, though less complicated, situation exists for the H(div) element.

Numerical quadrature on product cells
Many operations in the finite element method involve integrals of piecewise-polynomial functions over some part of the global mesh. This is typically implemented by performing quadrature in the reference cell.
The basis functions constructed in (2.25) are products of polynomials in distinct sets of variables. This suggests only using a quadrature rule capable of integrating, for example, terms of the form rather than a more expensive quadrature rule capable of integrating (2.44) Suppose we are given existing quadrature rules capable of integrating a polynomial f of degree p over K A ⊂ R n , and of integrating a polynomial g of degree q over K B ⊂ R m , say ( 2.45) where the w i and ξ i are quadrature weights and points, respectively. Then a quadrature rule capable of integrating a polynomial h defined on K A × K B ⊂ R n+m with maximal degree p in x 1 , . . . , x n and q in x n+1 , . . . , x n+m , as in (2.43), is

Product complexes in finite element exterior calculus
This section presents the results of subsection 2.5 and subsection 2.6 in the language of differential forms, which can be considered a generalisation of scalar and vector fields. The reader unfamiliar with finite element differential forms may find the relevant material in, for example, Arnold et al. (2006Arnold et al. ( , 2010. The following definitions can be found in , which contains a far more thorough overview than is given here. In R n , both 0-forms and n-forms can be identified with scalar fields. 1-forms and (n-1)-forms can be identified with vector fields. In three dimensions, 0-forms and 3-forms are identified with scalar fields, while 1-forms and 2-forms are identified with vector fields. In two dimensions, 0forms and 2-forms are identified with scalar fields. 1-forms are identified with vector fields, but this can be done in two different ways since 1-forms and (n-1)-forms coincide in two dimensions. This results in two possible vector fields, which differ by a 90-degree rotation; this was alluded to in the discussion of finite element spaces in two dimensions in subsection 2.4. In one dimension, both 0-forms and 1-forms are conventionally identified with scalar fields.
where each U k is a space of differential k-forms on K A , and a de Rham subcomplex on K B , where each V k is a space of differential k-forms on K B , the tensor product of the complexes (2.49) and (2.50) is a de Rham subcomplex on K A × K B : for k = 0, 1, . . . , n + m.
Note that (U ⊗V ) k is a space of k-forms on K A ⊗ K B , and can hence be interpreted as a scalar or vector field if we are in 2 or 3 spatial dimensions. It can be easily verified that the definitions in (2.51) and (2.52) give rise to (2.30) and (2.34) in two dimensions, and (2.38) in three. The discussion in subsubsection 2.5.2 and subsubsection 2.5.4 on the product of function spaces can be summarised by the definition of ⊗ on the right-hand side of (2.52), along with the definition of the standard wedge product of different forms.
It is clear that much of the apparent complexity of the HDiv and HCurl operators introduced in subsection 2.5 arises from trying to work with only scalars and vectors without introducing differential forms.

Implementation
The product elements, constructed in the previous section, have been implemented within Firedrake . Firedrake is an "automated system for the portable solution of partial differential equations using the finite element method".
Firedrake has several dependencies. Some of these are components of the FEniCS Project (Logg et al., 2012a): FIAT FInite element Automatic Tabulator (Kirby, 2004(Kirby, , 2012, for the construction and tabulation of finite element basis functions UFL Unified Form Language (Alnaes et al., 2014;Alnaes, 2012), a domain-specific language for the specification of finite element variational forms FFC FEniCS Form Compiler (Kirby and Logg, 2006;Logg et al., 2012b), for the generation of low-level C kernels from UFL forms; Firedrake uses a modified FFC that generates abstract syntax trees Dependencies unrelated to the FEniCS Project include: COFFEE (Luporini et al., 2014), a domain-specific compiler which optimises the abstract syntax trees generated by FFC to improve instruction-level parallelism PyOP2 (Rathgeber et al., 2012), a framework for carrying out parallel computations on unstructured meshes PETSc (Balay et al., 2014), which provides mesh objects and an extensive set of linear and nonlinear solvers The main Firedrake module provides a public interface and links the above dependencies together as required.
The changes required to effect the generation of product elements were largely confined to FIAT, UFL and FFC. We therefore begin this section with more detailed expositions on these components of Firedrake. We then discuss the implementation of product cells in subsection 3.4 and the implementation of product finite elements in subsection 3.5. We talk about the resulting algebraic structure in subsection 3.6. We discuss exact quadrature on product cells in subsection 3.7 and the new regions of integration in subsection 3.8. The necessary extensions to facilitate strong boundary conditions are discussed in subsection 3.9. We finish this section by listing some of the finite element families we are capable of generating, in subsection 3.10.

FIAT
FIAT, the FInite element Automatic Tabulator (Kirby, 2004(Kirby, , 2012, is responsible for computing finite element basis functions for a wide range of finite element families, including all those mentioned in subsection 2.4. To do this, it works with an abstraction based on Ciarlet's definition of a finite element, as given in subsection 2.1.
The reference cell K is defined as a set of vertices with specified coordinates together with higher-dimensional geometrical objects which are defined as sets of vertices. The polynomial space P is defined implicitly through a prime basis: typically, an orthonormal set of polynomials, possibly with additional polynomials for elements in which P is not just some P k . The set of nodes N is also defined; this implies the existence of a nodal basis for P, as explained previously.
The nodal basis, which is important in calculations, can be expressed as linear combinations of prime basis functions. This is done by evaluating the known prime basis at the nodes, and inverting the resultant generalised Vandermonde matrix. The main method of interacting with FIAT is through requesting the values of the nodal basis functions tabulated at a set of points inside K. FIAT also stores the geometric decomposition of nodes relative to the topological entities of K.
We have treated the tabulation of basis functions for existing finite elements as a black-box which we did not modify. However, the new product elements have been implemented within FIAT; tabulation of basis functions makes use of the tabulation functions for existing elements. The geometric decomposition of nodes is also calculated from the geometric decompositions of existing elements, as discussed in subsubsection 2.5.3.

UFL
UFL, the Unified Form Language (Alnaes et al., 2014), is a domain-specific language, embedded in Python, for representing weak formulations of partial differential equations. It is centred around expressing multilinear forms: maps from the product of some set of function spaces {V j } ρ j=1 into the real numbers which are linear in each argument. ρ is typically 0, 1 or 2; these are labelled functionals, linear forms and bilinear forms respectively. For ρ ≥ 1, V 1 is referred to as the test space. For ρ = 2, V 2 is referred to as the trial space. The form may additionally be parameterised over one or more coefficient functions, in which case the form is a mapping from the product of coefficient spaces {W k } n k=1 and argument spaces: the form is not required to be linear in the coefficient functions.
We can assume that the function spaces are finite element spaces; in UFL, these are represented by the FiniteElement class. This takes in three arguments: a string representing the family, the geometric cell, and the polynomial degree. A limited amount of symbolic manipulation on FiniteElement objects could previously be done: the UFL EnrichedElement class is used to represent the ⊕ operator discussed in subsection 2.2, and the MixedElement class is used when solving for more than one function simultaneously.
It is assumed that the form can be decomposed into integrals over cells, exterior facets and interior facets. In UFL, these integration regions, or measures, are denoted by dx, ds and dS, respectively. UFL also supports differential operators, symbolic differentiation, discontinuous Galerkin operators and conditionals. UFL has particularly rich support for tensor algebra operations and index notation.
The main extensions that were made to UFL were the addition of operators to allow the creation of products of FiniteElement objects. As discussed in subsubsection 2.5.6, we require an operator for taking the product of elements, called OuterProductElement, and operators for generating an H(curl) or H(div) element, called HCurl and HDiv. Another significant change was the addition of new integration measures: on simplices, all facets are completely equivalent, but on product cells, this is not the case. The new integration measures allow selective integration over the different types of facets that are present.

FFC
FFC, the FEniCS form compiler (Kirby and Logg, 2006;Logg et al., 2012b), was mainly responsible for converting the UFL representation of a form into C++ code for the computation of the local element tensor. It was also responsible for the generation of the local-to-global mappings. Within Firedrake, the role of FFC is somewhat reduced: it no longer handles the local-to-global mappings, and it converts the UFL representation of a form into a lower-level abstract syntax tree (rather than a C++ string), which can be further optimised by COFFEE before compilation.
Computation of the local element tensor, in Firedrake's FFC, is done via numerical quadrature. The form is evaluated at the quadrature points in the reference cell; the necessary quadrature points and weights are given by a combination of FFC and FIAT. To evaluate the form at quadrature points, it is necessary to have the values of the basis functions evaluated at the quadrature points; this information is obtained from FIAT. FFC can also use appropriate Piola transforms to map values between reference and physical space (Rognes et al., 2009).
The main extensions to FFC included adding functionality for the generation of quadrature points and weights for product cells, particularly when integrating forms over integration regions specific to product-cells. Previously, UFL EnrichedElement objects were handled by FFC without the 'sum element' being built in FIAT. To have full support for 'products of sums', this functionality was recreated in FIAT and removed from FFC.

Implementation of product cells
To implement a product element, the geometric product reference cell must be constructed. The mathematical details were described in subsubsection 2.5.1. Separate changes were made in UFL and FIAT.
In UFL, we defined an OuterProductCell class, which inherits from the base class Cell. The OuterProductCell takes two existing cells as arguments into its constructor. Two essential properties of UFL cells are the topological dimension and the geometric dimension, which we will refer to as tdim and gdim, respectively. The topological dimension is the number of spatial dimensions of the reference cell, while the geometric dimension is the number of dimensions of the space in which the physical mesh lies. These are typically the same, but may differ -for example, if a PDE is solved on the surface of a sphere, the topological dimension is 2 but the geometric dimension is 3. The tdim = gdim cases were implemented in Rognes et al. (2013).
The tdim of an OuterProductCell is trivially the sum of the tdims of the two constituent cells. The definition of gdim is not so obvious, since it relates to properties of the global mesh, while our product of elements is local to a cell. We therefore defined the default gdim to be the maximum of the tdim of the product cell and the gdims of the constituent cells. We allow this to be overridden during construction of the product cell: our definition only amounts to a heuristic, and a higher-level application is more likely to have the necessary information about the global mesh at its disposal. Two UFL OuterProductCells compare A product cell has two distinct types of facets. A facet can be expressed as the geometric product of one constituent cell's facet with the other constituent cell -for example, the facets of a triangular prism (the product of a triangle with an interval) are squares (the product of a triangle's facet with an interval) or triangles (the product of a triangle with an interval's facet). This can naturally be simplified if one of the product constituents is a vertex! The OuterProductCell class stores both types of facet, as either OuterProductCells or ordinary Cells.
A representation of a product cell was also created within FIAT. The coordinates of a vertex of the product cell are given by concatenating the coordinates of a vertex of the first cell with the coordinates of the vertex of the second cell. This must be done for every possible pairing in order to obtain the full set of vertices for the product cell. Let the constituent cells be K 1 and K 2 with n 1 and n 2 vertices, respectively. Assuming existing orderings for the vertices of K 1 and K 2 , with the numbering beginning from 0, there is a natural ordering for the vertices of the product cell: vertex i of K 1 and vertex j of K 2 are concatenated to produce vertex i · n 2 + j of K 1 × K 2 . This naturally runs from 0 to (n 1 n 2 − 1).
FIAT also stores information on the topological structure of cells; this is essential for the later association of nodes with topological entities. For existing simplicies, the topology was represented as a (Python) dictionary whose keys are the dimension of the topological entity. Each possible dimension then contains a sub-dictionary, in which the keys are simply 0 up to the number of entities of that dimension, and each value is a list of vertices corresponding to that entity. For example, the representations for a triangle and an interval are given in Listing 1 and Listing 2, respectively.
The topology of the product cell is generated from the topology of the constituent cells: every entity of the product cell is the product of entities of the constituent cell. As mentioned in subsubsection 2.5.1, it is no longer sensible to index by dimension alone, since entities of the same dimension are no longer fully interchangeable. Instead, we index by a pair of indices; sample Python code is given in Listing 3. If we take the product of the triangle and interval given above, the topology of the product cell is shown in Listing 4.  The letters v and c denote vertices and cells, respectively. These correspond to entities of dimension 0 and 1.
Listing 3: Algorithm producing the representation of a product reference cell in FIAT

Implementation of product finite elements
To implement product finite elements, additions to UFL and FIAT were required. The relevant numerical calculations are performed in FIAT; the UFL changes are purely symbolic.
As discussed in subsubsection 2.5.6, we implement several new classes in UFL: the base OuterProductElement operator, and the 'modifiers' HCurl and HDiv. The existing UFL FiniteElement classes are remarkably lightweight; they store little more than the degree and the value shape of the finite element. The degree is the maximal degree of any polynomial basis function -storing this information allows determination of an appropriate quadrature degree when integrating an entire form. By default, exact quadrature is used where possible, though this can be manually overridden. The value shape represents whether the element is scalar-valued or vector-valued and, if applicable, the dimension of the vector in physical space.
For OuterProductElements, we define the degree to be a Python tuple of length two: as discussed in subsection 2.7, the basis functions are products of polynomials in distinct sets of variables and so it is reasonable to store the degrees separately. The degree can then be represented by (A degree, B degree) where A degree and B degree are the degrees of the constituent elements. The value shape definition is simple: if both elements are scalarvalued then so is their product. If one of the elements is vector-valued then, since the other element must be scalar-valued, their product is vector-valued with the same dimension -recall that OuterProductElement does not do padding.
For HCurl and HDiv elements, the degree is identical to the degree of the underlying OuterProductElement. The value shape, on the other hand, requires modification. In Listing 5: Algorithm to produce the geometric decomposition of product nodes in FIAT physical space, these vector-valued elements have dimension equal to the geometric dimension gdim -of the product cell.
The FIAT implementation of elements serves two main purposes: it stores a representation of the geometric decomposition of nodes between topological entities of the cell, and it allows the corresponding basis functions to be evaluated at arbitrary points in the reference cell. The former is important for enforcement of inter-cell continuity via local and global numbering, while the latter facilitates numerical quadrature by providing numerical arrays to be used in low-level kernels.
For product elements, the geometric decomposition of nodes was described in subsubsection 2.5.3. An algorithm for producing this from the geometric decompositions of the constituent elements is given in Listing 5; note the similarity with the algorithm for producing the product cell. The actual nodes exist as objects within FIAT, although they have almost no purpose in calculation. Nevertheless, this was implemented to conform with the rules derived in subsubsection 2.5.5.
The main role of FIAT is to tabulate finite element basis functions, and derivatives thereof, at specified points in the reference cell. The tabulate method of a FIAT finite element accepts two arguments: the maximal order of derivatives to tabulate, and a list of points. The first step is to generate corresponding lists of points at which to tabulate the constituent elements. By tabulating the constituent finite elements at these points up to the same order of derivatives, it is possible to manipulate the resulting arrays to generate the tabulation of the product element. We reduced most of the work to the manipulation of NumPy arrays (Walt et al., 2011) to replicate the discussion in subsubsection 2.5.4. We remark that, when derivatives of basis functions are required, there is no need to explicitly use the product rule for differentiation since the derivative acts on only one of the constituent elements.

Algebraic structure
The extensions described in subsection 3.5 enable sophisticated manipulation of finite elements within UFL. For example, recall the complex on triangles as given in (2.13): Suppose we wish to take the product of this with some complex on intervals, such as This generates a complex on triangular prisms: where W 0 := (P 2 ⊕ B 3 ) ⊗ P 1 , (3.5) W 1 := HCurl((P 2 ⊕ B 3 ) ⊗ P 1 ) ⊕ HCurl(BDFM 2 ⊗ P 1 ), (3.6) W 2 := HDiv(BDFM 2 ⊗ DP 0 ) ⊕ HDiv(DP 1 ⊗ P 1 ), (3.7) W 3 := DP 1 ⊗ DP 0 . (3.8) We see that W 1 is a sum of HCurls of product elements, one of which already includes a sum element! Following our extensions to UFL, the product complex may be constructed as shown in Listing 6.
There is also partial support for recursive products, as in P 2 ⊗ (DP 1 ⊗ P 2 ), (3.9) which is necessary to produce elements on hexahedra.

Implementation of quadrature
As discussed in subsection 2.7, Firedrake uses numerical quadrature to evaluate integrals. The number and location of quadrature points is sufficient to integrate a polynomial of a particular degree; where possible, this is the actual polynomial degree of the form. If this is not possible (for example, if the form involves a division, or includes a trigonometric function), a heuristic degree, based on the operators in the form, is used for quadrature purposes. On simplices, a single polynomial degree is used to describe the form, since the basis functions have certain symmetry properties and there is no preferred direction. On product cells, while a single polynomial degree could be used, this is inefficient for the reasons given previously.
UFL automatically calculates the polynomial degree of a given form by traversing the abstract syntax tree representing the form. The rules applied are straightforward: for example, if the abstract syntax tree contains a multiplication operator, the degrees of the two subexpressions are added together. To make this compatible with product cells, the main change is to modify the traversal of the abstract syntax tree to operate on polynomial degrees represented by tuples as well as integers.
There is one small subtlety: on a simplex cell, a differential operator would unconditionally decrement the polynomial degree. On a product cell, it is clearly incorrect to decrement both entries of the degree tuple. Only one entry of the tuple needs to be decremented, depending on the direction in which the differential operator acts; the situation is somewhat complicated by the fact that the differential operators act in physical space rather than reference space. Ideally, there would be support for, for example, vector quantities whose components have separate degrees. Since this does not exist in UFL, we decided to simply prevent differential operators from decrementing the degree of the form when product cells are used. This does not affect correctness; it merely means that more quadrature points may be used to evaluate an integral than are strictly necessary. Since there are rarely more than two derivatives multiplied together in a single form, we believe that the negative impact is not serious. It is always possible to override the quadrature degree manually; this is often used for complicated forms anyway.

Support for new integration regions
On simplicial cells, Firedrake supports three types of integrals: integrals over cells, integrals over exterior facets and integrals over interior facets. Integrals over exterior facets are typically used to apply boundary conditions weakly, while integrals over interior facets are used to couple neighbouring cells when discontinuous function spaces are present. The implementation of the different types of integral is quite elegant: the only difference between integrating a function over the interior of the cell and over a single facet is the choice of quadrature points and quadrature weights.
The cell integrals are implemented using (2.46), with the generation of quadrature points and quadrature weights occuring within FIAT. These are generated from the quadrature points and weights of the constituent elements, following (2.47) and (2.48).
As we have mentioned several times, on simplices, all facets are equivalent. On product cells, this is not the case -for example, on triangular prisms, some facets are triangles while other facets are quadrilaterals. However, all facets (and indeed all entities) of a product cell can be considered as a product of entities on the constituent cells. We can therefore apply (2.46) to integrate over facets of the product cell, making use of the existing quadrature rules for constituent cells and facets thereof.

Support for boundary conditions
Boundary conditions may be implemented in a weak sense by using appropriate facet integrals. However, strong boundary conditions are implemented by removing the appropriate degrees of freedom from the global linear system. This requires the identification of degrees of freedom on (a subset of) the boundary of the global mesh. Due to details specific to Firedrake's implementation of strong boundary conditions, this requires the identification of nodes on facets of the local cell. In an early implementation, this was done "on-the-fly" by taking the correct products of nodes of the constituent cells. However, this was later updated to a single call to the entity closure dofs method, which lists nodes associated with the closure of each topological entity.
Firedrake therefore supports strong boundary conditions on the 'top', 'bottom', or 'sides' of an extruded mesh. Naturally, Firedrake also supports integration over these same regions, for the application of weak boundary conditions.

Example product finite element families
We give explicit constructions of several element families on product cells. Many of these have well-established names in the literature, which we attempt to use.
The Brezzi-Douglas-Marini face element on triangles contains all polynomials up to a given degree, and has continuous normal component across edges. An analogous, though perhaps naïve, product element may be built on quadrilaterals. However, the name BDM is already reserved for a different family of H(div) elements on quadrilaterals which do not have a product structure (Brezzi et al., 1985). The desired element can alternatively be thought of as a 2D reduction of the 3D 'Nédélec second-kind' elements introduced in Nédélec (1986), and so we will use the name Ndtwo. The construction is given in Listing 11.
The "Nédélec elements of the second kind" (Nédélec, 1986) contain all polynomials up to a given degree in each variable separately. Constructions are given in Listing 16 and Listing 17.
As for hexahedra, there are also Nédélec elements of the second kind (Nédélec, 1986); constructions are given in Listing 22 and Listing 23.
Tests are performed in both two and three spatial dimensions. We make use of Firedrake's ExtrudedMesh functionality. In two dimensions, the cells are quadrilaterals, usually squares. In three dimensions, the cells are triangular prisms.

Vector Laplacian
We seek a solution to in a domain Ω, with boundary conditions on ∂ Ω, where n is the outward normal. A naïve discretisation can lead to spurious solutions, especially on non-convex domains, but an accurate discretisation can be obtained by introducing an auxiliary variable (see, for example, Arnold et al. (2010)): be finite element spaces. A suitable formulation of (4.1), written in weak form, is then: find σ ∈ V 0 , u ∈ V 1 such that for all τ ∈ V 0 , v ∈ V 1 . For brevity, we have used angled brackets to denote the standard L 2 inner product: The boundary conditions have been implicitly applied, in a weak sense, through neglecting the surface terms when integrating by parts.

2D
We take Ω to be the unit square [0, 1] 2 . The choice f = π 2 (k 2 + l 2 ) sin(kπx) cos(lπy) cos(kπx) sin(lπy) (4.9) produces the solution u = sin(kπx) cos(lπy) cos(kπx) sin(lπy) ,  Figure 4: The L 2 error between the computed and 'analytic' solution is plotted against ∆x for the 2D problem described in subsubsection 4.1.1. The dotted lines are proportional to ∆x n , for n from 1 to 4, and are merely to aid comprehension. The convergence rates are as expected: the Q 1 -RTE 1 case demonstrates second-order convergence for u and first-order convergence for σ . Recall that, in our numbering convention, RTE 1 does not contain all linear functions, and so we can only expect first-order convergence. The superconvergence of u is unsurprising, due to the regular nature of the mesh. The higher-order cases Q 2 -RTE 2 and Q 3 -RTE 3 converge at correspondingly faster rates.
To discretise this problem, we subdivide Ω into squares with side length ∆x. We use Q n for the H 1 space, and RTE n for the H(curl) space, as defined in subsection 3.10, for n from 1 to 3. We take k and l to be 1 and 2, respectively. We approximate f by 'interpolating' the analytic expression onto a vector-valued function in Q n+1 . The L 2 errors between the calculated and 'analytic' solutions for varying ∆x are plotted in Figure 4. This is done for both u and σ ; the socalled analytic solutions are also approximations which are formed by interpolating the genuine analytic solution onto nodes of Q n+1 .

3D
We take Ω to be the unit cube [0, 1] 3 . The choice f = π 2   (k 2 + l 2 ) sin(kπx) cos(lπy) (l 2 + m 2 ) sin(lπy) cos(mπz) (k 2 + m 2 ) sin(mπz) cos(kπx) To discretise this problem, we subdivide Ω into triangular prisms whose base is a right-angled triangle with short sides of length ∆x and whose height is ∆x. We use the Q n prism element for the H 1 space, and the N1E n prism element for the H(curl) space, for n from 1 to 3. These are defined in subsection 3.10. We take k, l and m to be 1, 2 and 3, respectively. As before, we approximate f by 'interpolating' the analytic expression onto a vector-valued function in Q n+1 . The L 2 errors between the calculated and 'analytic' solutions for varying ∆x are plotted in Figure 5. This is done for both u and σ ; the so-called analytic solutions are also approximations which are formed by interpolating the genuine analytic solution onto nodes of Q n+1 .

Gravity wave
A simple model of the atmosphere is given by 14) ∂ p ∂t = −c 2 ∇ · u, (4.15) along with the boundary condition u · n = 0, (4.16) where n is a unit normal vector. The prognostic variables are the velocity, u, the pressure perturbation, p, and the buoyancy perturbation, b. The constants N and c represent the Brunt-Väisälä frequency and the speed of sound respectively, whileẑ represents a unit vector opposite to the direction of gravity. These are a reduction of, for example, the linear equations (17)- (21) in Skamarock and Klemp (1994), in which we have neglected the constant background velocity and the Coriolis term, and rescaled θ by θ 0 /g to produce b.
Given some product complex  Figure 5: The L 2 error between the computed and 'analytic' solution is plotted against ∆x for the 3D problem described in subsubsection 4.1.2. The dotted lines are proportional to ∆x n , for n from 1 to 4, and are merely to aid comprehension. The convergence rates are identical to the rates in the earlier 2D problem.
as in (2.38), we can seek a solution with u ∈ W 0 2 , b ∈ W v 2 and p ∈ W 3 . W 0 2 is the subspace of W 2 whose normal component vanishes on the boundary of the domain -this enforces the boundary condition (4.16). We have used W v 2 to represent the "vertical" part of W 2 : if we write W 2 as a sum of two product elements HDiv(U 1 ⊗ V 1 ) and HDiv(U 2 ⊗ V 0 ) then W v 2 is the scalar-valued product U 2 ⊗ V 0 . This combination of finite element spaces for u and b is analogous to the Charney-Phillips staggering of variables in the vertical direction (Charney and Phillips, 1953).
We can then write (4.13)-(4.15) in semi-discrete form as follows: for all w ∈ W 0 2 , γ ∈ W v 2 , φ ∈ W 3 . We have integrated by parts to produce (4.18); this avoids taking a derivative of the discontinuous pressure field.
It can be easily verified that (4.13)-(4.15) along with the boundary condition (4.16) lead to conservation of the quantity This quantity can be interpreted as a total perturbation energy, and the three terms in (4.21) can be interpreted as kinetic energy (KE), potential energy (PE) and internal energy (IE) respectively. The semi-discretisation given in (4.18)-(4.20) conserves total perturbation energy. If we discretise in time using the implicit midpoint rule, which preserves quadratic invariants (see, for example, Leimkuhler and Reich (2005)), then the fully discrete system will conserve energy as well. We verify this in two and three dimensions.

2D
We take the domain to be a rectangle of width 6000km and height 10km, with periodic boundaries in the horizontal direction. This is divided into rectangular cells of width 20km and height 1km. We use a Brunt-Väisälä frequency of 10 −2 s −1 and a sound speed of 300ms −1 . We start the simulation from rest, with a buoyancy perturbation together with a vertically balancing pressure field where H is the height of the domain, a is a horizontal length-scale, which we take to be 100km, and x 0 is the horizontal midpoint of the domain. We use a timestep of 480s, and run for a total of 120,000s.
To discretise this problem, we take W 2 to be RTF 2 , and W 3 to be DQ 1 . The initial conditions are interpolated into the buoyancy and pressure fields. The energy is calculated at every time step; the results are plotted in Figure 6. The total energy is conserved to roughly one part in 5 × 10 6 , which is comparable to the linear solver tolerances.

3D
We take the domain to be a spherical annulus, centred at the origin, whose inner radius is approximately 6371km, with a thickness of 10km. This is divided into triangular prism cells of side-length approximately 1000km and height 1km. We again use a Brunt-Väisälä frequency of 10 −2 s −1 and a sound speed of 300ms −1 . We start the simulation from rest, with a buoyancy where H is the height of the domain, a is the inner radius, and L is a horizontal length-scale, which we take to be 500km. We use a timestep of 1920s, and run for a total of 480,000s.
To discretise this problem, we use the product complex formed from the BDFM 2 complex on triangles (2.13) and the P 2 -DP 1 complex on intervals. The initial conditions are interpolated into the buoyancy and pressure fields. The energy is calculated at every time step; the results are plotted in Figure 7. The total energy is conserved to roughly one part in 1.4 × 10 8 , which is comparable to the linear solver tolerances.

DG advection
The advection of a scalar field q by a known, divergence-free velocity field u 0 can be described by the equation ∂ q ∂t + ∇ · ( u 0 q) = 0. (4.26) Let φ be an arbitrary smooth test-function. Multiplying (4.26) by φ and integrating over the domain Ω gives Now, let us decompose Ω into a number of cells. Suppose we take q to be in a discontinuous function space, V , and use test-functions from V . The expression ∇ · ( u 0 q), in (4.27), is then troublesome because it involves the derivative of a discontinuous quantity. To avoid this, we instead integrate (4.26) over a single cell, e: (4.28) We now apply integration by parts to move the derivative onto the test-functions φ , which are only non-zero on a single cell: (4.29) where ∂ e is the cell boundary, n is an outward-pointing normal, and dS is an appropriate integration measure. Note that q is not well-defined on cell boundaries since it is in a discontinuous function space; we choose to use the upwind value, which we denote by q. Summing (4.29) over all cells in the mesh then gives where the integrals over cell boundaries in (4.29) have been amalgamated into an integral over all exterior mesh facets and an integral over all interior mesh facets. Note that the integral on each interior facet has two contributions, corresponding to the two neighbouring cells; the sign of the test function depends on whether the corresponding outward cell normal is aligned or anti-aligned with the canonical facet normal.

2D
We take Ω to be the unit square [0, 1] 2 . Our initial condition will be a cosine hill with radius r 0 = 0.15, centred at x 0 = (0.25, 0.5). The prescribed velocity field is u 0 ( x,t) = cos πt T sin(πx) 2 sin(2πy) − sin(πy) 2 sin(2πx) , (4.36) as in LeVeque (1996). This gives a reversing, swirling flow field which vanishes on the boundaries of Ω. The initial condition should be recovered at t = T . Following LeVeque (1996), we take T = 3 2 . To discretise this problem, we subdivide Ω into squares with side length ∆x. We use DQ n for the discontinuous function space, for n from 0 to 2. We initialise q by interpolating the expression given in (4.35) into the appropriate space. We approximate u 0 by interpolating the expression given in (4.36) onto a vector-valued function in Q 2 . The L 2 errors between the initial and final q fields for varying ∆x are plotted in Figure 8.

3D
We take Ω to be the unit cube [0, 1] 3 . We use the same, z-independent, cosine hill for the initial condition: Again, this gives a reversing, swirling flow field which vanishes on the boundaries of Ω. The initial condition should be recovered at t = T . We take T = 1 4 . To discretise this problem, we subdivide Ω into right-angled triangular prisms whose base has short side ∆x, and whose height is 16∆x. We use the discontinuous DQ n space for q, where n ranges from 0 to 2. We initialise q by interpolating the expression given in (4.37) into the appropriate space. We approximate u 0 by interpolating the expression given in (4.38) onto a vector-valued function in Q 2 . The L 2 errors between the initial and final q fields for varying ∆x are plotted in Figure 9.  Figure 8: The L 2 error between the computed and 'analytic' solution is plotted against ∆x for the 2D problem described in subsubsection 4.3.1. The dotted lines are proportional to ∆x and ∆x 2 , and are merely to aid comprehension. The DQ 0 simulations converge at first-order for sufficiently small values of ∆x. The DQ 1 simulations converge at second-order, as expected. The cosine bell initial condition has a discontinuous second derivative, which inhibits the DG 2 simulations from exceeding a second-order rate of convergence.
10 -3 10 -2 10 -1 ∆x 10 -5 10 -4 10 -3 10 -2 10 -1 L 2 norm of error DQ 0 DQ 1 DQ 2 Figure 9: The L 2 error between the computed and 'analytic' solution is plotted against ∆x for the 3D problem described in subsubsection 4.3.2. The dotted lines are proportional to ∆x and ∆x 2 , and are merely to aid comprehension. The DQ 0 simulations converge at firstorder for sufficiently small ∆x. The DQ 1 simulations converge at second-order, as expected. The cosine bell initial condition has a discontinuous second derivative, which inhibits the DG 2 simulations from exceeding a second-order rate of convergence.

Limitations and extensions
There are several limitations of the current implementation, which leaves scope for future work. The most obvious is that the quadrature calculations are relatively inefficient, particularly at high order. The product structure of the basis functions can be exploited to generate a more efficient implementation of numerical quadrature. This can be done using the sum-factorisation method, which lifts invariant terms out of the innermost loop. In the very simplest cases, direct factorisation of the integral may be possible. It would have been possible to implement these operations within the FEniCS Form Compiler, FFC. However, this would mask the underlying issue -that interaction with FIAT, which is wholly responsible for producing the relevant finite elements, is essentially limited to passing in a set of points, and receiving an array of basis functions evaluated at this set of points. Any underlying basis function structure is lost in this process. Work is underway on a more sophisticated layer of software that returns an algorithm for performing a given operation on a finite element, rather than merely an array of numbers.
Another limitation is that, in the current version of FFC, the Jacobian of the coordinate mapping from reference to physical space is assumed to be constant in each cell. This is satisfactory for simplices, since the physical and reference cells can always be linked by an affine transformation. However, this is not the case for general quadrilateral, triangular prism and hexahedral meshes. The Jacobian used with product cells is currently calculated using a (strict) subset of the cell's vertices. It could be argued that it would be more desirable to use the Jacobian at the midpoint of the cell, instead. Indeed, effecting this would be a trivial change. However, it would clearly not resolve the underlying 'issue' -that the Jacobian will, in general, vary across the cell. For similar reasons, the normal vector to a quadrilateral side facet of a triangular prism is calculated using a subset of the facet vertices, and therefore will be inexact if the facet is not planar.
A planned feature is the introduction of curved, or 'bendy', cells, in which the coordinate transformation is quadratic or higher-order. This allows, for example, more faithful representations of a sphere or a spherical annulus. Clearly, the coordinate transformation to produce a curved cell also has a spatially-varying Jacobian. It is therefore envisaged that the previous limitation will be removed when support for bendy cells is added. Another planned Firedrake feature is full support for general unstructured quadrilateral meshes; the subsequent extrusion of an quadrilateral base mesh produces a mesh of hexahedral cells.

Conclusion
This paper presented the extensions to the automated code generation pipeline of Firedrake to facilitate the use of finite element spaces on non-simplex cells, in two and three dimensions. All numerical experiments given in this paper were performed with the latest version of Firedrake, currently at commit 091d3a5. The corresponding commits for other components are PyOP2 3d1416f, FFC 52f5e65, UFL 8c11f0e and FIAT bf85181. The examples made extensive use of the recently-added extruded mesh functionality in Firedrake; a related paper detailing the implementation of extruded meshes is in preparation.