# Foundations of Space-Time Finite Element Methods: Polytopes, Interpolation,
and Integration

Cory V. Frontin  Department of Aeronautics and Astronautics, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139  Gage S. Walters
Department of Mechanical Engineering, The Pennsylvania State University,
University Park, Pennsylvania 16802  Freddie D. Witherden  Department of Ocean
Engineering, Texas A&M University, College Station, Texas 77843  Carl W. Lee
Department of Mathematics, University of Kentucky, Lexington, Kentucky 40506
David M. Williams  Department of Mechanical Engineering, The Pennsylvania
State University, University Park, Pennsylvania 16802
<EMAIL_ADDRESS>David L. Darmofal  Department of Aeronautics and
Astronautics, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139

###### Abstract

The main purpose of this article is to facilitate the implementation of space-
time finite element methods in four-dimensional space. In order to develop a
finite element method in this setting, it is necessary to create a numerical
foundation, or equivalently a numerical infrastructure. This foundation should
include a collection of suitable elements (usually hypercubes, simplices, or
closely related polytopes), numerical interpolation procedures (usually
orthonormal polynomial bases), and numerical integration procedures (usually
quadrature rules). It is well known that each of these areas has yet to be
fully explored, and in the present article, we attempt to directly address
this issue. We begin by developing a concrete, sequential procedure for
constructing generic four-dimensional elements (4-polytopes). Thereafter, we
review the key numerical properties of several canonical elements: the
tesseract, tetrahedral prism, and pentatope. Here, we provide explicit
expressions for orthonormal polynomial bases on these elements. Next, we
construct fully symmetric quadrature rules with positive weights that are
capable of exactly integrating high-degree polynomials, e.g. up to degree 17
on the tesseract. Finally, the quadrature rules are successfully tested using
a set of canonical numerical experiments on polynomial and transcendental
functions.

###### keywords:

space-time; finite element methods; quadrature; tesseract; tetrahedral prism;
pentatope

###### MSC:

[2010] 52B11, 65D05, 65D32, 74S05, 76M10

††journal: Applied Numerical Mathematics

##  1 Introduction

Space-time finite element methods are a promising approach for producing
highly accurate and efficient solutions to partial differential equations
(PDEs). More specifically, these methods can adaptively redistribute their
degrees of freedom in both space and time in order to optimally resolve
critical features of the solution. Despite their potential, the application of
these methods to real-world engineering problems remains an open area of
research. While one may simulate space-time problems in 1$d+t$ (two-
dimensional space) or 2$d+t$ (three-dimensional space) by leveraging well-
established finite element techniques, many of these techniques do not carry
over to 3$d+t$ (four-dimensional space). In some sense, four-dimensional space
should be approached as a new frontier, and we must recognize the unique
geometric characteristics which arise in this context.

There are several specific challenges for space-time finite element methods in
the areas of meshing, interpolation, and numerical integration (i.e.
quadrature/cubature). In what follows, we will briefly review these
challenges, along with a short discussion of the associated literature. Our
focus will be limited to state-of-the-art challenges associated with
polynomial-based finite element methods. Note: for a more fundamental
discussion of the methods themselves, including their basic construction, one
may consult the textbooks of Ern and Guermond [1], and Johnson [2], (see
Chapters 6 and 8, respectively). In addition, for a discussion of non-uniform
rational B-spline (NURB)-based methods, please consult the recent work on
isogeometric analysis (IGA) [3, 4, 5, 6, 7, 8, 9, 10, 11].

Our first challenge is to develop unstructured space-time meshes in order to
simulate problems with complex boundaries. In this setting, one may extend the
traditional notion of an unstructured mesh by defining ‘fully-unstructured’
meshes as grids which are unstructured in both space and time, and ‘partially-
unstructured’ meshes as grids which are unstructured in space only. The latter
meshes are usually generated by extruding unstructured spatial meshes in the
temporal direction. Evidently, the partially-unstructured meshes are less
optimal than the fully-unstructured meshes, as all elements are required to
take the same time-step. Nevertheless, fully-unstructured meshes are more
difficult to generate. The problem of mesh generation may become significantly
easier if hybrid meshes are employed, however then it becomes necessary to
select a suitable conforming set of elements (4-polytopes) for populating the
mesh. In three-dimensional space, one may select combinations of hexahedron,
triangular prism, square pyramid, and tetrahedron elements. However, in four-
dimensional space there are significantly more elements to consider (as we
will demonstrate). In order to avoid these complications, previous researchers
have used non-hybrid meshes composed entirely by one of the following standard
elements: tesseracts, tetrahedral prisms, or pentatopes. For example,
partially-unstructured meshes of tesseract elements were used by van der Vegt
and coworkers [12, 13, 14] along with discontinuous Galerkin (DG) methods to
solve laminar compressible flow problems. In addition, Diosady and Murman [15,
16, 17, 18] have used a similar approach to solve turbulent compressible flow
problems. In these studies, the primary objective has been to exploit the
tensor-product structure of the tesseract and maximize computational
efficiency. A contrasting approach is exemplified by Tezduyar and coworkers
[19, 20, 21, 22], who used partially-unstructured meshes of tetrahedral prism
elements along with Galerkin least-squares (GLS) in space/DG in time methods
to solve fluid-structure interaction (FSI) problems. In these studies, the
primary objective was to leverage the flexibility of unstructured tetrahedral
grids for treating complicated geometries. Yang et al. [23] have used a
similar approach along with a hybrid finite-difference/finite-element method
to solve the _fractional_ Bloch-Torrey equations for biomedical flows. (Note:
for detailed discussions of fractional PDEs, including the Fredholm,
Schrödinger, and shallow water equations, one may consult the writings of Al-
Smadi et al. [24, 25, 26, 27], and Karniadakis and coworkers [28, 29]). In an
alternative fashion, Behr [30] has shown that tetrahedral prisms can be
subdivided in order to generate pentatope elements. The resulting meshes have
been used by Behr and coworkers [31, 32, 33] in conjunction with GLS in
space/DG in time methods to solve two-phase incompressible flow problems and
single-phase compressible flow problems. In addition, Wang [34] has used these
meshes in conjunction with space-time DG methods to solve compressible,
viscous flow problems with rotating bodies. Furthermore, Lehrenfeld [35] has
used these meshes in conjunction with enriched space-time DG methods (XFEM-DG
methods) to solve convection-diffusion problems with moving bodies. These
efforts (above) have all focused on partially-unstructured meshes. In order to
move towards fully-unstructured meshes, Mont [36] used a tent-pitching
algorithm [37] to generate unstructured pentatopic meshes that locally satisfy
a causality constraint. In addition, Neumüller and Steinbach [38, 39]
generated isotropic refinement strategies for fully-unstructured meshes of
pentatope elements. Similar work has been subsequently performed by Foteinos
et al. [40] and Belda-Ferrín et al. [41]. Most recently, Caplan and coworkers
[42, 43, 44, 45] have developed highly-anisotropic, fully-unstructured meshes
of pentatope elements. In principle, these meshes possess the greatest
flexibility of any meshes generated to-date for discretizing space-time
problems. Despite this flexibility, there are still opportunities to further
improve the generality of space-time meshes by considering a broader class of
elements beyond the classical tesseract, tetrahedral prism, and pentatope
elements; and furthermore by constructing hybrid meshes (as discussed
earlier).

The next challenge is to develop a family of modal basis functions and/or
nodal basis functions. For example, one may construct a modal basis on the
reference tesseract $[-1,1]^{4}$ by taking tensor products of one-dimensional
Legendre polynomials, and a nodal basis can be constructed in a similar
fashion by using Lagrange polynomials. Unfortunately, this simplicity-of-
construction does not carry over to other element types. For example, the
modal basis on the tetrahedral prism is usually obtained from a tensor product
of the one-dimensional Legendre basis and the three-dimensional PKDO basis
(originally derived in two dimensions by Proriol [46], Koornwinder [47],
Dubiner [48], and Owens [49]). Once this modal basis is obtained, a nodal
basis can be numerically computed using a generalized Vandermonde matrix,
following the procedure on p. 410 of [50] and p. 124 of [51]. However, the
resulting nodal basis functions do not have explicit expressions, and the
conditioning of the Vandermonde matrix can be sensitive to the polynomial
order and the placement of nodal points. In practice, the latter issue is not
especially serious for low-to-moderate polynomial orders ($p\leq 6$) [50].
Note that it is possible to completely avoid this issue by explicitly
constructing nodal basis functions using the Lagrange-mapping procedure of
Hughes ([52] pp. 164-174). However, this procedure is not general, and must be
tailored for each polynomial order. Therefore, the Vandermonde procedure of
[50] is usually preferred. A similar set of issues arise for pentatope basis
functions. In this case, the modal basis is chosen to be the four-dimensional
PKDO basis (defined in [53]), and the nodal basis can be obtained using the
aforementioned Vandermonde procedure. In summary, the construction of modal
and nodal basis functions for tesseract, tetrahedral prism, and pentatope
elements is well-established, with some caution required for nodal basis
functions at high polynomial orders. Nevertheless, there does not appear to be
a well-known resource which aggregates these basis function definitions in one
place. In addition, there does not appear to be any research on the
construction of modal or nodal basis functions for non-standard element types.

The last challenge is to numerically integrate polynomial functions and
transcendental functions on each element. One can easily construct quadrature
rules for the standard tesseract element by taking tensor products of one-
dimensional Gauss-Legendre rules. The resulting rules are ‘good’ in the
following sense: i) they exactly integrate polynomials of degree $\leq p$,
where $p=2M-1$ and $M$ is the number of quadrature points in one dimension;
ii) all quadrature points reside strictly within the tesseract; iii) all
quadrature points are arranged symmetrically; and iv) all quadrature weights
are positive. Unfortunately, while tensor-product rules are ‘good’, they are
not optimal as they possess an excessive number of quadrature points
($N_{p}=M^{4}$ points). Furthermore, they cannot be mapped to tetrahedral
prism or pentatope elements without violating symmetry. As a result, it is
important to generate good, non-tensor-product rules which are optimized for
each element type. In what follows, we briefly review the best known
quadrature rules of this kind on the tesseract and pentatope. Note: to the
authors’ knowledge, no such rules have been constructed for the tetrahedral
prism. On the tesseract, Stroud [54], and Sørevik and Espelid [55] constructed
rules which exactly integrate polynomials of degrees $p=3$ and 5, with 8 and
24 points, respectively. Thereafter, Majorana et al. [56], and Sørevik and
Espelid [57] constructed a rule which exactly integrates a polynomial of
degree $p=7$ with 57 points. Finally, Sørevik [58] discovered a rule which
exactly integrates a polynomial of degree $p=9$ with 160 points. In a similar
fashion, many researchers have constructed quadrature rules for the pentatope
(see the recent review in [59]). In particular, Hammer and Stroud [60, 61]
constructed rules which exactly integrate polynomials of degrees $p=2$ and 3,
with 5 and 15 points, respectively (although the latter rule has points on the
boundary). Thereafter, Gusev et al. [62] have constructed rules which exactly
integrate polynomials of degrees $4,5,6,7,$ and 8 with 20, 30, 56, 76, and 110
points, respectively. In summary, there is a critical absence of good
quadrature rules capable of exactly integrating polynomials with $p\gtrsim 9$
on tesseract and pentatope elements. Furthermore, only tensor-product based
rules are known to exist on the tetrahedral prism.

The purpose of this article is to address some of the key shortcomings
discussed above. In particular, we will: i) identify new types of space-time
elements, and ii) review/expand the interpolation and integration procedures
associated with tesseract, tetrahedral prism, and pentatope elements. The
format of this article is as follows. In section 2, we describe a systematic
procedure for generating sequences of space-time elements. In section 3, we
define the standard tesseract, tetrahedral prism, and pentatope reference
elements, and review their essential properties, including their symmetry
groups and orthonormal basis functions. In section 4, we describe a
computational strategy for constructing quadrature rules on these standard
elements. Thereafter, in sections 5 and 6, we summarize the resulting
quadrature rules, and demonstrate their effectiveness on a suite of numerical
experiments. Finally, in section 7, we summarize the main conclusions of the
article.

##  2 Identification of Element Sequences

In this section, our objective is to construct sequences of space-time
elements. Evidently, there are an infinite number of possible sequences, and
within each sequence there are an unlimited number of polytopes which can be
considered. However, for the sake of convenience, we narrow our focus to
finite sequences of convex polytopes which can be formed by degenerating the
vertices of the tesseract. These polytopes are most amenable to practical
applications, as in principle, a polynomial basis can be formed on each
element by constructing a suitable mapping from the tesseract to the
associated subpolytope. In what follows, we will generate sequences of
elements in this fashion, and introduce formal mathematical notation and
terminology to aid in this process.

The unit $d$-cube ($d$-dimensional cube) is the $d$-polytope $I^{d}$, the
$d$-fold prism of the unit interval $I=[0,1]$. Its vertices are all points of
the form $(x_{1},\ldots,x_{d})$ such that each $x_{i}$ equals either 0 or 1. A
0/1-polytope ${P}$ is the convex hull of a nonempty subset $V$ of the vertices
of the $d$-cube; in other words, the vertex set of ${P}$, ${{\rm vert}}({P})$,
will be denoted by $V$. If $V$ consists of the points $(0,\ldots,0)$,
$(1,0,\ldots,0)$, $(0,1,0,\ldots,0)$, …, $(0,0,\ldots,0,1)$, we will call
${P}$ the standard $d$-simplex. See [63, 64, 65, 66] for a description of some
of the properties of 0/1-polytopes.

In the remainder of this section, we describe the construction of two
sequences of $d$-elements ($d$-dimensional elements). These sequences are
noteworthy, as for $d=2$ and $d=3$ they contain standard elements which are
frequently used to construct 2D and 3D meshes. Naturally, here we are
interested in the generalization of these sequences to 4D, and possibly higher
dimensions.

First, let us consider the construction of the following sequence

${P}_{0},{P}_{1},\ldots,{P}_{2^{d}-2d},\ldots,{P}_{2^{d}-d-1}$

together with a sequence of 0/1-points

${v}_{0},\ldots,{v}_{2^{d}-d-2}$

with the following properties:

  * 1.

The element ${P}_{i}$ is a $d$-dimensional 0/1-polytope,
$i=0,\ldots,2^{d}-d-1$.

  * 2.

The point ${v}_{i}$ is a vertex of ${P}_{i}$, $i=0,\ldots,2^{d}-d-2$.

  * 3.

The element ${P}_{i+1}$ is the convex hull of ${{\rm
vert}}({P}_{i})\setminus{v}_{i}$, $i=0,\ldots,2^{d}-d-2$. Thus ${P}_{i+1}$ is
the result of degenerating the vertex ${v}_{i}$ of ${P}_{i}$.

  * 4.

The first element of the sequence ${P}_{0}$ is the unit $d$-cube.

  * 5.

The intermediate element of the sequence ${P}_{2^{d}-2d}$ is a prism formed by
extruding the $(d-1)$-simplex into one dimension higher. In other words, the
element ${P}_{2^{d}-2d}$ is a prism over the standard $(d-1)$-simplex.

  * 6.

The last element of the sequence ${P}_{2^{d}-d-1}$ is the standard
$d$-simplex.

Note in particular that when $d=4$, then the $3$-simplex is a $3$-dimensional
tetrahedron, and ${P}_{8}$ is the tetrahedral prism obtained by extruding the
tetrahedron in the direction of the fourth coordinate.

The construction of the sequence is recursive on $d=1,2,\ldots$. Naturally,
when $d=1$, ${P}_{0}={P}_{2^{d}-2d}={P}_{2^{d}-d-1}=[0,1]$ and the vertex
sequence is empty. When $d=2$, ${v}_{0}=(1,1)$, ${P}_{0}={P}_{2^{d}-2d}$ (the
unit square), and ${P}_{1}={P}_{2^{d}-d-1}$ (the standard $2$-simplex). For
$d\geq 3$, the vertex sequence ${v}^{\prime}_{i}$ in dimension $d$ is
constructed from the vertex sequence ${v}_{0},\ldots,{v}_{2^{d-1}-(d-1)-2}$ in
dimension $d-1$ by using the following steps.

  1. 1.

For each vertex ${v}_{i}$ create two new vertices
${v}^{\prime}_{2i}=({v}_{i},1)$ and ${v}^{\prime}_{2i+1}=({v}_{i},0)$ by
appending to ${v}_{i}$ first a 1 and second a 0, for
$i=0,\ldots,2^{d-1}-(d-1)-2$. The result is
${v}^{\prime}_{0},\ldots,{v}^{\prime}_{2^{d}-2d-1}$.

  2. 2.

Complete the sequence by setting
${v}^{\prime}_{2^{d}-2d}=(1,0,0,0,\ldots,0,1)$,
${v}^{\prime}_{2^{d}-2d+1}=(0,1,0,0,\ldots,0,1)$,
${v}^{\prime}_{2^{d}-2d+2}=(0,0,1,0,\ldots,0,1)$,
…,${v}^{\prime}_{2^{d}-d-2}=(0,0,0,0,\ldots,1,1)$.

If we regard the $d$-cube as a prism over the $(d-1)$-cube, the above sequence
of degenerations alternates by removing a vertex of the ‘upper’ $(d-1)$-cube,
followed by removing the corresponding vertex of the ‘lower’ $(d-1)$-cube,
thus replicating the $(d-1)$-dimensional degeneration sequence on both the top
and the bottom. At the end of phase 1 above, we will have a prism over the
standard $(d-1)$-simplex. (In fact, ${P}^{\prime}_{2i}$ will be a prism over
${P}_{i}$, $i=0,\ldots,2^{d-1}-d-1$.) From that point on, the sequence in
phase 2 strips away vertices on the top until a single point remains,
resulting in a standard $d$-simplex.

Table 1 shows the vertex sequences for dimensions 2, 3, and 4. (For
convenience we use an abbreviated form of the vertex coordinates.)

$\begin{array}[]{ccc}d=2&d=3&d=4\\\ \hline\cr 11&111&1111\\\ &110&1110\\\
&101&1101\\\ &011&1100\\\ &&1011\\\ &&1010\\\ &&0111\\\ &&0110\\\ &&1001\\\
&&0101\\\ &&0011\end{array}$ Table 1: Sequences of vertices in low dimensions

For $d=3$ the resulting elements $P_{i}$ are displayed in Figure 1.

|  |
---|---|---
${P}_{0}$ | ${P}_{1}$ | ${P}_{2}$
|  |
${P}_{3}$ | ${P}_{4}$ |
Figure 1: Sequence of elements in dimension 3

For $d=4$ the resulting elements $P^{\prime}_{i}$ are displayed in Figure 2.

|  |
---|---|---
${P}^{\prime}_{0}$ | ${P}^{\prime}_{1}$ | ${P}^{\prime}_{2}$
|  |
${P}^{\prime}_{3}$ | ${P}^{\prime}_{4}$ | ${P}^{\prime}_{5}$
|  |
${P}^{\prime}_{6}$ | ${P}^{\prime}_{7}$ | ${P}^{\prime}_{8}$
|  |
${P}^{\prime}_{9}$ | ${P}^{\prime}_{10}$ | ${P}^{\prime}_{11}$
Figure 2: Sequence of elements in dimension 4

Here, the depictions are Schlegel diagrams—central projections of 4-polytopes
onto a 3-dimensional hyperplane with respect to a point placed close to a
facet (3-face) of the 4-polytope; see [67]. In particular, $P^{\prime}_{0}$ is
the tesseract ($4$-cube), $P^{\prime}_{8}$ is a tetrahedral prism (prism over
the standard $3$-simplex), and $P^{\prime}_{11}$ is the pentatope
($4$-simplex). Note that
$P^{\prime}_{0},P^{\prime}_{2},P^{\prime}_{4},P^{\prime}_{6}$ are prisms over,
respectively, the 3-dimensional elements $P_{0},P_{1},P_{2},P_{3}$ of Figure
1.

The second vertex sequence of interest is the one obtained by regarding the
vertices of the $d$-cube as numbers expressed in binary notation. One may list
them in decreasing order, but omit the vertices of the standard $d$-simplex.
The corresponding low dimensional sequences are shown in Table 2.

$\begin{array}[]{ccc}d=2&d=3&d=4\\\ \hline\cr 11&111&1111\\\ &110&1110\\\
&101&1101\\\ &011&1100\\\ &&1011\\\ &&1010\\\ &&1001\\\ &&0111\\\ &&0110\\\
&&0101\\\ &&0011\end{array}$ Table 2: Decreasing binary sequences of vertices
in low dimensions

When $d=3$ this sequence matches the earlier one, but they are different when
$d\geq 4$. In particular, this sequence does not contain the tetrahedral
prism. For $d=4$ the resulting elements $P^{\prime}_{i}$ are displayed in
Figure 3.

|  |
---|---|---
${P}^{\prime}_{0}$ | ${P}^{\prime}_{1}$ | ${P}^{\prime}_{2}$
|  |
${P}^{\prime}_{3}$ | ${P}^{\prime}_{4}$ | ${P}^{\prime}_{5}$
|  |
${P}^{\prime}_{6}$ | ${P}^{\prime}_{7}$ | ${P}^{\prime}_{8}$
|  |
${P}^{\prime}_{9}$ | ${P}^{\prime}_{10}$ | ${P}^{\prime}_{11}$
Figure 3: Second sequence of elements in dimension 4

These two types of sequences are far from exhaustive, since the number of
combinatorially inequivalent 0/1 $d$-polytopes increases doubly-exponentially
with $d$—see [63]. For example, for $d=4$ there are 60,879 such 0/1
$d$-polytopes [68]. Nevertheless, we have highlighted the two sequences above
because for $d=2$ and $d=3$ they contain elements which are frequently used in
conventional finite element methods (as we mentioned previously). Furthermore,
for $d=4$ they contain new element types. Finally, the first sequence also
contains all of the standard 4D element types such as the tesseract, the
tetrahedral prism, and the pentatope. In the next few sections, we will focus
on developing numerical interpolation and integration procedures for these
standard 4D elements.

##  3 Polynomial Bases and Symmetry Groups

###  3.1 Overview

In the previous section, we constructed sequences of elements based on
degenerations of the $4$-cube $\left[0,1\right]^{4}$, and we recovered the
standard tesseract, tetrahedral prism, and pentatope elements. Our choice of
the interval $\left[0,1\right]$ was consistent with standard mathematical
conventions, however, for engineering purposes it is common practice to use
the interval $\left[-1,1\right]$ instead. For certain operations (such as
interpolation and integration) this interval is more convenient, as it is
symmetrical and centered at the origin. Of course, formally speaking, the
mathematical and engineering intervals are topologically equivalent, as they
are related by a linear map. Therefore, the intervals are interchangeable, and
we have a free choice of which one to use.

In what follows, we switch to the engineering convention, and construct a new
set of reference elements for the tesseract, tetrahedral prism, and pentatope.
In addition, we discuss the essential properties of these elements which are
useful for implementation purposes. In particular, integration bounds are
introduced in order to define the domains of our reference elements, basis
functions on each reference element are constructed for interpolation and
quadrature purposes, and lastly we define symmetry groups for each reference
element in order to facilitate the generation of fully symmetric quadrature
rules.

###  3.2 Preliminaries

Throughout this section, whenever we develop basis functions, we construct an
orthonormal basis, i.e.,

$\int_{\Omega}\psi_{ijkq}\left(\boldsymbol{x}\right)\psi_{rstv}\left(\boldsymbol{x}\right)d\boldsymbol{x}=\delta_{ir}\delta_{js}\delta_{kt}\delta_{qv},$

where $\delta_{ir}$ is the Kronecker delta. The simplest orthonormal basis
functions of degree $p$ have the form

$\psi_{ijkq}\left(\boldsymbol{x}\right)=\zeta_{ijkq}\,\hat{P}_{i}^{(\alpha_{1},\beta_{1})}\left(a\right)\hat{P}_{j}^{(\alpha_{2},\beta_{2})}\left(b\right)\hat{P}_{k}^{(\alpha_{3},\beta_{3})}\left(c\right)\hat{P}_{q}^{(\alpha_{4},\beta_{4})}\left(d\right),$

where $0\leq i+j+k+q\leq p$, $a=a\left(x_{1},x_{4}\right)$,
$b=b\left(x_{2},x_{4}\right)$, $c=c\left(x_{3},x_{4}\right)$, and
$d=d\left(x_{4}\right)$ are functions depending on the element type,
$\zeta_{ijkq}$, $\alpha_{1},\ldots,\alpha_{4}$ and
$\beta_{1},\ldots,\beta_{4}$ are constants, and $\hat{P}_{n}^{(\alpha,\beta)}$
is the 1D orthonormal Jacobi polynomial defined as

$\hat{P}_{n}^{(\alpha,\beta)}\left(x\right)=\frac{P_{n}^{(\alpha,\beta)}(x)}{\sqrt{\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}\frac{\left(n+\alpha\right)!\left(n+\beta\right)!}{n!\left(n+\alpha+\beta\right)!}}}.$

Here, $P_{n}^{(\alpha,\beta)}$ is the standard orthogonal (but not
orthonormal) Jacobi polynomial.

###  3.3 Tesseract

Figure 4: Tesseract reference element.

The first reference element is the tesseract with edge length $2$ and volume
$16$ shown in Figure 4. It is straightforward to define bounds of integration
as follows

$\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}dx_{1}dx_{2}dx_{3}dx_{4}=16.$

Trivially, the orthonormal basis inside of the reference tesseract is given by

$\psi_{ijkq}\left(\boldsymbol{x}\right)=\hat{P}_{i}\left(a\right)\hat{P}_{j}\left(b\right)\hat{P}_{k}\left(c\right)\hat{P}_{q}\left(d\right),$

where $0\leq i,j,k,q\leq p$, $a=x_{1}$, $b=x_{2}$, $c=x_{3}$, and $d=x_{4}$.
There are a total of $N_{\text{dof}}=\left(p+1\right)^{4}$ basis functions.

The tesseract has a total of $384$ symmetries. We can apply these symmetries
to a point $(\alpha,\beta,\gamma,\delta)$ with the constraints
$0\leq\alpha,\beta,\gamma,\delta\leq 1$. This will yield a set
$\chi(\alpha,\beta,\gamma,\delta)$ containing its counterparts noting that
$\chi$ denotes the permutations. The cardinality of $\chi$ depends on if any
of the symmetries give rise to identical points. By enumerating this process
with the conditions above we obtain twelve symmetry orbits (groups)

$\begin{aligned} S_{1}&=(0,0,0,0),\\\ S_{2}(\alpha)&=\chi(\alpha,0,0,0),\\\
S_{3}(\alpha)&=\chi(\alpha,\alpha,0,0),\\\
S_{4}(\alpha,\beta)&=\chi(\alpha,\beta,0,0),\\\
S_{5}(\alpha)&=\chi(\alpha,\alpha,\alpha,0),\\\
S_{6}(\alpha,\beta)&=\chi(\alpha,\alpha,\beta,0),\\\
S_{7}(\alpha,\beta,\gamma)&=\chi(\alpha,\beta,\gamma,0),\\\
S_{8}(\alpha)&=\chi(\alpha,\alpha,\alpha,\alpha),\\\
S_{9}(\alpha,\beta)&=\chi(\alpha,\alpha,\alpha,\beta),\\\
S_{10}(\alpha,\beta)&=\chi(\alpha,\alpha,\beta,\beta),\\\
S_{11}(\alpha,\beta,\gamma)&=\chi(\alpha,\alpha,\beta,\gamma),\\\
S_{12}(\alpha,\beta,\gamma,\delta)&=\chi(\alpha,\beta,\gamma,\delta),\end{aligned}\qquad\qquad\begin{aligned}
|S_{1}|&=1,\\\ |S_{2}|&=8,\\\ |S_{3}|&=24,\\\ |S_{4}|&=48,\\\ |S_{5}|&=32,\\\
|S_{6}|&=96,\\\ |S_{7}|&=192,\\\ |S_{8}|&=16,\\\ |S_{9}|&=64,\\\
|S_{10}|&=96,\\\ |S_{11}|&=192,\\\ |S_{12}|&=384.\end{aligned}$

###  3.4 Tetrahedral Prism

Figure 5: Tetrahedral prism reference element.

Our second reference element is the tetrahedral prism which is created by an
extrusion of a standard tetrahedron into four dimensional space; the element
itself is shown in Figure 5. (Note: the standard tetrahedron is the three
dimensional extension of the right triangle.) The volume of the tetrahedral
prism is $\frac{8}{3}$. The integration bounds for the tetrahedral prism are
similar to the standard tetrahedron while accounting for the extrusion in the
fourth dimension,

$\int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{-x_{3}}\int_{-1}^{-1-x_{2}-x_{3}}dx_{1}dx_{2}dx_{3}dx_{4}=\frac{8}{3}.$

Furthermore, the orthonormal polynomial basis for the tetrahedral prism is
found by extending the tetrahedron basis in [53, 50, 69],

$\psi_{ijkq}\left(\boldsymbol{x}\right)=\sqrt{8}\hat{P}_{i}\left(a\right)\hat{P}_{j}^{(2i+1,0)}\left(b\right)\hat{P}_{k}^{(2i+2j+2,0)}\left(c\right)\hat{P}_{q}\left(d\right)\left(1-b\right)^{i}\left(1-c\right)^{i+j},$

where $0\leq i+j+k\leq p$, $0\leq q\leq p$, $a=-2\frac{1+x_{1}}{x_{2}+x_{3}}$
\- 1, $b=2\frac{1+x_{2}}{1-x_{3}}-1$, $c=x_{3}$, $d=x_{4}$, and
$N_{\text{dof}}=\left(p+1\right)^{2}\left(p+2\right)\left(p+3\right)/6$.

For this element, it is convenient to work in barycentric coordinates due to
the two tetrahedral facets. The barycentric coordinates are specified as
follows

$\boldsymbol{\lambda}=\left(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},x_{4}\right)^{T}\quad\text{where}\
0\leq\lambda_{i}\leq
1,\quad\lambda_{1}+\lambda_{2}+\lambda_{3}+\lambda_{4}=1.$

They are related to Cartesian coordinates by

$\boldsymbol{x}=\begin{pmatrix}-1&1&-1&-1&0\\\ -1&-1&1&-1&0\\\ -1&-1&-1&1&0\\\
0&0&0&0&1\end{pmatrix}\boldsymbol{\lambda}.$

The symmetry orbits of the tetrahedral prism are similar to those of the
tetrahedron, with an additional parameter that accounts for the extension of
the prism in the fourth dimension. There are 48 total symmetries for the
tetrahedral prism, and the symmetry orbits are written as follows

$\begin{aligned}
S_{1}&=(\tfrac{1}{4},\tfrac{1}{4},\tfrac{1}{4},\tfrac{1}{4},0),\\\
S_{2}(\delta)&=(\tfrac{1}{4},\tfrac{1}{4},\tfrac{1}{4},\tfrac{1}{4},\pm\delta),\\\
S_{3}(\alpha)&=\chi(\alpha,\alpha,\alpha,1-3\alpha,0),\\\
S_{4}(\alpha,\delta)&=\chi(\alpha,\alpha,\alpha,1-3\alpha,\pm\delta),\\\
S_{5}(\alpha)&=\chi(\alpha,\alpha,\tfrac{1}{2}-\alpha,\tfrac{1}{2}-\alpha,0),\\\
S_{6}(\alpha,\delta)&=\chi(\alpha,\alpha,\tfrac{1}{2}-\alpha,\tfrac{1}{2}-\alpha,\pm\delta),\\\
S_{7}(\alpha,\beta)&=\chi(\alpha,\alpha,\beta,1-2\alpha-\beta,0),\\\
S_{8}(\alpha,\beta,\delta)&=\chi(\alpha,\alpha,\beta,1-2\alpha-\beta,\pm\delta),\\\
S_{9}(\alpha,\beta,\gamma)&=\chi(\alpha,\beta,\gamma,1-\alpha-\beta-\gamma,0),\\\
S_{10}(\alpha,\beta,\gamma,\delta)&=\chi(\alpha,\beta,\gamma,1-\alpha-\beta-\gamma,\pm\delta),\end{aligned}\qquad\qquad\begin{aligned}
|S_{1}|&=1,\\\ |S_{2}|&=2,\\\ |S_{3}|&=4,\\\ |S_{4}|&=8,\\\ |S_{5}|&=6,\\\
|S_{6}|&=12,\\\ |S_{7}|&=12,\\\ |S_{8}|&=24,\\\ |S_{9}|&=24,\\\
|S_{10}|&=48,\end{aligned}$

where $\alpha$, $\beta$, and $\gamma$ are constrained to ensure that
$0\leq\lambda_{i}\leq 1$, $\sum_{i}\lambda_{i}=1$ and $0\leq\delta\leq 1$.

###  3.5 Pentatope

Figure 6: Pentatope reference element.

The final reference element shown in Figure 6 is the pentatope or standard
four-dimensional simplex which has a volume of $\frac{2}{3}$. The bounds of
integration are

$\int_{-1}^{1}\int_{-1}^{-x_{4}}\int_{-1}^{-1-x_{3}-x_{4}}\int_{-1}^{-2-x_{2}-x_{3}-x_{4}}dx_{1}dx_{2}dx_{3}dx_{4}=\frac{2}{3}.$

The orthonormal polynomial basis for the pentatope is obtained from the
generalized simplex formula presented in [53]

$\displaystyle\psi_{ijkq}\left(\boldsymbol{x}\right)$
$\displaystyle=4\hat{P}_{i}\left(a\right)\hat{P}_{j}^{(2i+1,0)}\left(b\right)\hat{P}_{k}^{(2i+2j+2,0)}\left(c\right)\hat{P}_{q}^{(2i+2j+2k+3,0)}\left(d\right)$
$\displaystyle\times\left(1-b\right)^{i}\left(1-c\right)^{i+j}\left(1-d\right)^{i+j+k},$

where $0\leq i+j+k+q\leq p$, $a=-2\frac{x_{1}+1}{x_{2}+x_{3}+x_{4}+1}-1$,
$b=-2\frac{1+x_{2}}{x_{3}+x_{4}}-1$, $c=2\frac{1+x_{3}}{1-x_{4}}-1$,
$d=x_{4}$, and
$N_{\text{dof}}=\left(p+1\right)\left(p+2\right)\left(p+3\right)\left(p+4\right)/24$.
Similar to the tetrahedron, it is convenient to work in barycentric
coordinates. For the pentatope these are

$\boldsymbol{\lambda}=\left(\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4},\lambda_{5}\right)^{T}\quad\text{where}\
0\leq\lambda_{i}\leq
1,\quad\lambda_{1}+\lambda_{2}+\lambda_{3}+\lambda_{4}+\lambda_{5}=1,$

and are related to Cartesian coordinates by

$\boldsymbol{x}=\begin{pmatrix}-1&1&-1&-1&-1\\\ -1&-1&1&-1&-1\\\
-1&-1&-1&1&-1\\\ -1&-1&-1&-1&1\end{pmatrix}\boldsymbol{\lambda}.$

Note that the columns of the matrix correspond to the vertices of the
pentatope in Cartesian coordinates. The pentatope has a total of 120
symmetries as shown in [70], which take the following form

$\begin{aligned}
S_{1}&=(\tfrac{1}{5},\tfrac{1}{5},\tfrac{1}{5},\tfrac{1}{5},\tfrac{1}{5}),\\\
S_{2}(\alpha)&=\chi(\alpha,\alpha,\alpha,\alpha,1-4\alpha),\\\
S_{3}(\alpha)&=\chi(\alpha,\alpha,\alpha,\tfrac{1}{2}-\tfrac{3}{2}\alpha,\tfrac{1}{2}-\tfrac{3}{2}\alpha),\\\
S_{4}(\alpha,\beta)&=\chi(\alpha,\alpha,\alpha,\beta,1-3\alpha-\beta),\\\
S_{5}(\alpha,\beta)&=\chi(\alpha,\alpha,\beta,\beta,1-2\alpha-2\beta),\\\
S_{6}(\alpha,\beta,\gamma)&=\chi(\alpha,\alpha,\beta,\gamma,1-2\alpha-\beta-\gamma),\\\
S_{7}(\alpha,\beta,\gamma,\delta)&=\chi(\alpha,\beta,\gamma,\delta,1-\alpha-\beta-\gamma-\delta),\end{aligned}\qquad\qquad\begin{aligned}
|S_{1}|&=1,\\\ |S_{2}|&=5,\\\ |S_{3}|&=10,\\\ |S_{4}|&=20,\\\ |S_{5}|&=30,\\\
|S_{6}|&=60,\\\ |S_{7}|&=120,\end{aligned}$

where $\alpha$, $\beta$, $\gamma$, and $\delta$ are constrained to ensure that
$0\leq\lambda_{i}\leq 1,\sum_{i}\lambda_{i}=1$.

##  4 Solution Strategy for Quadrature Rules

###  4.1 Our approach

In this section, we describe a strategy for generating symmetric quadrature
rules on space-time elements. Note that this strategy leverages the numerical
infrastructure which was developed in the previous section.

The task of identifying a quadrature rule of strength $\mathcal{P}$ inside of
a domain $\Omega$ can be framed as a non-linear optimization problem where we
seek a set of $N_{p}$ abscissa $\\{\bm{x}_{i}\\}$ and associated weights
$\\{\omega_{i}\\}$ such that

$\int_{\Omega}g_{rstv}(\bm{x})\,d\bm{x}=\sum_{i=1}^{N_{p}}\omega_{i}g_{rstv}(\bm{x}_{i})\quad\text{for
}\quad 0\leq r+s+t+v\leq\mathcal{P},$ (4.1)

where $g_{rstv}(\bm{x})=x_{1}^{r}x_{2}^{s}x_{3}^{t}x_{4}^{v}$ and $r,s,t,v\geq
0$. When we evaluate the left hand side of this equation, the resulting system
is a _non-linear least squares problem_ ; albeit one where we seek a solution
whose residual is zero. Moreover, as formulated in Eq. (4.1), the system is
ill-conditioned and fails to enforce the requirements that the resulting rules
are symmetric, have positive weights, and have points which are contained
within the domain of integration.

The first issue involving ill-conditioning can be addressed by taking

$g_{rstv}(\bm{x})=\psi_{rstv}(\bm{x}),$

where $\psi_{rstv}(\bm{x})$ corresponds to the orthonormal polynomial basis
for $\Omega$ as defined in Section 3. Moreover, on account of the
orthogonality of the basis, it follows that

$\int_{\Omega}\psi_{rstv}(\bm{x})\,d\bm{x}=c\delta_{0r}\delta_{0s}\delta_{0t}\delta_{0v},$

where $c=1/\psi_{0000}$. This enables us to forego having to directly evaluate
integrals inside of our domain in order to formulate our non-linear least
squares problem.

The second issue involving symmetry can be resolved by decomposing the $N_{p}$
points into a set of symmetry groups (orbits). For example, consider taking
$\Omega$ to be the pentatope and $N_{p}$ to be $61$. One possible symmetric
arrangement of these points involves a combination of an $S_{1}$ orbit, an
$S_{3}$ orbit, an $S_{4}$ orbit, and an $S_{5}$ orbit. Summing the number of
points generated by each orbit we find

$1+10+20+30=61,$

as required. By symmetry, it clearly must be the case that all points
belonging to an orbit must share the same quadrature weight. Noting that an
$S_{1}$ orbit has no parameters, $S_{3}$ orbits have one parameter, and
$S_{4}$ and $S_{5}$ orbits both have two parameters, our non-linear least
squares problem has a total of $9$ parameters; $5$ from the orbits and $4$
from the associated quadrature weights. This is a substantial reduction
compared to the $61\times 4+61=305$ parameters in our
original—asymmetric—problem. This greatly reduces the difficulty associated
with obtaining an exact solution to the non-linear least squares problem. The
caveat here is that, in general, we do not know precisely which of the many
distinct orbital decompositions is the correct one. For example, in the case
where $N_{p}=61$ in a pentatope, there are $24$ unique decompositions,
including the one we’ve already presented. Others include a decomposition with
one $S_{1}$ orbit, four $S_{2}$ orbits, one $S_{3}$ orbit, and an $S_{5}$
orbit for a total of $61$ points defined by $14$ parameters, and another
decomposition with one $S_{1}$ orbit and 12 $S_{2}$ orbits for a total of $61$
points defined by $25$ parameters. Given that there is no _a priori_ means of
ascertaining which orbit is most likely to yield a rule, it is generally
necessary to try _all_ of these decompositions. Moreover, we remark here that
the number of decompositions grows rapidly with $N_{p}$. When $N_{p}=600$
there are some $37,457$ decompositions inside of the pentatope.

Symmetry also enables us to reduce the number of basis functions which must be
included in our non-linear least squares problem. Specifically, the ability of
a rule to exactly integrate one particular basis function oftentimes implies
that it is also capable of exactly integrating other basis functions, too. For
example, inside of a tesseract, if a rule is capable of integrating
$\psi_{rstv}(\bm{x})$ then, by symmetry, it must also be able to integrate
$\psi_{srtv}(\bm{x})$ and all other permutations thereof. This enables us to
substantially reduce the size of the non-linear least squares problem.

The third issue involving positivity of the quadrature weights can be
addressed in a variety of ways. One possibility is to compute them using a
_non-negative least squares_ routine. However, the performance of such
routines is typically not competitive with those of ordinary linear least
squares solvers. As such, we have found it more efficient to introduce a
penalty parameter to the non-linear least squares problem. Specifically, we
add an additional equation into the non-linear system which vanishes if all
the weights are positive

$\varepsilon=\sum_{i=1}^{N_{p}}\frac{\omega_{i}-|\omega_{i}|}{2}.$

The fourth issue involving confinement of quadrature points to $\Omega$ can
also be addressed in a variety of ways. One is to employ a non-linear least
squares solver with support for linear equality constraints. However, such
solvers are less widespread than unconstrained non-linear least squares
solvers. An alternative approach (used here), which foregoes the need for
explicit constraints, is to simply _clamp_ or cut off the relevant orbital
parameters before passing them to the $S_{i}$ orbital functions which yield
the Cartesian points. This clamping discourages the solver from attempting to
take orbital parameters outside of their respective ranges.

Another final improvement we can make to our baseline algorithm is to observe
that the non-linear least squares problem is _separable_. If the orbital
parameters—which define the abscissa—are known then our non-linear least
squares problem in Eq. (4.1) reduces down to a _linear_ least squares problem
for the $\omega_{i}$. Thus, it is possible to treat the quadrature weights as
_dependent_ variables. This enables us to further reduce the size of the non-
linear least squares problem; albeit at the cost of having to solve a linear
least squares problem at each iteration.

###  4.2 Alternative approaches

We remark here that in the case of the tesseract that it is possible to
generate symmetric rules through a tensor-production construction of a
suitable one-dimensional rule. If a Gauss–Legendre rule with $l+1$ points is
employed, it follows that the resulting tesseract rule has $(l+1)^{4}$ points
and is capable of exactly integrating $g_{rstv}(\bm{x})$ for $0\leq
r,s,t,v\leq 2l+1$. This property leads to the rules being extremely efficient
for integrating functions with a tensor-product structure. Moreover, when
paired with sum-factorization they can enable computationally efficient
spectral element discretizations [71].

It is also possible to apply these tensor product rules to tetrahedral prisms
and pentatopes by employing the well-known Duffy transformation [72]. In what
follows, we briefly review the Duffy-based procedure for generating quadrature
rules on the pentatope. One is encouraged to consult the appendix of [35] for
more details about the origins of the Duffy transformation.

We begin by introducing
$\bm{y}=\left(y_{1},y_{2},y_{3}\right)\in\mathbb{R}^{3}$ and
$\bm{x}=\left(\bm{y},t\right)\in\mathbb{R}^{4}$. Next, we construct the
following linear map

$\displaystyle\left(\bm{y},t\right)\rightarrow\left(\frac{2\bm{y}+\left(1+t\right)\bm{1}}{1-t},t\right)=\left(\bm{y}^{\prime},t\right),$

where $\bm{1}=\left(1,1,1\right)$. Thereafter, this map can be used to
reformulate a generic integral on the reference pentatope as follows

$\displaystyle\int_{-1}^{1}\int_{-1}^{-t}\int_{-1}^{-1-t-y_{3}}\int_{-1}^{-2-t-y_{3}-y_{2}}f\left(\bm{y},t\right)dy_{1}\,dy_{2}\,dy_{3}\,dt$
$\displaystyle=\frac{1}{8}\int_{-1}^{1}\left(1-t\right)^{3}\int_{-1}^{1}\int_{-1}^{-y_{3}^{\prime}}\int_{-1}^{-1-y_{3}^{\prime}-y_{2}^{\prime}}f\left(\frac{\left(1-t\right)\bm{y}^{\prime}-\left(1+t\right)\bm{1}}{2},t\right)dy_{1}^{\prime}\,dy_{2}^{\prime}\,dy_{3}^{\prime}\,dt$
$\displaystyle=\frac{1}{8}\int_{-1}^{1}\left(1-t\right)^{3}\int_{-1}^{1}\int_{-1}^{-y_{3}^{\prime}}\int_{-1}^{-1-y_{3}^{\prime}-y_{2}^{\prime}}g\left(\bm{y}^{\prime},t\right)d\bm{y}^{\prime}dt$
$\displaystyle=\frac{1}{8}\int_{-1}^{1}\left(1-t\right)^{3}G\left(t\right)dt,$
(4.2)

where on the last two lines we have introduced the following quantities

$\displaystyle g\left(\bm{y}^{\prime},t\right)$
$\displaystyle=f\left(\frac{\left(1-t\right)\bm{y}^{\prime}-\left(1+t\right)\bm{1}}{2},t\right),$
$\displaystyle G\left(t\right)$
$\displaystyle=\int_{-1}^{1}\int_{-1}^{-y_{3}^{\prime}}\int_{-1}^{-1-y_{3}^{\prime}-y_{2}^{\prime}}g\left(\bm{y}^{\prime},t\right)d\bm{y}^{\prime}.$
(4.3)

The key observation is that the integral in Eq. (4.3) is a three-dimensional
integral over the standard tetrahedron. Therefore, one may utilize a generic
three-dimensional integration formula to evaluate Eq. (4.3), in conjunction
with a 1D Gauss-Legendre rule for evaluating Eq. (4.2). In this fashion, the
Duffy transformation replaces a four-dimensional integration procedure with a
combination of three- and one-dimensional procedures. Furthermore, the
transformation can be repeated recursively until only 1D Gauss-Legendre rules
are required. In this case, the Duffy transformation effectively maps a Gauss-
Legendre tensor-product rule from the standard tesseract on to the standard
pentatope.

Somewhat unsurprisingly, none of the rules which arise from the Duffy
transformation are optimal. In addition, although such rules are functional
and easy to generalize, they are asymmetric and suffer from an undesirable
clustering of points.

##  5 Summary of Quadrature Rules

Using the optimisation procedure outlined in Section 4, we obtained a series
of symmetric quadrature rules on the tesseract, tetrahedral prism, and
pentatope. The implementation was carried out by modifying the polyquad code,
an open-source, parallel, highly robust quadrature-generating code written in
C++ (see [69] for details). The resulting quadrature rules, which span orders
2 to 16, are summarised in Table 3. In addition, Table 4 contains an example
of a new rule with strength $\mathcal{P}=9$ and number of points $N_{p}=151$
on the pentatope. All of the rules—which are provided as electronic
supplementary material—have positive weights with none of the points being
outside of the domain. Looking at Table 3, there are several points of note.
Firstly, due to symmetry all of the tesseract rules are capable of exactly
integrating _any_ odd function—including odd-degree polynomials. A consequence
of this is that all of the even rules pick up an extra order of accuracy.
Secondly, we also note that the rules in the tetrahedral prism are limited to
$\mathcal{P}=14$, whereas those inside of the tesseract and the pentatope go
up to $\mathcal{P}=16$. This is due to the fact that the number of unique
topological arrangements of points is an order of magnitude higher in the
tetrahedral prism than in the other elements. As such, it is necessary to
consider a greater number of non-linear least squares problems, which in turn,
translates into an increased computational cost.

| $N_{p}$
---|---
$\mathcal{P}$ | Tesseract | Pentatope | Tetrahedral prism
2 | 16 | 5 | 6
3 | 16 | 15 | 12
4 | 24 | 20 | 20
5 | 24 | 30 | 27
6 | 57 | 56 | 61
7 | 57 | 70 | 72
8 | 160 | 105 | 114
9 | 160 | 151 | 159
10 | 272 | 210 | 259
11 | 272 | 281 | 322
12 | 512 | 445 | 468
13 | 512 | 555 | 608
14 | 728 | 725 | 921
15 | 728 | 905 |
16 | 1384 | 1055 |
Table 3: Number of points $N_{p}$ required for a rule of strength $\mathcal{P}$ inside of each element type. Note that each newly discovered quadrature rule has been underlined. In addition, the explicit parameters (weights and abscissa) for each rule are provided in the electronic supplemental material, in quadruple precision. Orbit | Abscissa | Weight
---|---|---
$S_{1}$ | 0.20000000000000000000000000000000 | 0.026283450664919790544554931007647
$S_{2}$ | 0.033396077194956640848502386475095 | 0.00095166545512843622669550284918655
$S_{2}$ | 0.24389669924344346067467117680514 | 0.010657347396678176509183437341993
$S_{3}$ | 0.0050791921061062008978905902787788 | 0.00031411345650965132958719133177732
$S_{3}$ | 0.2967495659603128481315951580135 | 0.011334029741446587013061345750027
$S_{4}$ | 0.039703540947493234138821823418935 | 0.0034567258020126958766999595741076
| 0.18556651270749874050254409289523 |
$S_{4}$ | 0.098399070565063719667340340181978 | 0.0089501016561121758766577213589523
| 0.25309277814594493439811072010579 |
$S_{4}$ | 0.15177815296591115799210754308869 | 0.0042400257059868015877992812892785
| 0.011737624799814566830560054028264 |
$S_{5}$ | 0.0088561456625601727465721736167059 | 0.0015526370637512027432010914510579
| 0.21639650331182284793855392295174 |
$S_{5}$ | 0.085781201932866919896943973150337 | 0.0028780181522793959965686558635525
| 0.41353118973888920995638429357091 |
Table 4: The $\mathcal{P}=9$ quadrature rule on the pentatope with $N_{p}=151$
points. The remaining rules are provided in the electronic supplemental
material, in quadruple precision.

##  6 Numerical Experiments

In order to validate the quadrature rules developed in Section 5, we will now
perform a set of numerical experiments on the rules of strength 6-12. First,
we will validate that the quadrature rules of degree $\mathcal{P}$ exactly
integrate polynomials of total order $p\leq\mathcal{P}$, to within machine
precision. Thereafter, we will evaluate the ability of the quadrature rules to
integrate transcendental functions reliably, with well-behaved approximation
capabilities. All experiments in the following section are performed using
quad-precision numerical computations to avoid precision issues.

###  6.1 Integration of polynomial functions on solitary elements

Let us consider the integration of polynomials of total order $p$, which can
be given by

$f_{\mathrm{poly}}(\bm{x};p)=\sum_{(r,s,t,v)=0}^{r+s+t+v\leq
p}C_{rstv}g_{rstv}(\bm{x}),$

where $C_{rstv}$ is a constant factor. In this work, $C_{rstv}$ is generated
procedurally by choosing for each $(r,s,t,v)$ a value from a standard normal
distribution (see the detailed discussion in [59]).

For each of the tesseract, tetrahedral prism, and pentatope elements, we have
taken the quadrature rules of order $\mathcal{P}$ and used them to estimate
the integrals of polynomials $f_{\mathrm{poly}}$ of order $p$ as follows

$J=\sum_{i=1}^{N_{p}}\omega_{i}f_{\mathrm{poly}}(\bm{x}_{i}).$

For each monomial term, $g_{rstv}$, we can write the exact integral
$I_{g_{rstv}}$ on an element as follows

$I_{g_{rstv}}=\int_{\kappa}g_{rstv}\left(\bm{x}\right)\,d\bm{x}=\begin{cases}\frac{1}{(r+1)(s+1)(t+1)(v+1)}&\text{$\kappa$
is a tesseract}\\\ \frac{s!t!v!}{(r+1)(s+t+v+3)!}&\text{$\kappa$ is a
tetrahedral prism}\\\ \frac{r!s!t!v!}{(r+s+t+v+4)!}&\text{$\kappa$ is a
pentatope}\end{cases}$

Thus, the exact integral of a polynomial that is being approximated by $J$ can
be given analytically by

$J_{\infty}(\bm{x};p)=\sum_{(r,s,t,v)=0}^{r+s+t+v\leq p}C_{rstv}I_{g_{rstv}}.$

In Figures 7, 8, and 9, the percent error produced by the quadrature rules on
a single element, for each of the element types is shown as a function of the
polynomial order $p$.

Figure 7: Percent error in polynomial integration on a tesseract element due
to quadrature. Figure 8: Percent error in polynomial integration on a
tetrahedral prism element due to quadrature. Figure 9: Percent error in
polynomial integration on a pentatope element due to quadrature.

In these plots, it can be seen that the quadrature rules guarantee exactness
to machine precision for quadruple precision floating point numbers whenever
$\mathcal{P}\geq p$, as we expect.

###  6.2 Integration of transcendental functions on multi-element grids

Next, we seek to quantify the performance of the quadrature rules in
situations that are immediately relevant to finite element applications. We
will consider a domain $\Omega=[0,1]^{4}$, which we will tessellate with three
uniform triangulations of $\Omega$: $\mathcal{T}_{\mathrm{tesseract}}$,
$\mathcal{T}_{\mathrm{prism}}$, and $\mathcal{T}_{\mathrm{pentatope}}$. For
$\mathcal{T}_{\mathrm{tesseract}}$, each dimension of the 4-cube domain is
subdivided $m\in\\{1,2,3,\ldots,12\\}$ times, into $N_{\mathrm{elem}}=m^{4}$
tesseract elements. In the case of the prismatic and pentatope grids, each
tesseract element is further subdivided into either tetrahedral prisms or
pentatopes by using the Kuhn-Freudenthal triangulation on the first three or
all four dimensions of the tesseract element, respectively. This results in
$N_{\mathrm{elem}}=3!m^{4}$ tetrahedral prism or $N_{\mathrm{elem}}=4!m^{4}$
pentatope elements respectively. In each of these three cases, the integral of
a function $f\in L_{1}\left(\Omega\right)$ is then computed by

$J_{\infty}=\int_{\Omega}f(\bm{x})\mathrm{~{}d}\bm{x}\approx
J\equiv\sum_{\kappa\in\mathcal{T}}\sum_{i=1}^{N_{p}}\omega_{i}^{(\kappa)}f\left(\bm{x}_{i}^{(\kappa)}\right),$

where $\bm{x}_{i}^{(\kappa)}$ and $\omega_{i}^{(\kappa)}$ are transformed from
the appropriate reference element to each physical element.

In order to understand the performance of the quadrature rules, we will
evaluate their use on three transcendental functions, namely an exponential
function:

$f_{1}=\exp(x^{2}+2y^{3}+3z^{4}+4t^{5}),$

a sinusoidal function:

$f_{2}=\sin(x^{2}+2y^{3}+3z^{4}+4t^{5}),$

and an even function (also sinusoidal):

$f_{3}=\sin(x^{2}+y^{2}+z^{2}+t^{2}).$

In Figures 10, 11, and 12, the percent error induced by the quadrature
approximation for each function is plotted with respect to the grid
characteristic length

$h=\sqrt[4]{\frac{|\Omega|}{N_{\mathrm{elem}}}},$

with $|\Omega|$ the volume of the domain and $N_{\mathrm{elem}}$ the number of
elements.

Figure 10: Percent error in integration of $f_{1}$ on a tesseract grid due to
quadrature. Figure 11: Percent error in integration of $f_{2}$ on a tesseract
grid due to quadrature. Figure 12: Percent error in integration of $f_{3}$ on
a tesseract grid due to quadrature.

Similarly, in Figures 13, 14, and 15 we can see the results for the
tetrahedral prism elements, and in Figures 16, 17, and 18, the results for the
pentatope elements. The conclusion from these collected results is that
quadrature error seems to be well behaved for each of the included
transcendental functions, and approximations of integrals of similar functions
by these quadrature rules can be expected to be accurate. In each case,
quadrature rules of design order $\mathcal{P}$ produce convergence rates that
are at least $h^{\mathcal{P}}$ and frequently better ($h^{\mathcal{P}+1}$ or
$h^{\mathcal{P}+2}$), which makes them well suited for finite element methods.

We note that, in a few cases lower degree rules outperform higher degree
rules. This is likely due to a coincidental and beneficial alignment between
the symmetry groups of a lower degree rule and the topology of one of our
transcendental functions. Evidently, this trend will not hold for arbitrary
transcendental functions, and the stronger quadrature rule will tend to
outperform the weaker rule in most cases.

Figure 13: Percent error in integration of $f_{1}$ on a prismatic grid due to
quadrature. Figure 14: Percent error in integration of $f_{2}$ on a prismatic
grid due to quadrature. Figure 15: Percent error in integration of $f_{3}$ on
a prismatic grid due to quadrature. Figure 16: Percent error in integration of
$f_{1}$ on a pentatope grid due to quadrature. Figure 17: Percent error in
integration of $f_{2}$ on a pentatope grid due to quadrature. Figure 18:
Percent error in integration of $f_{3}$ on a pentatope grid due to quadrature.

##  7 Conclusion

In this article, we have explored the numerical foundations of finite element
methods for real-world 3$d+t$ problems. We began our discussion by introducing
a rigorous approach for generating four-dimensional elements using sequences
of 0/1-polytopes. Thereafter, we demonstrated that at least one such sequence
can be naturally encoded in terms of four digit binary numbers. Next, we
narrowed our focus to consider the tesseract, tetrahedral prism, and pentatope
elements. For each element type, we presented a brief review of the basic
numerical properties, including explicit definitions of the reference domain,
limits of integration, orthonormal polynomial basis, and symmetry groups. In
addition, for each element, we constructed a new set of fully symmetric
quadrature rules. Finally, we demonstrated the effectiveness of our quadrature
rules by performing a detailed set of numerical experiments.

We anticipate that our work (above) will significantly encourage the
development of _general_ space-time methods for practical applications, and
will not merely be limited to finite element methods. In particular, the newly
developed element sequences will be useful in the construction of hybrid four-
dimensional meshes for space-time visualization and data manipulation
purposes. In addition, the newly developed quadrature rules will be
immediately useful for both space-time finite element _and_ finite volume
methods.

## Declaration of Competing Interests

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgements

The authors would like to thank Dr. Chuluunbaatar Ochbadrakh for helping
error-check the quadrature rules.

## Funding

This research did not receive any specific grant from funding agencies in the
public, commercial, or not-for-profit sectors.

## References

  * [1] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Vol. 159, Springer Science & Business Media, 2013.
  * [2] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Courier Corporation, 2012.
  * [3] U. Langer, S. E. Moore, M. Neumüller, Space–time isogeometric analysis of parabolic evolution problems, Computer Methods in Applied Mechanics and Engineering 306 (2016) 342–363.
  * [4] K. Takizawa, T. E. Tezduyar, T. Terahara, Ram-air parachute structural and fluid mechanics computations with the space–time isogeometric analysis (ST-IGA), Computers & Fluids 141 (2016) 191–200.
  * [5] K. Takizawa, T. E. Tezduyar, Y. Otoguro, T. Terahara, T. Kuraishi, H. Hattori, Turbocharger flow computations with the space–time isogeometric analysis (ST-IGA), Computers & Fluids 142 (2017) 15–20.
  * [6] Y. Otoguro, K. Takizawa, T. E. Tezduyar, Space–time VMS computational flow analysis with isogeometric discretization and a general-purpose NURBS mesh generation method, Computers & Fluids 158 (2017) 189–200.
  * [7] K. Takizawa, T. E. Tezduyar, T. Terahara, T. Sasaki, Heart valve flow computation with the integrated space–time VMS, slip interface, topology change and isogeometric discretization methods, Computers & Fluids 158 (2017) 176–188.
  * [8] K. Takizawa, T. E. Tezduyar, Y. Otoguro, Stabilization and discontinuity-capturing parameters for space–time flow computations with finite element and isogeometric discretizations, Computational Mechanics 62 (5) (2018) 1169–1186.
  * [9] Y. Otoguro, K. Takizawa, T. E. Tezduyar, K. Nagaoka, S. Mei, Turbocharger turbine and exhaust manifold flow computation with the space–time variational multiscale method and isogeometric analysis, Computers & Fluids 179 (2019) 764–776.
  * [10] U. Langer, S. Matculevich, S. Repin, Adaptive space-time isogeometric analysis for parabolic evolution problems, in: Space–Time Methods: Applications to Partial Differential Equations, edited by U. Langer and O. Steinbach, Radon Series on Computational and Applied Mathematics, Vol. 25, De Gruyter, 2019, pp. 141–184.
  * [11] T. Terahara, K. Takizawa, T. E. Tezduyar, A. Tsushima, K. Shiozaki, Ventricle-valve-aorta flow analysis with the Space–Time Isogeometric Discretization and Topology Change, Computational Mechanics (2020) 1–21.
  * [12] J. J. Van der Vegt, H. Van der Ven, Space–time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows: I. General formulation, Journal of Computational Physics 182 (2) (2002) 546–585.
  * [13] H. Van der Ven, J. J. van der Vegt, Space–time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows: II. Efficient flux quadrature, Computer Methods in Applied Mechanics and Engineering 191 (41-42) (2002) 4747–4780.
  * [14] C. M. Klaij, J. J. van der Vegt, H. van der Ven, Space–time discontinuous Galerkin method for the compressible Navier–Stokes equations, Journal of Computational Physics 217 (2) (2006) 589–611.
  * [15] L. T. Diosady, S. M. Murman, Higher-order methods for compressible turbulent flows using entropy variables, in: 53rd AIAA Aerospace Sciences Meeting, 2015\.
  * [16] L. T. Diosady, S. M. Murman, Tensor-product preconditioners for higher-order space–time discontinuous Galerkin methods, Journal of Computational Physics 330 (2017) 296–318.
  * [17] L. T. Diosady, S. M. Murman, A linear-elasticity solver for higher-order space-time mesh deformation, in: 2018 AIAA Aerospace Sciences Meeting, 2018.
  * [18] L. T. Diosady, S. M. Murman, Scalable tensor-product preconditioners for high-order finite-element methods: Scalar equations, Journal of Computational Physics 394 (2019) 759–776.
  * [19] T. E. Tezduyar, S. Sathe, R. Keedy, K. Stein, Space–time finite element techniques for computation of fluid–structure interactions, Computer Methods in Applied Mechanics and Engineering 195 (17-18) (2006) 2002–2027.
  * [20] T. E. Tezduyar, S. Sathe, Modelling of fluid–structure interactions with the space–time finite elements: Solution techniques, International Journal for Numerical Methods in Fluids 54 (6-8) (2007) 855–900.
  * [21] T. E. Tezduyar, K. Takizawa, C. Moorman, S. Wright, J. Christopher, Space–time finite element computation of complex fluid–structure interactions, International Journal for Numerical Methods in Fluids 64 (10-12) (2010) 1201–1218.
  * [22] T. E. Tezduyar, K. Takizawa, Space–time computations in practical engineering applications: a summary of the 25-year history, Computational Mechanics 63 (4) (2019) 747–753.
  * [23] Z. Yang, F. Liu, Y. Nie, I. Turner, An unstructured mesh finite difference/finite element method for the three-dimensional time-space fractional Bloch-Torrey equations on irregular domains, Journal of Computational Physics 408 (2020) 109284.
  * [24] M. Al-Smadi, A. Freihat, H. Khalil, S. Momani, R. Ali Khan, Numerical multistep approach for solving fractional partial differential equations, International Journal of Computational Methods 14 (03) (2017) 1750029.
  * [25] M. Al-Smadi, O. A. Arqub, Computational algorithm for solving fredholm time-fractional partial integrodifferential equations of dirichlet functions type with error estimates, Applied Mathematics and Computation 342 (2019) 280–294.
  * [26] M. Al-Smadi, O. A. Arqub, S. Hadid, An attractive analytical technique for coupled system of fractional partial differential equations in shallow water waves with conformable derivative, Communications in Theoretical Physics 72 (8) (2020) 085001.
  * [27] M. Al-Smadi, O. A. Arqub, S. Momani, Numerical computations of coupled fractional resonant Schrödinger equations arising in quantum mechanics under conformable fractional derivative sense, Physica Scripta 95 (7) (2020) 075218\.
  * [28] G. E. Karniadakis, J. S. Hesthaven, I. Podlubny, Special Issue on “Fractional PDEs: Theory, Numerics, and Applications”, Journal of Computational Physics 293 (2015) 1–3.
  * [29] E. Kharazmi, Z. Mao, G. Pang, M. Zayernouri, G. E. Karniadakis, Fractional calculus and numerical methods for fractional PDEs, in: First Congress of Greek Mathematicians: Proceedings of the Congress held in Athens, Greece, June 25–30, 2018, Walter de Gruyter, 2020, p. 91.
  * [30] M. Behr, Simplex space–time meshes in finite element simulations, International Journal for Numerical Methods in Fluids 57 (9) (2008) 1421–1434.
  * [31] V. Karyofylli, M. Frings, S. Elgeti, M. Behr, Simplex space-time meshes in two-phase flow simulations, International Journal for Numerical Methods in Fluids 86 (3) (2018) 218–230.
  * [32] V. Karyofylli, L. Wendling, M. Make, N. Hosters, M. Behr, Simplex space-time meshes in thermally coupled two-phase flow simulations of mold filling, Computers & Fluids 192 (2019) 104261.
  * [33] M. von Danwitz, V. Karyofylli, N. Hosters, M. Behr, Simplex space-time meshes in compressible flow simulations, International Journal for Numerical Methods in Fluids 91 (1) (2019) 29–48.
  * [34] L. Wang, Discontinuous galerkin methods on moving domains with large deformations, Ph.D. thesis, University of California, Berkeley (2015).
  * [35] C. Lehrenfeld, The Nitsche XFEM-DG space-time method and its implementation in three space dimensions, SIAM Journal on Scientific Computing 37 (1) (2015) A245–A270.
  * [36] A. D. Mont, Adaptive unstructured spacetime meshing for four-dimensional spacetime discontinuous galerkin finite element methods, Master’s thesis, University of Illinois at Urbana-Champaign (2012).
  * [37] A. Üngör, A. Sheffer, Tent-Pitcher: A meshing algorithm for space-time discontinuous Galerkin methods, in: In proc. 9th int’l. meshing roundtable, Citeseer, 2000, pp. 111–122.
  * [38] M. Neumüller, O. Steinbach, Refinement of flexible space–time finite element meshes and discontinuous Galerkin methods, Computing and Visualization in Science 14 (5) (2011) 189–205.
  * [39] M. Neumüller, Space-time methods: Fast solvers and applications, Ph.D. thesis, Graz University of Technology (2013).
  * [40] P. Foteinos, N. Chrisochoides, 4d space–time Delaunay meshing for medical images, Engineering with Computers 31 (3) (2015) 499–511.
  * [41] G. Belda-Ferrín, A. Gargallo-Peiró, X. Roca, Local bisection for conformal refinement of unstructured 4D simplicial meshes, in: International Meshing Roundtable, Springer, 2018, pp. 229–247.
  * [42] P. C. Caplan, R. Haimes, D. L. Darmofal, M. C. Galbraith, Anisotropic geometry-conforming d-simplicial meshing via isometric embeddings, Procedia Engineering 203 (2017) 141–153.
  * [43] P. Caplan, R. Haimes, D. Darmofal, M. Galbraith, Extension of local cavity operators to 3d+t space-time mesh adaptation, in: AIAA Scitech 2019 Forum, 2019\.
  * [44] P. C. Caplan, Four-dimensional anisotropic mesh adaptation for spacetime numerical simulations, Ph.D. thesis, Massachusetts Institute of Technology (2019).
  * [45] P. C. Caplan, R. Haimes, D. L. Darmofal, Four-dimensional anisotropic mesh adaptation, Computer-Aided Design (2020) 102915.
  * [46] J. Proriol, Sur une famille de polynomes á deux variables orthogonaux dans un triangle, Comptes Rendus Hebdomadaires des Seances de l Academie des Sciences 245 (26) (1957) 2459–2461.
  * [47] T. Koornwinder, Two-variable analogues of the classical orthogonal polynomials, in: Theory and Application of Special Functions, Elsevier, 1975, pp. 435–495.
  * [48] M. Dubiner, Spectral methods on triangles and other domains, Journal of Scientific Computing 6 (4) (1991) 345–390.
  * [49] R. G. Owens, Spectral approximations on the triangle, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 454 (1971) (1998) 857–872.
  * [50] J. S. Hesthaven, T. Warburton, Nodal discontinuous Galerkin methods: Algorithms, analysis, and applications, Springer Science & Business Media, 2007\.
  * [51] G. Karniadakis, S. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics, Oxford University Press, 2013.
  * [52] T. J. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Courier Corporation, 2012.
  * [53] T. Warburton, J. S. Hesthaven, On the constants in hp-finite element trace inverse inequalities, Computer Methods in Applied Mechanics and Engineering 192 (25) (2003) 2765–2773.
  * [54] A. H. Stroud, Approximate Calculation of Multiple Integrals, Prentice-Hall, 1971\.
  * [55] T. Sørevik, T. O. Espelid, Fully symmetric integration rules for the 4-cube, Tech. rep., University of Bergen, Norway (1987).
  * [56] C. Majorana, S. Odorizzi, R. Vitaliani, Shortened quadrature rules for finite elements, Advances in Engineering Software 4 (2) (1982) 52–57.
  * [57] T. Sørevik, T. O. Espelid, Fully symmetric integration rules for the 4-cube, BIT Numerical Mathematics 29 (1) (1989) 148–153.
  * [58] T. Sørevik, Reliable and efficient adaptive quadrature algorithms, Ph.D. thesis, University of Bergen, Norway (1988).
  * [59] D. M. Williams, C. V. Frontin, E. A. Miller, D. L. Darmofal, A family of symmetric, optimized quadrature rules for pentatopes, Computers & Mathematics with Applications 80 (5) (2020) 1405–1420.
  * [60] P. C. Hammer, A. H. Stroud, Numerical integration over simplexes, Mathematical Tables and Other Aids to Computation 10 (55) (1956) 137–139.
  * [61] A. Stroud, Some approximate integration formulas of degree 3 for an n-dimensional simplex, Numerische Mathematik 9 (1) (1966) 38–45.
  * [62] A. A. Gusev, V. P. Gerdt, O. Chuluunbaatar, G. Chuluunbaatar, S. I. Vinitsky, V. Derbov, A. Góźdź, P. Krassovitskiy, Symbolic-numerical algorithms for solving elliptic boundary-value problems using multivariate simplex lagrange elements, in: International Workshop on Computer Algebra in Scientific Computing, Springer, 2018, pp. 197–213.
  * [63] G. M. Ziegler, Lectures on $0/1$-polytopes, in: Polytopes—Combinatorics and Computation (Oberwolfach, 1997), Vol. 29 of DMV Sem., Birkhäuser, Basel, 2000, pp. 1–41.
  * [64] O. Aichholzer, Extremal properties of 0/1-polytopes of dimension 5, in: Polytopes—Combinatorics and Computation (Oberwolfach, 1997), Vol. 29 of DMV Sem., Birkhäuser, Basel, 2000, pp. 111–130.
  * [65] B. Baumeister, C. Haase, B. Nill, A. Paffenholz, On permutation polytopes, Advances in Mathematics 222 (2009) 431–452.
  * [66] A. Cihangir, On 0/1-polytopes with nonobtuse triangulations, Applied Mathematics and Computation 267 (2015) 17–27.
  * [67] B. Grünbaum, Convex polytopes, 2nd Edition, Vol. 221 of Graduate Texts in Mathematics, Springer-Verlag, New York, 2003.
  * [68] C. Zong, What is known about unit cubes?, Bulletin of the American Mathematical Society 42 (2) (2005) 181–211.
  * [69] F. D. Witherden, P. E. Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers & Mathematics with Applications 69 (10) (2015) 1232–1241.
  * [70] J. Maeztu, E. Sáinz de La Maza, Consistent structures of invariant quadrature rules for the $n$-simplex, Mathematics of Computation 64 (211) (1995) 1171–1192.
  * [71] S. A. Orszag, Spectral methods for problems in complex geometrics, in: Numerical Methods for Partial Differential Equations, Elsevier, 1979, pp. 273–305.
  * [72] M. G. Duffy, Quadrature over a pyramid or cube of integrands with a singularity at a vertex, SIAM Journal on Numerical Analysis 19 (6) (1982) 1260–1262.