5.2 The Fundamental Data Structures within StdRegions

In almost all object-oriented languages (which includes C + +), there exists the concepts of class attributes and object attributes. Class attributes are those attributes shared by all object instances (both immediate and derived objects) of a particular class definition, and object attributes (sometimes called data members) are those attributes whose values vary from object to object and hence help to characterize (or make unique) a particular object. In C + +, object attributes are specified a header file containing class declarations; within a class declaration, attributes are grouped by their accessibility: public attributes, protected attributes and private attributes. A detailed discussion of the nuances of these categories are beyond the scope of this guide; we refer the interested reader to the following books for further details: [5551]. For our purposes, the main thing to appreciate is that categories dictate access patters within the inheritance hierarchy and to the “outside” world (i.e. access from outside the object). We have summarized the relationships between the categories and their accessibility in Tables 5.1, 5.2 and 5.3 1.

Table 5.1: Accessibility in Public Inheritance

Accessibility private variables protected variables public variables

Accessibility from own class? yes yes yes
Accessibility from derived class? no yes yes
Accessibility from 2nd derived class? no yes yes

Table 5.2: Accessibility in Protected Inheritance

Accessibility private variables protected variables public variables

Accessibility from own class? yes yes yes
Accessibility from derived class? no yes yes
(inherited as
protected variable)
Accessibility from 2nd derived class? no yes yes

Table 5.3: Accessibility in Private Inheritance

Accessibility private variables protected variables public variables

Accessibility from own class? yes yes yes
Accessibility from derived class? no yes yes
(inherited as (inherited as
private variable) private variable)
Accessibility from 2nd derived class? no no no

Within the StdRegions directory of the library, there exists a class inheritance hierarchy designed to try to encourage re-use of core algorithms (while simultaneously trying to minimize duplication of code). We present this class hierarchy in Figure 5.5.


Figure 5.5: Class hierarchy derived from StdExpansion, the base class of the StdRegions Directory.

As is seen in Figure 5.5, the StdRegions hierarchy consists of three levels: the base level from which all StdRegion objects are derived is StdExpansion. This object is then specialized by dimension, yielding StdExpansion0D, StdExpansion1D, StdExpansion2D and StdExpansion3D. The dimension-specific objects are then specialized based upon shape.

The object attributes (variables) at various levels of the hierarchy can be understood in light of Figure 5.6. At its core, an expansion is a means of representing a function over a canonically-defined region of space evaluated at a collection of point positions. The various data members hold information to allow all these basic building blocks to be specified.


Figure 5.6: Diagram to help understand the various data members (object attributes) contained within StdRegions and how they connect with the mathematical representation presented earlier.

The various private, protected and public data members contained within StdRegions are provided in the subsequent sections.

5.2.1 Variables at the Level of StdExpansion

Private: There are private methods but no private data members within StdExpansion.


Public: There are public methods but no public data members within StdExpansion.

5.2.2 Variables at the Level of StdExpansion$D for various Dimensions

Private: There are private methods but no private data members within StdExpansion$D.


Public: There are public methods but no public data members within StdExpansion$D.

5.2.3 Variables at the Level of Shape-Specific StdExpansions




5.2.4 General Layout of the Basis Functions in Memory

5.2.5 General Layout

Basis functions are stored in a 1D array indexed by both mode and quadrature point. The fast index runs over quadrature points while the slow index runs over modes. This was done to match the memory access pattern of the inner product, which is the most frequently computed kernel for most solvers.

Bases are built from the tensor product of three different equation types (subsequently called Type 1, Type II and Type III respectively):

       ||{ 1-2z                  p = 0
ϕ (z) =  1+z                  p = 1
 p     ||( (21-z )(1+z)   1,1
           2     2   Pp-1(z)  2 ≤ p < P

         ( ϕ (z)                   p = 0       0 ≤ q < P
         ||{ ( q1-z)p
ϕpq(z) =   ( -2-) (    )           1 ≤ p < P   q = 0
         ||(   1-z p  1+z- P 2p- 1,1(z) 1 ≤ p < P,  1 ≤ q < P - p
              2      2    q- 1

         |||{  ϕ(qr )p+q                     p = 0      0 ≤ q < P      0 < r < P - q
ϕpqr(z) =     1-2z                         1 ≤ p < P, 0 ≤ q < P - p, r = 0
         |||  (1-z)p+q (1+z)  2p+2q-1,r
         (   -2-      -2-  Pr-1     (z)  1 ≤ p < P, 0 ≤ q < P - p, 1 ≤ r < P - p - qr.

Here, P is the polynomial order of the basis and Ppα,β are the pth order jacobi polynomial.

A Note Concerning Adjustments For C0 Continuity

Before going further it is worth reviewing the spatial shape of each node. The term 1+2z is an increasing function which is equal to zero at z = -1 and equal to one at z = 1. Similarly, 1-2z is a decreasing function which is equal to one at z = -1 and equal to zero at z = +1. These two functions are thus non-zero at one of each of the boundaries. If we need to maintain C0 continuity with surrounding elements (as we do in the continuous Galerkin method), then these local modes must be assembled together with the correct local modes in adjacent elements to create a continuous, global mode. For instance 1+z
 2 in the left element would be continuous with 1-2z in the right element. The union of these two modes under assembly form a single “hat” function. By contrast, functions of the form

  2    2

are zero at both end points z = 1. As a result, they are trivially continuous with any other function which is also equal to zero on the boundary. These “bubble” functions may be treated entirely locally and thus are used to construct the interior modes of a basis. Only bases with p > 1 have interior modes.

All of this holds separately in one dimension. Higher dimensional bases are constructed via the tensor product of 1D basis functions. As a result, we end up with a greater number of possibilities in terms of continuity. When the tensor product is taken between two bubble functions from different bases, the result is still a bubble function. When the tensor product is taken between a hat function and a bubble function we get “edge” modes (in both 2D and 3D). These are non-zero along one edge of the standard domain and zero along the remaining edges. When the tensor product is taken between two hat functions, they form a “vertex” mode which reaches its maximum at one of the vertices and is non-zero on two edges. The 3D bases are constructed similarly.

Based upon this convention, the 1D basis set consists of vertex modes and bubble modes. The 2D basis function set consists of vertex modes, edge modes and bubble modes. The 3D basis set contains vertex modes, edge modes, face modes and bubble modes.

5.2.6 2D Geometries

Quadrilateral Element Memory Organization

Quads have the simplest memory organization. The quadrilateral basis is composed of the tensor product of two Type I functions ϕp0,iq1,j). This function would then be indexed as

1basis0[p*nq0 + i] * basis1[q*nq1 + j]

where nq<b> is the number of quadrature points for the bth basis. Unlike certain mode orderings (e.g. Karniadakis and Sherwin [44]), the two hat functions are accessed as the first and second modes in memory with interior modes placed afterward. Thus,

1basis[i], basis[nq + i]

correspond to 1- z
-2- and 1+z
-2- respectively.

Triangle Element Memory Organization

Due to the use of collapsed coordinates, triangular element bases are formed via the tensor product of one basis function of Type I, and one of Type II, i.e. ϕp0,iϕpq(η1,j)). Since ϕp is also a Type I function, its memory ordering is identical to that used for quads. The second function is complicated by the mixing of ξ0 and ξ1 in the construction of η1.

In particular, this means that the basis function has two modal indices, p and q. While p can run all the way to P, The number of q modes depends on the value of the p index q index such that 0 ≤ q < P - p. Thus, for p = 0, the q index can run all the way up to P. When p = 1, it runs up to P - 1 and so on. Memory is laid out in the same way starting with p = 0. To access all values in order, we write

1mode = 0 
2for p in P: 
3    for q in P - p: 
4        out[mode*nq + q] = basis0[p*nq]*basis1[mode*nq + q] 
5    mode += P-p

Notice the use of the extra “mode” variable. Since the maximum value of q changes with p, basis1 is not simply a linearized matrix and instead has a triangular structure which necessitates keeping track of our current memory location.

The collapsed coordinate system introduces one extra subtlety. The mode


represents the top right vertex in the standard basis. However, when we move to the standard element basis, we are dealing with a triangle which only has three vertices. During the transformation, the top right vertex collapses into the top left vertex. If we naively construct an operators by iterating through all of our modes, the contribution from this vertex to mode Φ01 will not be included. To deal with this, we add its contribution as a correction when computing a kernel. The correction is Φ01 = ϕ0ϕ01 + ϕ1ϕ10 for a triangle.

5.2.7 3D Geometries

Hexahedral Element Memory Organization

The hexahedral element does not differ much from the quadrilateral as it is the simply the product of three Type I functions.

Φpqr = ϕp (ξ0)ϕq(ξ1)ϕr(ξ2).

Prismatic Element Memory Organization

Cross sections of a triangular prism yield either a quad or a triangle based chosen direction. The basis, therefore, looks like a combination of the two different 2D geometries.

Φpqr = ϕp(η0)ϕq(ξ1)ϕpr(η1).

Taking ϕpϕpr on its own produces a triangular face while taking ϕpϕq on its own produces a quadrilateral face. When the three basis functions are combined into a single array (as in the inner product kernel), modes are accessed in the order p,q,r with r being the fastest index and p the slowest. The access pattern for the prism thus looks like

1mode_pqr = 0 
2mode_pr = 0 
3for p in P: 
4    for q in Q: 
5        for r in P - p: 
6            out[mode_pqr*nq + r] = basis0[p*nq]*basis1[q*nq]*basis2[mode_pr + r] 
7        mode_pqr += P - p 
8    mode_pr += P - p

As with the triangle, we have to deal with complications due to collapsed coordinates. This time, the singular vertex from the triangle occurs along an entire edge of the prism. Our correction must be added to a collection of modes indexed by q

Φ0q1+ = ϕ1ϕqϕ10.

Tetrahedral Element Memory Organization

The tetrahedral element is the most complicated of the element constructions. It cannot simply be formed as the composition of multiple triangles since η2 is constructed by mixing three coordinate directions. We thus need to introduce our first Type III function.

Φpqr(η0,η1,η2) = ϕp(η0)ϕpq(η1)ϕpqr(η2).

The r index is constrained by both p and q indices. It runs from P - p - q to 1 in a similar manner to the Type II function. Our typical access pattern is thus

1mode_pqr = 0 
2mode_pq = 0 
3for p in P: 
4    for q in P - p: 
5        for r in P - p - r: 
6            out[mode_pqr*nq + r] = basis0[p*nq]*basis1[mode_pq + q]*basis2[mode_pqr + r] 
7        mode_pqr += (P - p - r) 
8    mode_pq += (P - p)

The tetrahedral element also has to add a correction due to collapsed coordinates. Similar to the prism, the correction must be applied to an entire edge indexed by r

Φ01r+ =  ϕ1ϕ1ϕ11r.

Pyramidic Element Memory Organization

Like the tetrahedral element, a pyramid contains a collapsed coordinate direction which mixes three standard coordinates from the standard region. Unlike the tetrahedra, the collapse only occurs along one axis. Thus it is constructed from two Type I functions and one Type III function

Φpqr = ϕp (η1)ϕq(η2)ϕpqr(η3).

The product ϕpϕq looks like the a quad construction which reflects the quad which serves as the base of the pyramid. A typical memory access looks like

1mode_pqr = 0 
2for p in P: 
3    for q in P - p: 
4        for r in P - p - r: 
5            out[mode_pqr*nq + r] = basis0[p*nq]*basis1[q*nq]*basis2[mode_pqr*nq + r] 
6        mode_pqr += (P - p - r)