The idea behind standard regions can be traced back to one of the fundamental principles we learn in calculus: differentiation and integrating on arbitrary regions are often difficult; however, differentiation and integration have been worked out on canonical (straightsided or planarsided), rightangled, coordinate aligned domains. In the case of Nektar++, we do not need to work on truly arbitrary domains, but rather on seven fundamental domains: segments in 1D, triangles and quadrilaterals in 2D, and tetrahedra, hexahedra, prisms and pyramids in 3D. Since Nektar++ deals with polynomial methods, the natural domain over which reference elements should be build is [1,1] as this is the interval over which Jacobi Polynomials are defined and over which Gaussian quadrature is defined [19]. We show in Figure 5.1 a pictorial representation of the standard segment, standard quadrilateral (often shorthanded quad) and standard triangle.
The standard quad and standard hexahedra (shorthanded ’hex’) are geometrically tensorproduct constructions defined on [1,1]^{d} for d = 2 and d = 3 respectively. The standard triangle is constructed by taking ξ_{1} ∈ [1,1] and taking ξ_{2} ≤ξ_{1}. The standard tetrahedra (shorthanded ’tet’) is built upon the standard triangle and has all four faces being triangles, with the two triangles along the coordinate directions looking like the standard triangle. The standard prism consists of a standard triangle along the ξ_{1} ξ_{2} plane extruded into the third direction (yielding three quadrilateral faces). The standard pyramid consists of a standard quadrilateral at the base with four triangular faces reaching up to its top vertex We show this pictorially in Figure 5.2.
Regardless of the particular element type, we use the fact one can build polynomial spaces over these geometric objects. In the case of triangles and tetrahedra, these are polynomial spaces in the mathematical sense, i.e. P(k) spaces. In the case of quadrilaterals and hexahedra, these are bi– and tripolynomial spaces, i.e. Q(k). In the case of prisms and pyramids, the spaces are more complicated and are a mixture in some sense of these spaces, but they are indeed polynomial in form. This allows us to define differentiation exactly and to approximate integrals exactly up to machine precision.
In what is to follow, we describe (1) how you can build differentiation and integration operators over standard regions by mapping them to a domain over which you can exploit index separability; and then (2) how you can build differentiation and integration operators natively over various standard region shapes. The former allows one to benefit from tensorproduct operators, and is at present the primary focus of all Nektar++ operators. The latter as implemented for nodal element types is currently an isa class definition extended from the (tensorproductbased) standard element definitions.
Differentiation and integration over the standard elements within Nektar++, in general, always try to map things back to tensorproduct (e.g. tensorcontraction indexing). This allows one to transform operators that would normally have iterators that go from i = 0,…,N^{d} where d is the dimension of the shape to operators constructed with d⋅N operations (i.e., the product of onedimensional operations).
As an example (to help the reader gain an intuitive understanding), consider the integration of a function f() over a hexahedral element. If we had a quadrature rule in 3D that needed N^{d} points to integrate this function exactly, we would express the operation as:

where E denotes our ddimensional element and the set Q = {_{i},ω_{i}} denotes the set of points and quadrature weights over the element E. Now, suppose that both our function f() and our Q were separable: that is, f() can be written as f_{1}(x_{1}) ⋅ f_{2}(x_{2}) ⋅ f_{3}(x_{3}) where = (x_{1},x_{2},x_{3})^{T } and Q can we written in terms of 1D quadrature: _{i} = z_{i1}^{(1)} ⋅ z_{i2}^{(2)} ⋅ z_{i3}^{(3)} with the index handled by an index map σ defined by i = σ(i_{1},i_{2},i_{3}) = i_{1} + N ⋅ i_{2} + N^{2} ⋅ i_{3}, and similarly for the weights ω_{i}. We can then rewrite the integral above as follows:
As you can see, when the functions are separable and the quadrature (or collocating points) can be written in separable form, we can take O(N^{d}) operators and transform them into O(d ⋅ N)) operators. The above discussion is mainly focusing on the mathematical transformations needed to accomplish this; in subsequent sections, we will point out the memory layout and index ordering (i.e. now σ is ordered and implemented) to gain maximum performance.
The discussion of tensorproduct operations on quadrilaterals and hexahedra may seem quite natural as the element construction is done with tensorproducts of the segment; however, how do you create these types of operations on triangles (and anything involving triangles such as tetrahedra, prisms and pyramids)? The main mathematical building block we use that allows such operators is a transformation due to Duffy [29], which was subsequently used by Dubiner [28] in the context of finite elements and was extended to general polyhedral types by Ainsworth [7, 6].
An image denoting the Duffy transformation is shown in Figure 5.3. On the left is an example of the rightsided standard triangle with coordinate system (ξ_{1},ξ_{2}), and on the right is the separable tensorproduct domain with coordinates (η_{1},η_{2}). We often refer to this transformation as “collapsed coordinates" as it appears as the “collapsing" of a quadrilateral domain to a triangle. Note that the edge along η_{1} = 1 on the quadrilateral collapses down to the single point (ξ_{1},ξ_{2}) = (1,1). This, in general, does not cause issues with our integration over volumes as this is a point of measure zero. However, special case will need be taken when doing edge integrals. We will make it a point to highlight those places in which special care is needed.
With the 2D Duffy transformation in place, we can actually create all four 3D element types from the starting point of a hexahedra (our base tensorproduct shape). One collapsing operation yields a prism. The second collapsing operation yields a pyramid. Lastly, a final collapsing operation yields a tetrahedra. A visual representation of these operations are shown in Figure 5.4.
This is meant merely to be a summary of tensorproduct operations, enough to allow us to discuss the basic data types and algorithms within Nektar++. If more details are needed, we encourage the reader to consult Karniadakis and Sherwin [46] and the references therein.
Although the primary backbone of Nektar++ is tensorproduct operations across all element types through the use of collapsed coordinates, we have implemented two commonly used nontensorial basis set definitions. These are based upon Lagrange polynomial construction from a nodal (collocating) point set. As derived element types (in the C + + sense), we have implemented NodalStdTriangle (and Tet) based upon the electrostatic points of Hesthaven [40] and the Fekete points of Taylor and Wingate [63, 64].
Although these types of elements, at first glance, might seem to be “suboptimal" (in terms of operations), there are various reasons why people might choice to use them. For instance, the point set you use for defining the collocation points dictates the interpolative projection operator (and correspondingly the properties of that operator) – an important issue when evaluating boundary conditions, initial conditions and for some nonlinear operator evaluations. Within the Nektar++ team, we started an examination of this within the paper by Kirby and Sherwin [48]. There has since been various papers in the literature (beyond the scope of this developer’s guide) that give the pros and cons of tensorproduct (separable) and nontensorproduct element types.