Automated generation and symbolic manipulation of tensor product finite elements Andrew T. T. McRae1,2,* , Gheorge-Teodor Bercea3 , Lawrence Mitchell2,3 , David A. Ham2,3 , and Colin J. Cotter2

arXiv:1411.2940v1 [math.NA] 11 Nov 2014

1 The

Grantham Institute, Imperial College London, London, SW7 2AZ, UK of Mathematics, Imperial College London, London, SW7 2AZ, UK 3 Department of Computing, Imperial College London, London, SW7 2AZ, UK * Correspondence to: [email protected]

2 Department

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. Keywords: tensor product finite element; finite element exterior calculus; automated code generation

1 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 (Alnæs et al.,

1

2014; Alnæs, 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, 2012). These local assembly kernels are executed on each cell1 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 (Rathgeber et al., 2014; Rathgeber, 2014). 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 1 Note on terminology:

throughout this paper, we use the term ‘cell’ to denote the geometric component of the mesh; we reserve the term ‘finite element’ to denote the space of functions on a cell and supplementary information related to global continuity.

2

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, 2010, 2014). 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.

2 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. 2.1 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 Rn , 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 = {n1 , . . . , ndim P } is a basis for the dual space P0 – 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,

3

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: ni (Φ j ) = δi j ;

(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. 2.2 Sum of finite elements

Suppose we have two finite elements U = (K, PA , NA ) and V = (K, PB , NB ) defined over the same reference cell K. If the intersection of PA and PB is trivial, we can define the direct sum U ⊕V to be the finite element (K, P, N), where P := PA ⊕ PB ≡ { fA + fB | fA ∈ PA , fB ∈ PB }

(2.2)

N = NA ∪ NB .

(2.3)

and Even though the intersection of PA and PB is trivial, it is not generally true that concatenating the existing nodal bases for PA and PB 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. 2.3 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.

4

We start by defining L2 (Ω) 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 L2 (Ω). 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 Z

f 2 + ∇ f · ∇ f dx < ∞,

(2.4)

Ω

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 facet2 between neighbouring cells. Given a function ~u in some vector-valued finite element space, we can also say ~u ∈ H 1 (Ω) (resp. L2 (Ω)) if each scalar component is in H 1 (Ω) (resp. L2 (Ω)), 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: Z ~u ·~u + (∇ ·~u)2 dx < ∞.

(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: Z

~u ·~u + (∇ ×~u) · (∇ ×~u) dx < ∞.

(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 L2 . 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 L2 , 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 6⊂ H 1 , i.e., we have made the strongest statement possible. In particular, we will use L2 to denote a total absence of continuity between cells. 2 In

three dimensions, the terms ‘vertices’, ‘edges’ and ‘facets’ refer to distinct entities. In fewer dimensions, some of these terms coincide.

5

2.3.2 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

6

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. L2 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 L2 , 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. 2.3.3 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 L2 function from the reference cell to the physical cell, the simplest mapping suffices. Let fˆ(~X) be a function defined over the reference cell. The corresponding function f (~x) defined over the physical cell is then f (~x) = fˆ ◦ φ −1 (~x).

(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: ~f (~x) = 1 J ~fˆ ◦ φ −1 (~x), det J

(2.8)

while H(curl) functions are mapped using the covariant Piola transform, which preserves tangential components: ~f (~x) = J −T ~fˆ ◦ φ −1 (~x). (2.9) Further details on the Piola transform can be found in, for example, Boffi et al. (2013).

7

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). 2.4 Simplicial finite element spaces within finite element exterior calculus

The work of Arnold et al. (2006, 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 d dx Pn −→ DPn−1 ,

(2.10)

where Pn represents the space of functions which are piecewise polynomials of degree n, continuous between cells, and DPn represents the space of functions which are piecewise polynomials of degree n, possibly discontinuous between cells – this is often denoted by PDG n in other literature. The P spaces, being continuous between cells, are subspaces of H 1 , while the DP spaces are only in L2 . 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) are ∇⊥

∇·

Pn −→ RTFn −→ DPn−1 and

∇⊥

∇·

Pn −→ BDMFn−1 −→ DPn−2 ,

(2.11)

(2.12)

where RTFn represents the Raviart–Thomas element (Raviart and Thomas, 1977) and BDMFn 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

8

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 Pn 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 d2 = 0. Another example, highlighted in Cotter and Shipton (2012) as being relevant for numerical weather prediction, is ∇⊥

∇·

P2 ⊕ B3 −→ BDFM2 −→ DP1 ,

(2.13)

where B3 represents a cubic ‘bubble’ function, vanishing at cell edges, and BDFM2 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 RTF1 , rather than RTF0 . Similarly, the BDFM2 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: ∇⊥ ·

∇

Pn −→ RTEn −→ DPn−1 and

∇⊥ ·

∇

Pn −→ BDMEn−1 −→ DPn−2 ;

(2.14)

(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 Pn element, and ∇⊥ · for the vector element, essentially zˆ · ∇×. As before, we have the identity ∇⊥ · ∇ = 0. On tetrahedra (three spatial dimensions), the examples given by Arnold et al. correspond to ∇

∇×

∇·

Pn −→ N1En −→ N1Fn −→ DPn−1 , ∇

∇×

∇·

Pn −→ N1En −→ N2Fn−1 −→ DPn−2 ,

(2.16) (2.17)

∇

∇×

∇·

(2.18)

∇

∇×

∇·

(2.19)

Pn −→ N2En−1 −→ N1Fn−1 −→ DPn−2 , Pn −→ N2En−1 −→ N2Fn−2 −→ DPn−3 ,

where N1E and N2E denote the N´ed´elec edge elements of the 1st and 2nd kind, and N1F and N2F denote the N´ed´elec face elements of the 1st and 2nd kind. These element families were

9

introduced in N´ed´elec (1980, 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 d2 = 0 is encoded in the identities ∇ × ∇ = 0 and ∇ · ∇× = 0. 2.5 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 := (KA , PA , NA ) and V := (KB , PB , NB ) 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. 2.5.1 Product cells

Given reference cells KA ⊂ Rn and KB ⊂ Rm , the reference product cell KA × KB can be defined straightforwardly as follows:  KA × KB := (x1 , . . . , xn+m ) ∈ Rn+m | (x1 , . . . , xn ) ∈ KA , (xn+1 , . . . , xn+m ) ∈ KB . (2.20) The topological entities of KA × KB correspond to products of topological entities of KA and KB . If we label the entities of a reference cell (in Rn , 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 KA × KB can be labelled as follows: (0, 0): vertices of KA × KB – the product of a vertex of KA with a vertex of KB (1, 0): edges of KA × KB – the product of an edge of KA with a vertex of KB (0, 1): edges of KA × KB – the product of a vertex of KA with an edge of KB

.. .

(n-1, m): facets of KA × KB – the product of a facet of KA with the cell of KB (n, m-1): facets of KA × KB – the product of the cell of KA with a facet of KB (n, m): cell of KA × KB – the product of the cell of KA with the cell of KB

It is important to distinguish between different types of entities, even those with the same dimension. For example, if KA is a triangle and KB an interval, the (2, 0) facets of the prism KA × KB are triangles while the (1, 1) facets are quadrilaterals.

10

2.5.2 Product spaces of functions – simple elements

Given spaces of functions PA and PB , the product space PA ⊗ PB can be defined as the span of products of functions in PA and PB : PA ⊗ PB := span { f · g | f ∈ PA , g ∈ PB } ,

(2.21)

where the product function f · g is defined so that ( f · g)(x1 , . . . , xn+m ) = f (x1 , . . . , xn ) · g(xn+1 , . . . , xn+m ).

(2.22)

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 PA ⊗ PB can be constructed from bases for PA and PB . If PA and PB have nodal bases n o n o (A) (A) (A) (B) (B) (B) Φ1 , Φ2 , . . . ΦN , Φ1 , Φ2 , . . . ΦM (2.23) respectively, a nodal basis for PA ⊗ PB is given by  Φi, j , i = 1, . . . , N, j = 1, . . . , M , where

(A)

(B)

Φi, j := Φi · Φ j ,

i = 1, . . . , N, j = 1, . . . , M;

(2.24)

(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 PA ⊗ PB . We will return to this after a brief description of product nodes. 2.5.3 Product nodes – geometric decomposition

Recall that the nodes are a basis for the dual space (PA ⊗ PB )0 , 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 PA0 and PB0 , there is a natural basis for (PA ⊗ PB )0 which (A) is essentially an outer product of the bases for PA0 and PB0 . Let ni, j denote a “product” of ni , (B) the i’th node in NA , with n j , the j’th node in NB – though we will refrain from actually defining (A)

the product of nodes until subsubsection 2.5.5. If ni is associated with an entity of KA of (A) dimension p and n j is associated with an entity of KB of dimension q then ni, j is associated with an entity of KA × KB with label (p, q). This geometric decomposition of nodes in the product element is used to motivate further manipulation of PA ⊗ PB to produce a more natural space of functions, particularly in the case of vector-valued elements.

11

2.5.4 Product spaces of functions – scalar- and vector-valued elements in 2D and 3D 2D

We take the reference cells KA and KB to be intervals, so the product cell KA × KB is twodimensional. The elements on intervals are scalar-valued and are either in H 1 or L2 . There are, therefore, three cases that we will consider: the product of two H 1 elements, the product of two L2 elements, and the product of an H 1 element with an L2 element (or vice-versa). A summary of the below is given in Table 2.1. H 1 × H 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. H 1 × L2 : The H 1 element is guaranteed to have nodes associated with the vertices of its reference cell, while the L2 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 L2 . 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.

12

Product (1D × 1D) H1 × H1 H 1 × L2 H 1 × L2 H 1 × L2 L2 × H 1 L2 × H 1 L2 × H 1 L2 × L2

Components f ×g f ×g f ×g f ×g f ×g f ×g f ×g f ×g

Modifier (none) (none) H(curl) H(div) (none) H(curl) H(div) (none)

Result fg fg (0, f g) ( f g, 0) fg ( f g, 0) (0, f g) fg

Mapping identity identity covariant Piola contravariant Piola identity covariant Piola contravariant Piola identity

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 L2 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. L2 ×L2 : 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 KA ⊂ R2 and KB is an interval, so the product cell KA × KB 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 L2 . The elements on a one-dimensional reference cell, as before, are scalar-valued and are either in H 1 or L2 . 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

13

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 × L2 : 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 L2 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.

14

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 R2 , 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 R3 , 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. L2 × 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.

15

H(div)/H(curl) × L2 : 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 R2 , 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.

L2 × L2 : 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.

2.5.5 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 facet – used 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 KA and KB are intervals; the elements on intervals are scalar-valued and are

16

Product (2D × 1D) H1 × H1 H 1 × L2 H 1 × L2 H(curl) × H 1 H(curl) × H 1 H(div) × H 1 H(div) × H 1 H(curl) × L2 H(curl) × L2 H(div) × L2 H(div) × L2 L2 × H 1 L2 × H 1 L2 × L2

Components f ×g f ×g f ×g ( fx , fy ) × g ( fx , fy ) × g ( fx , fy ) × g ( fx , fy ) × g ( fx , fy ) × g ( fx , fy ) × g ( fx , fy ) × g ( fx , fy ) × g f ×g f ×g f ×g

Modifier (none) (none) H(curl) (none) H(curl) (none) H(curl) (none) H(div) (none) H(div) (none) H(div) (none)

Result fg fg (0, 0, f g) ( fx g, fy g)† ( fx g, fy g, 0) ( fx g, fy g)† (− fy g, fx g, 0) ( fx g, fy g)† ( fy g, − fx g, 0) ( fx g, fy g)† ( fx g, fy g, 0) fg (0, 0, f g) fg

Mapping identity identity covariant Piola * covariant Piola * covariant Piola * contravariant Piola * contravariant Piola identity contravariant Piola identity

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 L2 . All nodes are simply 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 1D element nodes. L2 × L2 : 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 1D element nodes. H 1 × L2 : If the scalar-valued element is used, all nodes remain point evaluation. If the element is padded into an H(div) element, the nodes on facets are the evaluation of the normal component of the vector-valued function. The interior nodes become the evaluation of one component – the ‘horizontal’ component – of the vector-valued function. If the element is padded into an H(curl) element, the nodes on facets are the evaluation of the tangential component of the vector-valued function; the interior nodes are now the evaluation of the ‘vertical’ component of the vectorvalued function. 3D

The reference cell KA ⊂ R2 , while KB is an interval. If the element on KA is in H 1 or L2 , 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

17

the interior nodes are evaluation of the two components of the function at given points. The elements on KB 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. H 1 × L2 : For the scalar-valued product element, all nodes remain point evaluation. If the product element is padded into a vector-valued H(curl) element, the nodes on (0, 1) edges become evaluation of the tangential component. The nodes on (1, 1) facets become evaluation of the ‘vertical’ tangential component, and interior nodes become evaluation of the ‘vertical’ component of the vector-valued function. H(div)/H(curl) × H 1 : Assume the product element, taking values in R2 is padded into an H(curl) element. The nodes on (1, 0) edges become evaluation of the tangential component. The nodes on (2, 0) facets become evaluation of both tangential components. The nodes on (1, 1) facets become evaluation of the ‘horizontal’ tangential component. Interior nodes become evaluation of the two ‘horizontal’ components of the vector-valued function. L2 × H 1 : If the scalar-valued element is used, all nodes remain point evaluation. If the product element is padded into a vector-valued H(curl) element, the nodes on (2, 0) facets become evaluation of the normal component, and interior nodes become evaluation of the ‘vertical’ component of the vector-valued function. H(div)/H(curl) × L2 : Assume the product element, taking values in R2 is padded into an H(div) element. The nodes on (1, 1) facets become evaluation of the normal component. The interior nodes become evaluation of the two ‘horizontal’ components of the vector-valued function. L2 × L2 : 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. 2.5.6 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 PA ⊗ PB , as described in subsubsection 2.5.2, but with no manipulation (e.g., padding with zeros). The basis for PA ⊗ PB 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.

18

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. 2.6 Product finite elements within finite element exterior calculus

Arnold et al. (2014) 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 Qr and QDG r , and various well-known vector-valued spaces as introduced in Brezzi et al. (1985); N´ed´elec (1980, 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. 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 d dx U0 −→ U1 ,

(2.26)

where U0 ⊂ H 1 and U1 ⊂ L2 . In two dimensions, there were two types of complex: ∇⊥

∇·

U0 −→ U1 −→ U2 ,

(2.27)

where U0 ⊂ H 1 , U1 ⊂ H(div) and U2 ⊂ L2 , and ∇⊥ ·

∇

U0 −→ U1 −→ U2 ,

(2.28)

where U0 ⊂ H 1 , U1 ⊂ H(curl) and U2 ⊂ L2 . As mentioned previously, this arises due to two possible identifications of differential 1-forms with vector fields. In three dimensions, the complex takes the form ∇

∇×

∇·

U0 −→ U1 −→ U2 −→ U3 ,

(2.29)

where U0 ⊂ H 1 , U1 ⊂ H(curl), U2 ⊂ H(div) and U3 ⊂ L2 . If we are given two one-dimensional complexes (U0 ,U1 ) and (V0 ,V1 ), we can generate two possible product complexes: ∇⊥

∇·

W0 −→ W1 −→ W2 ,

19

(2.30)

where W0 := U0 ⊗V0 , W1 := HDiv(U0 ⊗V1 ) ⊕ HDiv(U1 ⊗V0 ),

(2.31)

W2 := U1 ⊗V1 ,

(2.33)

(2.32)

with W0 ⊂ H 1 ,W1 ⊂ H(div),W2 ⊂ L2 (compare the complex (2.27)), and ∇⊥ ·

∇

W0 −→ W1 −→ W2 ,

(2.34)

W0 := U0 ⊗V0 , W1 := HCurl(U0 ⊗V1 ) ⊕ HCurl(U1 ⊗V0 ),

(2.35)

W2 := U1 ⊗V1 ,

(2.37)

where (2.36)

with W0 ⊂ H 1 ,W1 ⊂ H(curl),W2 ⊂ L2 (compare the complex (2.28)). 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(U0 ⊗V1 ) has nodes on (0, 1) facets, while HDiv(U1 ⊗V0 ) has nodes on (1, 0) facets. Any interior nodes of HDiv(U0 ⊗V1 ) are evaluation of the first component of the vectorvalued function at certain points, while any interior nodes of HDiv(U1 ⊗ V0 ) are evaluation of the second component. A similar statement can be made for the H(curl) element. If we are given a two-dimensional complex (U0 ,U1 ,U2 ) and a one-dimensional complex (V0 ,V1 ), we can generate a product complex in three dimensions: ∇

∇×

∇·

W0 −→ W1 −→ W2 −→ W3 ,

(2.38)

W0 := U0 ⊗V0 , W1 := HCurl(U0 ⊗V1 ) ⊕ HCurl(U1 ⊗V0 ),

(2.39)

W2 := HDiv(U1 ⊗V1 ) ⊕ HDiv(U2 ⊗V0 ), W3 := U2 ⊗V1 ,

(2.41)

where (2.40) (2.42)

with W0 ⊂ H 1 ,W1 ⊂ H(curl),W2 ⊂ H(div),W3 ⊂ L2 (compare the complex (2.29)). 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(U0 ⊗V1 ) has nodes on (0, 1) edges, while HCurl(U1 ⊗V0 ) has nodes on (1, 0) edges. HCurl(U0 ⊗V1 ) may have nodes on (1, 1) facets; if so, these are evaluation of one of the two tangential components. HCurl(U1 ⊗V0 ) may have nodes on (1, 1) facets; if so, these are evaluation of the other tangential component. HCurl(U1 ⊗ V0 ) may also have nodes on (2, 0) facets; if so, these are evaluation of both tangential components. HCurl(U0 ⊗ V1 ) may have interior nodes; if so, these are evaluation of one of the three vector components. HCurl(U1 ⊗ V0 ) 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.

20

2.7 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 q1 q2 qm x1p1 x2p2 . . . xnpn · xn+1 xn+2 . . . xn+m ,

n

m

∑ pi ≤ p, ∑ qi ≤ q

i=1

(2.43)

i=1

rather than a more expensive quadrature rule capable of integrating n

qm q2 q1 . . . xn+m , xn+2 x1p1 x2p2 . . . xnpn · xn+1

m

∑ pi + ∑ qi ≤ p + q.

i=1

(2.44)

i=1

Suppose we are given existing quadrature rules capable of integrating a polynomial f of degree p over KA ⊂ Rn , and of integrating a polynomial g of degree q over KB ⊂ Rm , say Z KA

Np

f (x1 , . . . , xn ) dx = ∑

(A) (A) wi f (ξi ),

i=1

Nq

Z KB

(B)

(B)

g(x1 , . . . , xm ) dx = ∑ wi g(ξi ),

(2.45)

i=1

where the wi and ξi are quadrature weights and points, respectively. Then a quadrature rule capable of integrating a polynomial h defined on KA × KB ⊂ Rn+m with maximal degree p in x1 , . . . , xn and q in xn+1 , . . . , xn+m , as in (2.43), is Z KA ×KB

h(x1 , . . . , xn+m ) dx = ∑ wi, j h(ξi, j ),

(2.46)

1≤i≤N p 1≤ j≤Nq

where the quadrature weights are (A) (B)

wi, j := wi w j

(2.47)

and the quadrature points are   (A) (A) (B) (B) ξi, j := (ξi )1 , . . . , (ξi )n , (ξ j )1 , . . . , (ξ j )m , (A)

a concatenation of the coordinates of ξi

(2.48)

(B)

and ξ j .

2.8 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. (2006, 2010, 2014). The following definitions can be found in Arnold et al. (2014), which contains a far more thorough overview than is given here.

21

In Rn , 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. Let KA ⊂ Rn , KB ⊂ Rm be domains. Given a de Rham subcomplex on KA , d

d

d

U0 −→ U1 −→ · · · −→ Un ,

(2.49)

where each Uk is a space of differential k-forms on KA , and a de Rham subcomplex on KB , d

d

d

V0 −→ V1 −→ · · · −→ Vm ,

(2.50)

where each Vk is a space of differential k-forms on KB , the tensor product of the complexes (2.49) and (2.50) is a de Rham subcomplex on KA × KB : d

d

d

(U ⊗V )0 −→ (U ⊗V )1 −→ · · · −→ (U ⊗V )n+m ,

(2.51)

where (U ⊗V )k := (U0 ⊗Vk ) ⊕ (U1 ⊗Vk−1 ) ⊕ · · · ⊕ (Uk ⊗V0 ) ≡

M

(Ui ⊗V j ),

(2.52)

i+ j=k

for k = 0, 1, . . . , n + m. Note that (U ⊗V )k is a space of k-forms on KA ⊗ KB , 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.

3 Implementation The product elements, constructed in the previous section, have been implemented within Firedrake (Rathgeber et al., 2014; Rathgeber, 2014). 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):

22

FIAT FInite element Automatic Tabulator (Kirby, 2004, 2012), for the construction and tabula-

tion of finite element basis functions UFL Unified Form Language (Alnæs et al., 2014; Alnæs, 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 syn-

tax trees generated by FFC to improve instruction-level parallelism PyOP2 (Rathgeber et al., 2012), a framework for carrying out parallel computations on unstruc-

tured 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. 3.1 FIAT

FIAT, the FInite element Automatic Tabulator (Kirby, 2004, 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 Pk . 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

23

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. 3.2 UFL

UFL, the Unified Form Language (Alnæs 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, V1 is referred to as the test space. For ρ = 2, V2 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 {Wk }nk=1 and argument spaces: a : W1 ×W2 × . . .Wn ×Vρ × . . . ×V1 → R;

(3.1)

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.

24

3.3 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. 3.4 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 6= 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

25

Listing 1: FIAT representation of a triangular reference cell triangle_topology = {0: {0: [0], 1: [1], 2: [2]}, 1: {0: [1, 2], 1: [0, 2], 2: [0, 1]}, 2: {0: [0, 1, 2]}}

Listing 2: FIAT representation of an interval reference cell interval_topology = {0: {0: [0], 1: [1]}, 1: {0: [0, 1]}}

equal if their constituent cells compare equal and, since gdim can be overridden at creation, their geometric dimensions compare equal. 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 K1 and K2 with n1 and n2 vertices, respectively. Assuming existing orderings for the vertices of K1 and K2 , with the numbering beginning from 0, there is a natural ordering for the vertices of the product cell: vertex i of K1 and vertex j of K2 are concatenated to produce vertex i · n2 + j of K1 × K2 . This naturally runs from 0 to (n1 n2 − 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.

26

v2

e0

e1 c0

v0

e2

v1

Figure 1: FIAT reference triangle with topological entities labelled consistent with Listing 1. The letters v, e and c denote vertices, edges and cells, respectively. These correspond to entities of dimension 0, 1 and 2. v0

c0

v1

Figure 2: FIAT reference interval with topological entities labelled consistent with Listing 2. 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 product_topology = {} for dim_A in A_topology: for dim_B in B_topology: product_topology[(dim_A, dim_B)] = {} tmp = 0 for idx_A in A_topology[dim_A]: for idx_B in B_topology[dim_B]: product_topology[(dim_A, dim_B)][tmp] = \ [i*nvert_B + j for i in A_topology[dim_A][idx_A] for j in B_topology[dim_B][idx_B]] tmp += 1

Listing 4: FIAT representation of a triangular prism product reference cell product_topology = {(0, 0): {0: 4: (0, 1): {0: (1, 0): {0: 3: (1, 1): {0: 2: (2, 0): {0: (2, 1): {0:

[0], 1: [1], 2: [2], 3: [3], [4], 5: [5]}, [0, 1], 1: [2, 3], 2: [4, 5]}, [2, 4], 1: [3, 5], 2: [0, 4], [1, 5], 4: [0, 2], 5: [1, 3]}, [2, 3, 4, 5], 1: [0, 1, 4, 5], [0, 1, 2, 3]}, [0, 2, 4], 1: [1, 3, 5]}, [0, 1, 2, 3, 4, 5]}}

27

v5 v1

(1,0)

(1,0)

e1

e3

(2,0)

v3

f1

(1,0)

e5 (0,1)

e2

(1,1)

(0,1) e0

(0,1) e1

f0

(1,1)

f1

(1,1)

f2 (1,0) e2

v4

v0

v2

(1,0)

e0 (1,0)

e4

(2,0)

f0

Figure 3: FIAT reference triangular prism product cell with topological entities labelled consistent with Listing 4. For clarity, the vertices, edges and faces are shown on separate images. Left: the 6 vertices are labelled. Middle: the 9 edges are labelled; note that the 3 ‘vertical’ (0, 1) edges and the 6 ‘horizontal’ (1, 0) edges are treated separately. Right: the 5 faces are labelled; note that the 3 ‘vertical’ (1, 1) faces and the 2 ‘horizontal’ (2, 0) faces are treated separately. 3.5 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

28

Listing 5: Algorithm to produce the geometric decomposition of product nodes in FIAT product_nodes = {} for dim_A in A_nodes: for dim_B in B_nodes: product_nodes[(dim_A, dim_B)] = {} tmp = 0 for idx_A in A_nodes[dim_A]: for idx_B in B_nodes[dim_B]: product_nodes[(dim_A, dim_B)][tmp] = \ [i*nnodes_B + j for i in A_nodes[dim_A][idx_A] for j in B_nodes[dim_B][idx_B]] tmp += 1

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.

29

Listing 6: Construction of a complicated product complex in UFL U0_0 U0_1 U0 = U1 = U2 =

= FiniteElement("P", triangle, 2) = FiniteElement("B", triangle, 3) EnrichedElement(U0_0, U0_1) FiniteElement("BDFM", triangle, 2) FiniteElement("DP", triangle, 1)

V0 = FiniteElement("P", interval, 1) V1 = FiniteElement("DP", interval, 0) W0 = W1_h W1_v W1 = W2_h W2_v W2 = W3 =

OuterProductElement(U0, V0) = HCurl(OuterProductElement(U1, V0)) = HCurl(OuterProductElement(U0, V1)) EnrichedElement(W1_h, W1_v) = HDiv(OuterProductElement(U1, V1)) = HDiv(OuterProductElement(U2, V0)) EnrichedElement(W2_h, W2_v) OuterProductElement(U2, V1)

3.6 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): ∇⊥

∇·

P2 ⊕ B3 −→ BDFM2 −→ DP1 .

(3.2)

Suppose we wish to take the product of this with some complex on intervals, such as d dx DP0 . P1 −→

(3.3)

This generates a complex on triangular prisms: ∇

∇×

∇·

W0 −→ W1 −→ W2 −→ W3 ,

(3.4)

W0 := (P2 ⊕ B3 ) ⊗ P1 , W1 := HCurl((P2 ⊕ B3 ) ⊗ P1 ) ⊕ HCurl(BDFM2 ⊗ P1 ),

(3.5)

W2 := HDiv(BDFM2 ⊗ DP0 ) ⊕ HDiv(DP1 ⊗ P1 ), W3 := DP1 ⊗ DP0 .

(3.7)

where (3.6) (3.8)

We see that W1 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 P2 ⊗ (DP1 ⊗ P2 ), which is necessary to produce elements on hexahedra.

30

(3.9)

3.7 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. 3.8 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).

31

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. 3.9 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. 3.10 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. 3.10.1 Quadrilaterals

Consider the following complex on intervals: d dx Pn −→ DPn−1 .

(3.10)

Taking the product of (3.10) with itself produces the following complexes of well-known elements on quadrilaterals: ∇⊥

∇·

∇

∇⊥ ·

Qn −→ RTFn −→ DQn−1 , or

Qn −→ RTEn −→ DQn−1 ;

(3.11) (3.12)

Qn represents the continuous, scalar-valued element consisting of polynomials of degree at most n in each variable separately, DQn (alternatively QDG n ) is the discontinuous counterpart

32

Listing 7: Construction of the continuous bicubic element Q3 on quadrilaterals U = FiniteElement("P", interval, 3) Q3_quad = OuterProductElement(U, U)

Listing 8: Construction of the div-conforming element RTF3 on quadrilaterals U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 2) RTF3_quad_x = HDiv(OuterProductElement(U, V)) RTF3_quad_y = HDiv(OuterProductElement(V, U)) RTF3_quad = EnrichedElement(RTF3_quad_x, RTF3_quad_y)

to Qn , and RTF and RTE represent the Raviart–Thomas families of ‘face’ and ‘edge’ elements on quadrilaterals (Raviart and Thomas, 1977). Constructions of these elements are given in Listing 7–Listing 10. 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´ed´elec second-kind’ elements introduced in N´ed´elec (1986), and so we will use the name Ndtwo. The construction is given in Listing 11. 3.10.2 Hexahedra

Although Firedrake does not yet have full support for meshes of hexahedral cells, the appropriate finite element families can still be expressed in UFL and tabulated in FIAT. Taking the product of three copies of the one-dimensional complex given in (3.10), or the product of either (3.11)

Listing 9: Construction of the curl-conforming element RTE3 on quadrilaterals U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 2) RTE3_quad_x = HCurl(OuterProductElement(V, U)) RTE3_quad_y = HCurl(OuterProductElement(U, V)) RTE3_quad = EnrichedElement(RTE3_quad_x, RTE3_quad_y)

33

Listing 10: Construction of the discontinuous biquadratic element DQ2 on quadrilaterals V = FiniteElement("DP", interval, 2) DQ2_quad = OuterProductElement(V, V)

Listing 11: Construction of the div-conforming bicubic element N2F3 on quadrilaterals U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 3) N2F3_quad_x = HDiv(OuterProductElement(U, V)) N2F3_quad_y = HDiv(OuterProductElement(V, U)) N2F3_quad = EnrichedElement(N2F3_quad_x, N2F3_quad_y)

or (3.12) with (3.10), gives a complex of well-known elements on hexahedra: ∇

∇×

∇·

Qn −→ N1En −→ N1Fn −→ DQn−1 .

(3.13)

Here, Qn represents the continuous, scalar-valued element consisting of polynomials of degree up to n in each of the three coordinates separately, while DQn would be the discontinuous counterpart. N1E and N1F represent curl- and div-conforming elements which are similar to the “N´ed´elec elements of the first kind” on tetrahedra that were introduced in N´ed´elec (1980). Constructions of these elements are given in Listing 12–Listing 15. The “N´ed´elec elements of the second kind” (N´ed´elec, 1986) contain all polynomials up to a given degree in each variable separately. Constructions are given in Listing 16 and Listing 17. 3.10.3 Triangular prisms

The RT complex on triangles is ∇⊥

∇·

∇

∇⊥ ·

Pn −→ RTFn −→ DPn−1 or

Pn −→ RTEn −→ DPn−1 .

(3.14) (3.15)

Taking the product of (3.14) or (3.15) with the interval complex (3.10) gives a complex on triangular prisms: ∇ ∇× ∇· Qn −→ N1En −→ N1Fn −→ DQn−1 . (3.16) Listing 12: Construction of the continuous tricubic element Q3 on hexahedra U = FiniteElement("P", interval, 3) Q3_quad = OuterProductElement(U, U) Q3_hex = OuterProductElement(Q3_quad, U)

34

Listing 13: Construction of the curl-conforming element N1E3 on hexahedra U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 2) Q3_quad = OuterProductElement(U, U) RTE3_quad_x = HCurl(OuterProductElement(V, U)) RTE3_quad_y = HCurl(OuterProductElement(U, V)) N1E3_hex_x = HCurl(OuterProductElement(RTE3_quad_x, U)) N1E3_hex_y = HCurl(OuterProductElement(RTE3_quad_y, U)) N1E3_hex_z = HCurl(OuterProductElement(Q3_quad, V)) N1E3_hex_xy = EnrichedElement(N1E3_hex_x, N1E3_hex_y) N1E3_hex = EnrichedElement(N1E3_hex_xy, N1E3_hex_z)

Listing 14: Construction of the div-conforming element N1F3 on hexahedra U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 2) RTF3_quad_x = HDiv(OuterProductElement(U, V)) RTF3_quad_y = HDiv(OuterProductElement(V, U)) DQ2_quad = OuterProductElement(V, V) N1F3_hex_x = HDiv(OuterProductElement(RTF3_quad_x, V)) N1F3_hex_y = HDiv(OuterProductElement(RTF3_quad_y, V)) N1F3_hex_z = HDiv(OuterProductElement(DQ2_quad, U)) N1F3_hex_xy = EnrichedElement(N1F3_hex_x, N1F3_hex_y) N1F3_hex = EnrichedElement(N1F3_hex_xy, N1F3_hex_z)

Listing 15: Construction of the discontinuous triquadratic element DQ2 on hexahedra V = FiniteElement("DP", interval, 2) DQ2_quad = OuterProductElement(V, V) DQ2_hex = OuterProductElement(DQ2_quad, V)

35

Listing 16: Construction of the curl-conforming tricubic element N2E3 on hexahedra U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 3) Q3_quad = OuterProductElement(U, U) N2E3_quad_x = HCurl(OuterProductElement(V, U)) N2E3_quad_y = HCurl(OuterProductElement(U, V)) N2E3_hex_x = HCurl(OuterProductElement(N2E3_quad_x, U)) N2E3_hex_y = HCurl(OuterProductElement(N2E3_quad_y, U)) N2E3_hex_z = HCurl(OuterProductElement(Q3_quad, V)) N2E3_hex_xy = EnrichedElement(N2E3_hex_x, N2E3_hex_y) N2E3_hex = EnrichedElement(N2E3_hex_xy, N2E3_hex_z)

Listing 17: Construction of the div-conforming tricubic element N2F3 on hexahedra U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 3) RTF3_quad_x = HDiv(OuterProductElement(U, V)) RTF3_quad_y = HDiv(OuterProductElement(V, U)) DQ3_quad = OuterProductElement(V, V) N2F3_hex_x = HDiv(OuterProductElement(RTF3_quad_x, V)) N2F3_hex_y = HDiv(OuterProductElement(RTF3_quad_y, V)) N2F3_hex_z = HDiv(OuterProductElement(DQ3_quad, U)) N2F3_hex_xy = EnrichedElement(N2F3_hex_x, N2F3_hex_y) N2F3_hex = EnrichedElement(N2F3_hex_xy, N2F3_hex_z)

36

Listing 18: Construction of the continuous element Q3 on triangular prisms P3_tri = FiniteElement("P", triangle, 3) U = FiniteElement("P", interval, 3) Q3_prism = OuterProductElement(P3_tri, U)

Listing 19: Construction of the curl-conforming element N1E3 on triangular prisms P3_tri = FiniteElement("P", triangle, 3) RTE3_tri = FiniteElement("RTE", triangle, 3) U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 2) N1E3_prism_xy = HCurl(OuterProductElement(RTE3_tri, U)) N1E3_prism_z = HCurl(OuterProductElement(P3_tri, V)) N1E3_prism = EnrichedElement(N1E3_prism_xy, N1E3_prism_z)

Here, we have used Qn to denote the continuous, scalar-valued element consisting of polynomials at most degree n in the first two spatial variables and the third spatial variable, separately. DQn represents the discontinuous counterpart, while N1E and N1F represent N´ed´elec elements of the first kind, analogous to the tetrahedral elements introduced in N´ed´elec (1980). Constructions of these elements are given in Listing 18–Listing 21. As for hexahedra, there are also N´ed´elec elements of the second kind (N´ed´elec, 1986); constructions are given in Listing 22 and Listing 23.

4 Numerical examples In this section, we give several examples to demonstrate both the range of functionality and the correctness of our implementation. Quantitative analysis is performed where possible, e.g. demonstration of convergence to a known solution at expected order with increasing mesh reso-

Listing 20: Construction of the div-conforming element N1F3 on triangular prisms RTF3_tri = FiniteElement("RTF", triangle, 3) DP2_tri = FiniteElement("DP", triangle, 2) U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 2) N1F3_prism_xy = HDiv(OuterProductElement(RTF3_tri, V)) N1F3_prism_z = HDiv(OuterProductElement(DP2_tri, U)) N1F3_prism = EnrichedElement(N1F3_prism_xy, N1F3_prism_z)

37

Listing 21: Construction of the discontinuous triquadratic element DQ2 on triangular prisms DP2_tri = FiniteElement("DP", triangle, 2) V = FiniteElement("DP", interval, 2) DQ2_prism = OuterProductElement(DP2_tri, V)

Listing 22: Construction of the curl-conforming element N2E3 on triangular prisms P3_tri = FiniteElement("P", triangle, 3) BDME3_tri = FiniteElement("BDME", triangle, 3) U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 3) N2E3_prism_xy = HCurl(OuterProductElement(BDME3_tri, U)) N2E3_prism_z = HCurl(OuterProductElement(P3_tri, V)) N2E3_prism = EnrichedElement(N2E3_prism_xy, N2E3_prism_z)

Listing 23: Construction of the div-conforming element N2F3 on triangular prisms BDMF3_tri = FiniteElement("BDMF", triangle, 3) DP3_tri = FiniteElement("DP", triangle, 3) U = FiniteElement("P", interval, 3) V = FiniteElement("DP", interval, 3) N2F3_prism_xy = HDiv(OuterProductElement(BDMF3_tri, V)) N2F3_prism_z = HDiv(OuterProductElement(DP3_tri, U)) N2F3_prism = EnrichedElement(N2F3_prism_xy, N2F3_prism_z)

38

lution. 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. 4.1 Vector Laplacian

We seek a solution to −∇(∇ ·~u) + ∇ × (∇ ×~u) = ~f

(4.1)

in a domain Ω, with boundary conditions ~u ·~n = 0,

(4.2)

(∇ ×~u) ×~n = 0

(4.3)

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)): σ = −∇ ·~u, ∇σ + ∇ × (∇ ×~u) = ~f .

(4.4) (4.5)

Let V0 ⊂ H 1 , V1 ⊂ H(curl) be finite element spaces. A suitable formulation of (4.1), written in weak form, is then: find σ ∈ V0 , ~u ∈ V1 such that hτ, σ i − h∇τ,~ui = 0,

(4.6)

h~v, ∇σ i + h∇ ×~v, ∇ ×~ui = h~v, ~f i,

(4.7)

for all τ ∈ V0 ,~v ∈ V1 . For brevity, we have used angled brackets to denote the standard L2 inner product: Z Z ha, bi := ab dx, h~a,~bi := ~a ·~b dx. (4.8) Ω

Ω

The boundary conditions have been implicitly applied, in a weak sense, through neglecting the surface terms when integrating by parts. 4.1.1 2D

We take Ω to be the unit square [0, 1]2 . The choice   ~f = π 2 (k2 + l 2 ) sin(kπx) cos(lπy) cos(kπx) sin(lπy)

(4.9)

produces the solution   sin(kπx) cos(lπy) ~u = , cos(kπx) sin(lπy)

39

(4.10)

101 100 10-1 L2 norm of error

10-2 10-3 10-4 10-5 10

-6

10-7 10-8 10-9 0 10

Q1 σ RTE1 u Q2 σ RTE2 u Q3 σ RTE3 u 10-1

∆x

10-2

10-3

Figure 4: The L2 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 ∆xn , for n from 1 to 4, and are merely to aid comprehension. The convergence rates are as expected: the Q1 –RTE1 case demonstrates second-order convergence for ~u and first-order convergence for σ . Recall that, in our numbering convention, RTE1 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 Q2 –RTE2 and Q3 –RTE3 converge at correspondingly faster rates. which satisfies the boundary conditions (4.2) and (4.3). To discretise this problem, we subdivide Ω into squares with side length ∆x. We use Qn for the H 1 space, and RTEn 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 Qn+1 . The L2 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 Qn+1 .

40

4.1.2 3D

We take Ω to be the unit cube [0, 1]3 . The choice  2 2  (k + l ) sin(kπx) cos(lπy) ~f = π 2  (l 2 + m2 ) sin(lπy) cos(mπz)  (k2 + m2 ) sin(mπz) cos(kπx)

(4.11)

produces the solution 

 sin(kπx) cos(lπy) ~u =  sin(lπy) cos(mπz)  sin(mπz) cos(kπx)

(4.12)

which satisfies the boundary conditions (4.2) and (4.3). 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 Qn prism element for the H 1 space, and the N1En 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 Qn+1 . The L2 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 Qn+1 . 4.2 Gravity wave

A simple model of the atmosphere is given by ∂~u = −∇p + bˆz, ∂t ∂b = −N 2~u · zˆ, ∂t ∂p = −c2 ∇ ·~u, ∂t

(4.13) (4.14) (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¨ais¨al¨a frequency and the speed of sound respectively, while zˆ 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 ∇

∇×

∇·

W0 −→ W1 −→ W2 −→ W3 ,

41

(4.17)

101

L2 norm of error

100 10-1 10-2 10-3 10-4 0 10

Q1 σ N1E1 u Q2 σ N1E2 u Q3 σ N1E3 u ∆x

10-1

Figure 5: The L2 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 ∆xn , 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.

42

as in (2.38), we can seek a solution with ~u ∈ W20 , b ∈ W2v and p ∈ W3 . W20 is the subspace of W2 whose normal component vanishes on the boundary of the domain – this enforces the boundary condition (4.16). We have used W2v to represent the “vertical” part of W2 : if we write W2 as a sum of two product elements HDiv(U1 ⊗ V1 ) and HDiv(U2 ⊗ V0 ) then W2v is the scalar-valued product U2 ⊗ V0 . 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: find ~u ∈ W20 , b ∈ W2v , p ∈ W3 such that   ∂~u ~w, − h∇ · ~w, pi − h~w, bˆzi = 0 (4.18) ∂t   ∂b γ, + N 2 hγ,~u · zˆi = 0 (4.19) ∂t   ∂p φ, + c2 hφ , ∇ ·~ui = 0 (4.20) ∂t for all ~w ∈ W20 , γ ∈ W2v , φ ∈ W3 . 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 1 2 1 1 |~u| + 2 b2 + 2 p2 dx. 2N 2c Ω2

Z

(4.21)

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. 4.2.1 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¨ais¨al¨a frequency of 10−2 s−1 and a sound speed of 300ms−1 . We start the simulation from rest, with a buoyancy perturbation b=

sin(πy/H) 1 + (x − x0 )2 /a2

(4.22)

together with a vertically balancing pressure field p=−

H cos(πy/H) , π 1 + (x − x0 )2 /a2

43

(4.23)

1e12 4.0 3.5 3.0 energy

2.5 2.0 1.5

PE KE IE total

1.0 0.5 0.0 0

20000

40000

60000 t

80000

100000

120000

Figure 6: Evolution of energy for the 2D simulation described in subsubsection 4.2.1. The components are the potential energy, PE, the kinetic energy, KE, and the internal energy, IE. The choice of spatial and temporal discretisations leads to exact conservation of total energy up to solver tolerances; this is indeed observed. The event at approximately t = 95,000s corresponds to the two gravity waves briefly recombining, as the domain is periodic. where H is the height of the domain, a is a horizontal length-scale, which we take to be 100km, and x0 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 W2 to be RTF2 , and W3 to be DQ1 . 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 × 106 , which is comparable to the linear solver tolerances. 4.2.2 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¨ais¨al¨a frequency of 10−2 s−1 and a sound speed of 300ms−1 . We start the simulation from rest, with a buoyancy

44

1e20 7 6

energy

5 4 3 2

PE KE IE total

1 0 0

100000

200000

t

300000

400000

500000

Figure 7: Evolution of energy for the 3D simulation described in subsubsection 4.2.2. The components are the potential energy, PE, the kinetic energy, KE, and the internal energy, IE. The choice of spatial and temporal discretisations leads to exact conservation of total energy up to solver tolerances; this is indeed observed. The event at approximately t = 320,000s corresponds to the zonally-symmetric gravity wave reaching the poles of the spherical domain. perturbation b=

sin(π(|~x| − a)/H) 1 + z2 /L2

(4.24)

together with a vertically balancing pressure field p=−

H cos(π(|~x| − a)/H) , π 1 + z2 /L2

(4.25)

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 BDFM2 complex on triangles (2.13) and the P2 –DP1 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 × 108 , which is comparable to the linear solver tolerances.

45

4.3 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 + ∇ · (~ u0 q) = 0. (4.26) ∂t Let φ be an arbitrary smooth test-function. Multiplying (4.26) by φ and integrating over the domain Ω gives Z Z ∂q dx + φ ∇ · (~ u0 q) dx = 0. (4.27) φ Ω Ω ∂t 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 ∇ · (~ u0 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: ∂q dx + φ ∇ · (~ u0 q) dx = 0, ∀ φ ∈ V. (4.28) e ∂t e We now apply integration by parts to move the derivative onto the test-functions φ , which are only non-zero on a single cell: Z

Z

φ

∂q φ dx − q~ u0 · ∇φ dx + φ qeu~0 ·~n dS = 0 ∀ φ ∈ V, (4.29) e ∂t e ∂e 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 qe. Summing (4.29) over all cells in the mesh then gives   Z Z ∂q = h∇φ , q~ u0 i − φ, φ qeu~0 ·~n ds − (φ+ − φ− )e qu~0 ·~n dS ∀ φ ∈ V, (4.30) ∂t Γext Γint Z

Z

Z

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. We will discretise (4.30) in time using a three-stage strong-stability-preserving Runge-Kutta scheme (Shu and Osher, 1988):   ∆q(1) = L q(n) (4.31)   ∆q(2) = L q(n) + ∆q(1) (4.32)   1 1 (4.33) ∆q(3) = L q(n) + ∆q(0) + ∆q(1) 4 4   1 (1) 1 (2) 2 (3) (n+1) (n) ∆q + ∆q + ∆q , (4.34) q = q + ∆t 6 6 3 where L is the appropriate time-evolution operator.

46

4.3.1 2D

We take Ω to be the unit square [0, 1]2 . Our initial condition will be a cosine hill   (  |~x−~ x0 | 1 1 + cos π , |~x − x~0 | < r0 r0 q= 2 0, otherwise, with radius r0 = 0.15, centred at x~0 = (0.25, 0.5). The prescribed velocity field is  πt   sin(πx)2 sin(2πy)  , u~0 (~x,t) = cos − sin(πy)2 sin(2πx) T

(4.35)

(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 = 23 . To discretise this problem, we subdivide Ω into squares with side length ∆x. We use DQn 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 Q2 . The L2 errors between the initial and final q fields for varying ∆x are plotted in Figure 8. 4.3.2 3D

We take Ω to be the unit cube [0, 1]3 . We use the same, z-independent, cosine hill for the initial condition:   √   p  1 1 + cos π (x−x0 )2 +(y−y0 )2 , (x − x0 )2 + (y − y0 )2 < r0 r0 q= 2 (4.37)  0, otherwise, where r0 = 0.15, x0 = 0.25 and y0 = 0.5. The prescribed velocity field is    πt  sin(πx)2 sin(2πy) sin(πz) − sin(πy)2 sin(2πx) sin(πz) , u~0 (~x,t) = cos T 0

(4.38)

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 = 14 . 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 DQn 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 Q2 . The L2 errors between the initial and final q fields for varying ∆x are plotted in Figure 9.

47

10-1

DQ0 DQ1 DQ2

L2 norm of error

10-2

10-3

10-4

10-5 -1 10

10-2

∆x

10-3

10-4

Figure 8: The L2 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 ∆x2 , and are merely to aid comprehension. The DQ0 simulations converge at first-order for sufficiently small values of ∆x. The DQ1 simulations converge at second-order, as expected. The cosine bell initial condition has a discontinuous second derivative, which inhibits the DG2 simulations from exceeding a second-order rate of convergence.

48

10-1

DQ0 DQ1 DQ2

L2 norm of error

10-2

10-3

10-4

10-5 10-1

10-2 ∆x

10-3

Figure 9: The L2 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 ∆x2 , and are merely to aid comprehension. The DQ0 simulations converge at firstorder for sufficiently small ∆x. The DQ1 simulations converge at second-order, as expected. The cosine bell initial condition has a discontinuous second derivative, which inhibits the DG2 simulations from exceeding a second-order rate of convergence.

49

5 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.

6 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

50

implementation of extruded meshes is in preparation.

Acknowledgements Andrew McRae wishes to acknowledge funding and other support from the Grantham Institute and Climate-KIC. This work was supported by the National Environmental Research Council [grant numbers NE/K006789/1, NE/K008951/1]; and an Engineering and Physical Sciences Research Council prize studentship.

References Martin S. Alnæs, Anders Logg, Kristian B. Ølgaard, Marie E. Rognes, and Garth N. Wells. Unified Form Language: A domain-specific language for weak formulations of partial differential equations. ACM Transactions on Mathematical Software, 40(2):9:1–9:37, 2014. ISSN 0098-3500. doi: 10.1145/2566630. Martin Sandve Alnæs. UFL: a finite element form language. In Automated Solution of Differential Equations by the Finite Element Method, volume 84 of Lecture Notes in Computational Science and Engineering, pages 303–338. Springer, 2012. doi: 10.1007/978-3-642-23099-8 17. Douglas N. Arnold and Gerard Awanou. Finite element differential forms on cubical meshes. Mathematics of Computation, 83(288):1551–1570, 2014. ISSN 0025-5718. doi: 10.1090/ S0025-5718-2013-02783-4. Douglas N. Arnold and Anders Logg. Periodic table of the finite elements. SIAM News, November 2014. URL http://femtable.org. Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. Finite element exterior calculus, homological techniques, and applications. Acta Numerica, 15:1–155, 2006. ISSN 1474-0508. doi: 10.1017/S0962492906210018. Douglas N. Arnold, Richard S. Falk, and Ragnar Winther. Finite element exterior calculus: from Hodge theory to numerical stability. Bulletin (New Series) of the American Mathematical Society, 47(2):281–354, 2010. ISSN 0273-0979. doi: 10.1090/S0273-0979-10-01278-4. Douglas N. Arnold, Daniele Boffi, and Francesca Bonizzoni. Finite element differential forms on curvilinear cubic meshes and their approximation properties. Numerische Mathematik, pages 1–20, 2014. ISSN 0029-599X. doi: 10.1007/s00211-014-0631-3. Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, and Hong Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.5, Argonne National Laboratory, 2014. URL http:// www.mcs.anl.gov/petsc.

51

Daniele Boffi, Michel Fortin, and Franco Brezzi. Mixed finite element methods and applications. Springer series in computational mathematics. Springer, Berlin, Heidelberg, 2013. ISBN 9783-642-36518-8. Franco Brezzi and Michel Fortin. Mixed and Hybrid Finite Element Methods. Springer Series in Computational Mathematics. Springer-Verlag, New York, 1991. ISBN 0-387-97582-9. Franco Brezzi, Jim Douglas, Jr., and L. D. Marini. Two families of mixed finite elements for second order elliptic problems. Numerische Mathematik, 47(2):217–235, 1985. ISSN 0029599X. doi: 10.1007/BF01389710. J. G. Charney and N. A. Phillips. Numerical integration of the quasi-geostrophic equations for barotropic and simple baroclinic flows. Journal of Meteorology, 10(2):71–99, 1953. ISSN 0095-9634. doi: 10.1175/1520-0469(1953)010h0071:NIOTQGi2.0.CO;2. Philippe G Ciarlet. The finite element method for elliptic problems. North-Holland, 1978. C. J. Cotter and J. Shipton. Mixed finite elements for numerical weather prediction. Journal of Computational Physics, 231:7076–7091, 2012. ISSN 0021-9991. doi: 10.1016/j.jcp.2012.05. 020. Robert C. Kirby. Algorithm 839: FIAT, a new paradigm for computing finite element basis functions. ACM Trans. Math. Softw., 30(4):502–516, 2004. ISSN 0098-3500. doi: 10.1145/ 1039813.1039820. Robert C. Kirby. FIAT: numerical construction of finite element basis functions. In Automated Solution of Differential Equations by the Finite Element Method, volume 84 of Lecture Notes in Computational Science and Engineering, pages 247–255. Springer, 2012. doi: 10.1007/ 978-3-642-23099-8 13. Robert C. Kirby and Anders Logg. A compiler for variational forms. ACM Trans. Math. Softw., 32(3):417–444, 2006. ISSN 0098-3500. doi: 10.1145/1163641.1163644. Benedict Leimkuhler and Sebastian Reich. Simulating Hamiltonian Dynamics, chapter 12. Cambridge University Press, 2005. ISBN 9780521772907. doi: 10.1017/CBO9780511614118. Randall J LeVeque. High-resolution conservative algorithms for advection in incompressible flow. SIAM Journal on Numerical Analysis, 33(2):627–665, 1996. ISSN 0036-1429. doi: 10.1137/0733033. Anders Logg, Kent-Andre Mardal, Garth N. Wells, et al. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012a. ISBN 978-3-642-23098-1. doi: 10.1007/978-3-642-23099-8. Anders Logg, Kristian B Ølgaard, Marie E Rognes, and Garth N Wells. FFC: the FEniCS form compiler. In Automated Solution of Differential Equations by the Finite Element Method, volume 84 of Lecture Notes in Computational Science and Engineering, pages 227–238. Springer, 2012b. doi: 10.1007/978-3-642-23099-8 11.

52

Fabio Luporini, Ana Lucia Varbanescu, Florian Rathgeber, Gheorghe-Teodor Bercea, J. Ramanujam, David A. Ham, and Paul H. J. Kelly. COFFEE: an optimizing compiler for finite element local assembly. In review, ACM Transactions on Architecture and Code Optimization, 2014. J. C. N´ed´elec. Mixed finite elements in R3 . Numerische Mathematik, 35(3):315–341, 1980. ISSN 0029-599X. doi: 10.1007/BF01396415. J. C. N´ed´elec. A new family of mixed finite elements in R3 . Numerische Mathematik, 50(1): 57–81, 1986. ISSN 0029-599X. doi: 10.1007/BF01389668. Florian Rathgeber. Productive and Efficient Computational Science Through Domain-specific Abstractions. PhD thesis, Imperial College London, July 2014. Florian Rathgeber, Graham R. Markall, Lawrence Mitchell, Nicolas Loriant, David A. Ham, Carlo Bertolli, and Paul H.J. Kelly. PyOP2: A high-level framework for performance-portable simulations on unstructured meshes. In High Performance Computing, Networking Storage and Analysis, SC Companion:, pages 1116–1123, Los Alamitos, CA, USA, 2012. IEEE Computer Society. ISBN 978-1-4673-3049-7. doi: 10.1109/SC.Companion.2012.134. Florian Rathgeber et al. Firedrake: automating the finite element method by composing abstractions. In preparation, 2014. P. A. Raviart and J. M. Thomas. A mixed finite element method for 2nd order elliptic problems. In Mathematical aspects of finite element methods, pages 292–315. Springer, 1977. doi: 10.1007/BFb0064470. M. E. Rognes, R. C. Kirby, and A. Logg. Efficient assembly of H(div) and H(curl) conforming finite elements. SIAM Journal on Scientific Computing, 31(6):4130–4151, 2009. doi: 10. 1137/08073901X. M. E. Rognes, D. A. Ham, C. J. Cotter, and A. T. T. McRae. Automating the solution of pdes on the sphere and other manifolds in FEniCS 1.2. Geoscientific Model Development, 6(6): 2099–2119, 2013. ISSN 1991-9603. doi: 10.5194/gmd-6-2099-2013. Chi-Wang Shu and Stanley Osher. Efficient implementation of essentially non-oscillatory shockcapturing schemes. Journal of Computational Physics, 77(2):439–471, 1988. ISSN 00219991. doi: 10.1016/0021-9991(88)90177-5. W. C. Skamarock and J. B. Klemp. Efficiency and accuracy of the Klemp-Wilhelmson timesplitting technique. Monthly Weather Review, 122(11):2623–2630, 1994. ISSN 0027-0644. doi: 10.1175/1520-0493(1994)122h2623:EAAOTKi2.0.CO;2. St´efan van der Walt, S. Chris Colbert, and Ga¨el Varoquaux. The NumPy array: A structure for efficient numerical computation. Computing in Science & Engineering, 13(2):22–30, 2011. doi: 10.1109/MCSE.2011.37.

53