∎

11institutetext: M. Petrov 22institutetext: Department of Technical Mechanics,
Technical University, 5300 Gabrovo, Bulgaria 33institutetext: T. Todorov
44institutetext: Department of Mathematics, Informatics and Natural Sciences,
Technical University, 5300 Gabrovo, Bulgaria 55institutetext: G. Walters
66institutetext: Department of Mechanical Engineering, The Pennsylvania State
University, University Park, Pennsylvania 16802, USA 77institutetext: D.
Williams 88institutetext: Department of Mechanical Engineering, The
Pennsylvania State University, University Park, Pennsylvania 16802, USA
88email<EMAIL_ADDRESS>99institutetext: F. Witherden
1010institutetext: Department of Ocean Engineering, Texas A&M University,
College Station, Texas 77843, USA

# Enabling four-dimensional conformal hybrid meshing with cubic pyramids

Miroslav S. Petrov  Todor D. Todorov  Gage S. Walters  David M. Williams
Freddie D. Witherden

(Received: date / Accepted: date)

###### Abstract

The main purpose of this article is to develop a novel refinement strategy for
four-dimensional hybrid meshes based on cubic pyramids. This optimal
refinement strategy subdivides a given cubic pyramid into a conforming set of
congruent cubic pyramids and invariant bipentatopes. The theoretical
properties of the refinement strategy are rigorously analyzed and evaluated.
In addition, a new class of fully symmetric quadrature rules with positive
weights are generated for the cubic pyramid. These rules are capable of
exactly integrating polynomials with degrees up to 12. Their effectiveness is
successfully demonstrated on polynomial and transcendental functions. Broadly
speaking, the refinement strategy and quadrature rules in this paper open new
avenues for four-dimensional hybrid meshing, and space-time finite element
methods.

###### Keywords:

4D hybrid meshes cubic pyramids bipentatopes pentatopes optimal refinement
strategy quadrature

###### MSC:

52B11 65N30 65N50 65D30 65D32

††journal: Numerical Algorithms

##  1 Introduction

The classical way of solving nonstationary boundary value problems applies
separate discretizations in space and time. One promising alternative to this
approach is the finite element method with full space-time discretization. In
principle, this method can be very effective for solving hyperbolic and
parabolic problems. Some basic examples of space-time methods can be found in
steinbach2015space ; steinbach2017algebraic ; toulopoulos2020space  and the
recent survey of Langer and Steinbach langer2019space . In addition, a
representative example of this work is provided by Dumont et al.
dumont2012space  who developed a space-time finite element method to solve
elastodynamics problems. Since the dynamics problem in dumont2012space
involves a three-dimensional spatial body, the authors have constructed four-
dimensional meshes. Dumont et al. have not developed general, unstructured
triangulations in four-dimensional domains. Rather, they firstly make a
spatial triangulation and then extend the domain in the temporal direction to
obtain a four-dimensional finite element triangulation. A similar approach for
generating four-dimensional simplex meshes has been employed by Neumüller and
Karabelas in neumuller2019generating . A detailed survey of four-dimensional
structured and unstructured simplicial meshing techniques is provided by
Caplan caplan2019four  and Frontin et al. frontin2020polytopes . Finally,
examples of refinement strategies for high-dimensional simplicial meshes have
been developed in the work of Brandts et al. brandts2007simplicial , and
Korotov and Křížek korotov2014red .

There have been virtually no efforts to develop four-dimensional hybrid
meshes. However, in the last few decades, there has been interest in the
development of hybrid _three-dimensional_ meshes bergot2010higher ; cao2012new
, with some important work being contributed by Yamakawa and coworkers
yamakawa2011subdivision ; yamakawa2003increasing ; yamakawa2009converting ;
yamakawa201088 ; yamakawa2011automatic . In yamakawa2011automatic , Yamakawa
and Shimada developed a procedure for creating hexahedron dominant meshes in
order to model thin-wall structures. Additionally, Yamakawa et al. described a
conforming coupling between hexahedron and tetrahedron meshes by using square
pyramid interface elements yamakawa2009converting , and using triangular prism
interface elements yamakawa2011subdivision . The conforming coupling of
hexahedra and tetrahedra via pyramids has been discussed independently by
Meshkat and Talmor meshkat2000generating , and Devloo et al. devloo2019high .
In addition, Pantano and Averill pantano2007penalty  have used a non-
conforming, penalty-based procedure to couple independently modeled three-
dimensional finite element meshes. Similar non-conforming techniques using
mortar elements have been explored by Maday et al. maday1988nonconforming  and
Seshaiyer and Suri seshaiyer2000hp . Furthermore, hybrid non-conforming
hexahedral-tetrahedral meshes have been examined by Reberol and Lévy
reberol2016low . Unfortunately, these non-conforming meshes are incompatible
with most standard finite element methods, with the exception of discontinuous
Galerkin methods. Evidently, the most convenient way to unify independently
created finite element meshes is through the conforming coupling of the
elements in the interface subdomain. This strongly motivates the importance of
transitional elements, such as pyramidal and prismatic elements.

Let us briefly review previous efforts to analyze transitional elements.
Evidently, the three-dimensional prismatic elements are relatively
straightforward to analyze, as they can be constructed by forming tensor
products between line segments and triangles, whereas, a similar procedure is
impossible for pyramids. As a result, prismatic elements inherit their
interpolation and quadrature procedures from well-established 1D/2D
procedures, and researchers have primarily focused their attention on
pyramids. The pioneering work on pyramids was performed by Bedrosian
bedrosian1992shape  in the early nineties. He was the first researcher who
explained the role of pyramidal elements for coupling different kinds of
meshes, and developed the first interpolation and integration procedures. A
number of researchers have followed in Bedrosian’s footsteps
ainsworth2017lowest ; bergot2010higher ; chan2015comparison ; chan2015hp ;
chan2016orthogonal ; chan2016short ; coulomb1997pyramidal ;
gillette2016serendipity . In particular, Bergot et al. bergot2010higher
performed a deep analysis of the interpolation properties of pyramidal
elements. Chan and Warburton extended this work, developing alternative sets
of basis functions chan2016orthogonal ; chan2016short , analyzing the
numerical stability of interpolation points chan2015comparison , and deriving
trace inequalities chan2015hp . Furthermore, Gillette gillette2016serendipity
recently used techniques from finite element exterior calculus
arnold2010finite ; arnold2006finite  to construct seredipity-type basis
functions on pyramids. Quadrature formulas on pyramidal elements have been
created and analyzed in bergot2010higher ; nigam2012numerical ;
witherden2015identification . The error estimates arising from the quadrature
rules of Bergot et al. bergot2010higher  have been improved by Nigam and
Phillips nigam2012numerical  in order to obtain the classical rate of
convergence for second-order boundary value problems. However, the best set of
rules to date (in terms of integration strength for a minimum number of
points) was identified by Witherden and Vincent in witherden2015identification
.

Now, let us return our attention to four-dimensional space. It turns out that
the three-dimensional definitions of transitional elements do not extend to
higher dimensions. For example, while the triangular prism provides a
conforming interface between the hexahedron and the tetrahedron in three
dimensions, its analog in four dimensions (the tetrahedral prism) does _not_
provide a conforming interface between the tesseract and the pentatope. This
immediately follows from the fact that the tetrahedral prism has two
tetrahedral facets and four triangular prismatic facets, and therefore it
cannot conformally interface with the tesseract. In a similar fashion, the
cubic pyramid (the four-dimensional analog of the square pyramid) has one
hexahedral facet and six square pyramidal facets, and therefore it cannot
conformally interface with the pentatope. Fortunately, the lack of
transitional elements between tesseracts and pentatopes does not completely
prevent the formulation of hybrid meshes in four dimensions. We will show that
it is still possible to build hybrid meshes using tesseract, cubic pyramid,
and pentatopal elements. We note that some basic geometric properties of these
elements are discussed in the work of Coxeter coxeter1940regular ;
coxeter1973regular , Sommerville mclaren1958introduction , and Zamboj
zamboj2018sections .

The main goals of this paper are: i) to develop conforming hybrid four-
dimensional meshes, ii) to create an optimal refinement strategy for these
meshes, and iii) to develop numerical integration procedures on the elements
of these meshes. The major contributions of the paper are summarized as
follows. First, we identify a four-dimensional conforming coupling between
tesseract elements and cubic pyramid elements. Next, the tesseract and cubic
pyramid elements are subdivided, while maintaining conformity. Evidently,
while the tesseracts can be uniformly subdivided into smaller tesseracts, the
cubic pyramids cannot be subdivided into smaller cubic pyramids without
leaving gaps. Therefore, in the refinement strategy for the cubic pyramids, we
use a combination of congruent cubic pyramids _and_ invariant bipentatopes. A
simple two-level refinement tree is obtained. The theoretical properties of
the refinement strategy are thoroughly analyzed. The proposed theoretical
results cannot be improved, otherwise the cubic pyramid would be invariant
with a one-level refinement tree, which is obviously not the case. Finally, we
conclude our work by developing numerical integration procedures for the cubic
pyramid elements. We note that, while such procedures are already well-
established for the tesseract and pentatope elements frontin2020polytopes ;
williams2020family , they have yet to be developed for cubic pyramid elements.

The format of this paper is as follows. In Section 2, we introduce some
standard notation and terminology. In Section 3, we formulate new hybrid
meshes of tesseract, cubic pyramid, and bipentatope elements, along with a
non-degenerate mesh refinement strategy. In Section 4, we develop theoretical
results which govern the hybrid meshes. In Section 5, we introduce a new set
of fully symmetric quadrature rules for cubic pyramid elements. Finally, in
Section 6, we summarize the main conclusions of the paper.

##  2 Preliminary concepts

###  2.1 Background and motivation

The four-dimensional space $\mathbb{R}^{4}$ has unique geometric properties
which cannot be found in other multidimensional Euclidean spaces excluding the
two-dimensional one. This fact can be illustrated by comparing the three- and
four-dimensional spaces. It is well-known that the 3-cube cannot be
triangulated by standard tetrahedra (i.e. cube corners) korotov2014red .
Rather, the cube is usually triangulated into five simplices, four of which
are standard tetrahedra and one of which is a regular tetrahedron.
Additionally, it is impossible to triangulate the standard tetrahedron with
standard tetrahedra. Instead, there exists a partition of the standard
tetrahedron into six standard tetrahedra and one regular tetrahedron
todorov2013optimal . Interestingly enough, these principles do not extend to
the four-dimensional case. Petrov and Todorov petrov2018stable  have proved
that each tesseract can be divided into standard pentatopes (i.e. tesseract
corners), and each standard pentatope can be partitioned into standard
pentatopes. The latter fact means that an arbitrary 4D _canonical_ domain can
be divided into standard pentatopes petrov2018stable , which are the most
convenient elements from a computational point of view. The authors have
tested six- and eight-dimensional hypercubes for such properties, but without
success. Therefore, given its useful properties, more attempts to explore and
analyze four-dimensional space are certainly warranted.

In what follows, we will introduce some useful definitions which facilitate
our subsequent discussions of $\mathbb{R}^{4}$.

###  2.2 Definitions

###### Definition 1

A simply connected domain in $\mathbb{R}^{4}$ is called canonical if there
exists a conforming partition of the domain into tesseracts (4-cubes).

###### Definition 2

Each $n$-dimensional hypercube can be divided into $2n$ $(n-1)$-hypercube
pyramids. These pyramids are called canonical pyramids.

###### Definition 3

Let $\Omega$ be a nondegenerated polytope in $\mathbb{R}^{4}$, and $T_{i}$,
$i=1,2,\ldots,k,$ be $4$-dimensional finite elements. The triangulation

$\displaystyle\tau=\left\\{T_{i}\subset\mathbb{R}^{4}\ |\
\Omega=\bigcup_{i=1}^{k}T_{i}\right\\},$

of the polytope $\Omega$ is said to be consistent if $T_{i}$ and $T_{j}$,
$1\leq i<j\leq k$, share nothing other than the empty set, a vertex, an edge
or an $\ell$-dimensional facet, $\ell=2,3$.

###### Definition 4

The polytope

$\displaystyle T_{-i}=[t_{1},t_{2},\dots,t_{i-1},t_{i+1},\dots,t_{n+1}],\quad
i=2,\dots,n,$

related to the vertex $t_{i}$ is obtained from the polytope
$T=[t_{1},t_{2},\dots,t_{n+1}]$, $n\in\mathbb{N}$ by removing the vertex
$t_{i}$. The polytopes $T_{-1}$ and $T_{-(n+1)}$ are defined in a similar
fashion as

$\displaystyle T_{-1}=[t_{2},t_{3},\dots,t_{n+1}],\quad
T_{-(n+1)}=[t_{1},t_{2},\dots,t_{n}].$

###### Definition 5

The polytopes $T_{1}$ and $T_{2}$ are congruent if one of them can be obtained
from the other by applying a linear transformation

$\displaystyle T_{2}=\underline{b}+cQT_{1},$

where $c\not=0$ is a scaling factor, $\underline{b}$ is a translation vector,
and $Q$ is an orthogonal matrix.

###### Definition 6

We say that two polytopes $T_{1}$ and $T_{2}$ are from the same class if they
are congruent.

The class

$\displaystyle[K]=\\{T\subset\mathbb{R}^{4}\ |\ T\cong K\\},$

consists of all equivalent pyramids to the pyramid $K$ in $\mathbb{R}^{4}$
with respect to the congruence relation.

###### Definition 7

A polytope $T$ is said to be invariant concerning a refinement strategy
$\mathcal{A}$ if all elements of $\mathcal{A}T$ belong to $[T]$.

###### Definition 8

The degeneracy measure of an arbitrary pyramid $K=[k_{1},k_{2},$ $\ldots,$
$k_{n+1}]$ is equal to

$\delta(K)=\frac{h(K)\cdot{\rm vol}(\partial K)}{8\cdot{\rm vol}(K)},$

where $h(K)$ is the diameter of $K$.

Here, we use ${\rm vol}(T)$ and ${\rm vol}(\partial T)$ instead of ${\rm
vol}_{4}(T)$ and ${\rm vol}_{3}(\partial T)$ in order to avoid complicated
notation.

###  2.3 Reference elements

In this section, we introduce a set of convenient reference elements in
$\mathbb{R}^{4}$. We review the important properties of these elements,
including their vertex locations, integration limits, and orthonormal
polynomial bases.

In order to prepare our discussion of the reference elements, let us briefly
review the definition of an orthonormal basis. Broadly speaking, the basis
maintains the following key property

$\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 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)f\left(\boldsymbol{x}\right),$

where $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)$,
$d=d\left(x_{4}\right)$, and $f=f\left(\boldsymbol{x}\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)}$ are the 1D orthonormal Jacobi polynomials
defined as

$\hat{P}_{n}^{(\alpha,\beta)}\left(x_{1}\right)=\frac{P_{n}^{(\alpha,\beta)}(x_{1})}{\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, the functions $P_{n}^{(\alpha,\beta)}$ are the well-known orthogonal
Jacobi polynomials, which themselves are _not_ orthonormal. The scale factor
under the square root operator provides the desired normalization.

It is convenient to omit the Jacobi polynomial superscripts when
$\alpha=\beta=0$. Consequently, we frequently write $\hat{P}_{n}$ in place of
$\hat{P}_{n}^{\left(0,0\right)}$.

Now, let us turn our attention to the definitions of the reference elements.

####  2.3.1 Tesseract

Consider the reference tesseract $T^{\ast}$ centered at the origin, having
edge length $l[T^{\ast}]=2$. The vertices are defined such that

$\displaystyle T^{\ast}=\Big{[}$ $\displaystyle
t_{1}\left(-1,-1,-1,-1\right),\;t_{2}\left(-1,1,-1,-1\right),\;t_{3}\left(-1,-1,-1,1\right),\;t_{4}\left(1,-1,-1,-1\right),$
$\displaystyle
t_{5}\left(-1,1,-1,1\right),\;t_{6}\left(1,-1,-1,1\right),\;t_{7}\left(1,1,-1,-1\right),\;t_{8}\left(1,1,-1,1\right),$
$\displaystyle
t_{9}\left(-1,-1,1,1\right),\;t_{10}\left(-1,1,1,1\right),\;t_{11}\left(-1,-1,1,-1\right),\;t_{12}\left(1,-1,1,1\right),$
$\displaystyle
t_{13}\left(-1,1,1,-1\right),\;t_{14}\left(1,-1,1,-1\right),\;t_{15}\left(1,1,1,1\right),\;t_{16}\left(1,1,1,-1\right)\Big{]},$

as shown in Figure 1.

Figure 1: The reference tesseract $T^{\ast}$.

The volume is $\text{vol}\left(T^{\ast}\right)=16$ and it is straightforward
to define bounds of integration

$\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 $a=x_{1}$, $b=x_{2}$, $c=x_{3}$, and $d=x_{4}$.

####  2.3.2 Cubic pyramid

Consider the reference cubic pyramid $K^{\ast}$ centered at the origin, having
edge length $l[K^{\ast}]=2$. The vertices are defined such that

$\displaystyle K^{\ast}=\Big{[}$ $\displaystyle
k_{1}\left(-1,-1,-1,-1\right),\;k_{2}\left(-1,1,-1,-1\right),\;k_{3}\left(-1,-1,1,-1\right),\;k_{4}\left(1,-1,-1,-1\right),$
$\displaystyle
k_{5}\left(-1,1,1,-1\right),\;k_{6}\left(1,-1,1,-1\right),\;k_{7}\left(1,1,-1,-1\right),\;k_{8}\left(1,1,1,-1\right),$
$\displaystyle k_{9}\left(0,0,0,0\right)\Big{]},$

as shown in Figure 2.

Figure 2: The reference cubic pyramid $K^{\ast}$.

The volume of the cubic pyramid is $1/8$th of the tesseract, or equivalently
$\text{vol}(K^{\ast})=\frac{Ah}{4}$ where $A=8$ is the area of the hexahedron
base and $h=1$ is the height of the cubic pyramid. Based on the vertices of
$K^{\ast}$, the integration bounds for the cubic pyramid are

$\int_{-1}^{0}\int_{x_{4}}^{-x_{4}}\int_{x_{4}}^{-x_{4}}\int_{x_{4}}^{-x_{4}}\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}\,\mathrm{d}x_{3}\,\mathrm{d}x_{4}=2.$

The orthonormal polynomial basis for the cubic pyramid is found to be

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

where $a=-\frac{x_{1}}{x_{4}}$, $b=-\frac{x_{2}}{x_{4}}$,
$c=-\frac{x_{3}}{x_{4}}$, $d=2x_{4}+1$, and $\mu_{ijk}=2i+2j+2k+3$.

####  2.3.3 Bipentatope

Consider the reference bipentatope $P^{\ast}$ having edge length
$l[P^{\ast}]=1$. This element is also frequently called the tetrahedral
bipyramid. The vertices are defined such that

$\displaystyle P^{\ast}=\Big{[}$ $\displaystyle
p_{1}\left(\tfrac{1}{2},\tfrac{1}{2},\tfrac{1}{2},-\tfrac{1}{2}\right),\;p_{2}\left(\tfrac{1}{2},\tfrac{1}{2},-\tfrac{1}{2},-\tfrac{1}{2}\right),\;p_{3}\left(-\tfrac{1}{2},\tfrac{1}{2},-\tfrac{1}{2},-\tfrac{1}{2}\right),\;p_{4}\left(-\tfrac{1}{2},\tfrac{1}{2},\tfrac{1}{2},-\tfrac{1}{2}\right),$
$\displaystyle
p_{5}\left(0,1,0,-1\right),\;p_{6}\left(0,0,0,-1\right)\Big{]},$

as shown in Figure 3.

Figure 3: The reference bipentatope $P^{\ast}$.

The reference bipentatope volume is
$\text{vol}\left(P^{\ast}\right)=\tfrac{1}{24}$.

In order to facilitate interpolation and integration, the reference
bipentatope can be divided into two pentatopes, as shown in Figure 4.

Figure 4: A subdivision of the reference bipentatope into two pentatopes.

There exists a convenient mapping from these pentatopes onto the standard
pentatope $S^{\ast}$ shown in Figure 5.

Figure 5: The standard pentatope $S^{\ast}$.

We note that $S^{\ast}$ has the following vertices

$\displaystyle S^{\ast}=$
$\displaystyle\Big{[}s_{1}\left(1,-1,-1,-1\right),s_{2}\left(-1,1,-1,-1\right),s_{3}\left(-1,-1,1,-1\right),s_{4}\left(-1,-1,-1,1\right),$
$\displaystyle s_{5}\left(-1,-1,-1,-1\right)\Big{]}.$

The standard pentatope volume is
$\text{vol}\left(S^{\ast}\right)=\tfrac{2}{3}$ and the bounds of integration
are

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

The orthonormal polynomial basis takes the following form

$\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},$
(2)

where $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$, and
$d=x_{4}$.

Lastly, we consider the following rescaled and shifted version of the standard
pentatope

$\displaystyle\mathcal{S}^{\ast}=\left[s_{1}^{*}(1,0,0,0),s_{2}^{*}(0,1,0,0),s_{3}^{*}(0,0,1,0),s_{4}^{*}(0,0,0,1),s_{5}^{*}(0,0,0,0)\right],$

where by inspection

$\displaystyle\mathcal{S}^{\ast}=\tfrac{1}{2}S^{\ast}+\left(1,1,1,1\right).$

This reference element has been used by many authors brandts2007simplicial ;
petrov2018stable  due to its mathematical simplicity and convenience. It is
frequently referred to as the ‘cube corner’. Its volume is
$\text{vol}\left(\mathcal{S}^{\ast}\right)=\tfrac{1}{24}$ and the bounds of
integration are

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

The orthonormal polynomial basis on the cube corner can be obtained by a
simple linear transformation of the basis in Eq. (2).

##  3 Hybrid meshes of tesseracts, cubic pyramids, and bipentatopes

In this section, we introduce a set of hybrid meshes by constructing the
following domain model. Let $\Omega$ be a four-dimensional canonical domain.
The domain $\Omega$ is partitioned into two subdomains $\Omega_{1}$ and
$\Omega_{2}$, which are also canonical. We suppose that $\check{\tau}_{0}$ and
$\hat{\tau}_{0}$ are tesseract triangulations of $\Omega_{1}$ and $\Omega_{2}$
such that all elements of both triangulations form a conforming tesseract
triangulation $\tau_{0}$ of the domain $\Omega$. We define the boundary layers
in the subdomains $\Omega_{1}$ and $\Omega_{2}$ by

$\displaystyle B_{1}$ $\displaystyle=\\{T\in\Omega_{1}\ |\ \dim\left(\partial
T\cap\Omega_{12}\right)=3\\},$ $\displaystyle B_{2}$
$\displaystyle=\\{T\in\Omega_{2}\ |\ \dim\left(\partial
T\cap\Omega_{12}\right)=3\\},\quad\Omega_{12}=\Omega_{1}\cap\Omega_{2}.$

Furthermore, we construct a hybrid mesh so that $\Omega_{1}$ is partitioned by
tesseract elements and $\Omega_{2}$ is partitioned by cubic pyramid elements.
Our main goal is to establish a conforming coupling between the different
kinds of elements in the boundary layers of both subdomains. Additionally, we
need a refinement strategy for both groups of elements with as small as
possible number of classes and as small as possible measure of degeneracy for
the pyramidal elements. We restrict ourselves to refinement strategies
dividing the edges of each element into two parts. We anticipate that such
strategies will be very useful in multigrid and/or adaptive-quadrature
procedures.

We will now define a refinement algorithm for creating nested hierarchical
hybrid triangulations in five steps:

  * $\left(r_{1}\right)$

The tesseract triangulations of $\Omega$, denoted by $\check{\tau}_{0}$,
$\hat{\tau}_{0}$ and $\tau_{0}=\check{\tau}_{0}\cup\hat{\tau}_{0}$ are
constructed. We denote this partition operator by $\mathcal{R}$.

  * $\left(r_{2}\right)$

Each element $T$ of $\check{\tau}_{0}$ is uniformly divided into sixteen
smaller tesseracts. We denote this partition operator by $\mathcal{E}$.
Evidently, this operator’s definition is straightforward, and does not need to
be explicitly stated.

  * $\left(r_{3}\right)$

Each element $T$ of $\hat{\tau}_{0}$ is divided into eight cubic pyramids by
connecting the cubic facets of each element to its centroid. We denote this
partition operator by $\mathcal{B}$. It is given explicitly by

$\displaystyle\mathcal{B}T^{\ast}=\Big{\\{}$ $\displaystyle
K_{1}[t_{1},t_{2},t_{3},t_{4},t_{5},t_{6},t_{7},t_{8},t_{0}],\
K_{2}[t_{1},t_{3},t_{4},t_{6},t_{9},t_{11},t_{12},t_{14},t_{0}],$
$\displaystyle
K_{3}[t_{9},t_{10},t_{11},t_{12},t_{13},t_{14},t_{15},t_{16},t_{0}],\
K_{4}[t_{2},t_{5},t_{7},t_{8},t_{10},t_{13},t_{15},t_{16},t_{0}],$
$\displaystyle
K_{5}[t_{4},t_{6},t_{7},t_{8},t_{12},t_{14},t_{15},t_{16},t_{0}],\
K_{6}[t_{1},t_{2},t_{3},t_{5},t_{9},t_{10},t_{11},t_{13},t_{0}],$
$\displaystyle
K_{7}[t_{1},t_{2},t_{4},t_{7},t_{11},t_{13},t_{14},t_{16},t_{0}],\
K_{8}[t_{3},t_{5},t_{6},t_{8},t_{9},t_{10},t_{12},t_{15},t_{0}]\Big{\\}}.$

One of the resulting cubic pyramids $K_{7}$ is shown in Figure 6. Note:
$K_{7}$ is equivalent to $K^{\ast}$, which is a canonical pyramid.

  * $\left(r_{4}\right)$

The elements of $\mathcal{B}\hat{\tau}_{0}$ are further subdivided by the
partition operator $\mathcal{L}$, which is defined below.

  * $\left(r_{5}\right)$

Steps 2–4 are performed repeatedly until the desired level of mesh refinement
is obtained. These refinement steps are denoted by the operator $\mathcal{H}$.

Figure 6: The cubic pyramid $K_{7}$ belonging to the decomposition
$\mathcal{B}T^{\ast}$.

Now, let us turn our attention to formulating an explicit definition for the
partition operator $\mathcal{L}$

$\displaystyle\mathcal{L}X=\begin{cases}\mathcal{D}X,&\mbox{if }X\;\mbox{is a
cubic pyramid}\\\ \mathcal{M}X,&\mbox{if }X\;\mbox{is a
bipentatope.}\end{cases}$

It remains for us to define the operators $\mathcal{D}$ and $\mathcal{M}$. The
operator $\mathcal{D}$ can be defined as follows

$\displaystyle\mathcal{D}K^{*}=\bigg{\\{}\Big{\\{}$
$\displaystyle\mathbb{K}_{1}[k_{1},k_{10},k_{11},k_{12},k_{13},k_{14},k_{15},k_{16},k_{21}],\;\mathbb{K}_{2}[k_{2},k_{10},k_{13},k_{14},k_{15},k_{17},k_{18},k_{19},k_{20}],$
$\displaystyle\mathbb{K}_{3}[k_{4},k_{12},k_{14},k_{15},k_{21},k_{26},k_{27},k_{28},k_{29}],\;\mathbb{K}_{4}[k_{7},k_{14},k_{15},k_{18},k_{19},k_{27},k_{28},k_{34},k_{35}],$
$\displaystyle\mathbb{K}_{5}[k_{3},k_{11},k_{13},k_{15},k_{21},k_{22},k_{23},k_{24},k_{25}],\;\mathbb{K}_{6}[k_{5},k_{13},k_{15},k_{17},k_{19},k_{22},k_{24},k_{30},k_{31}],$
$\displaystyle\mathbb{K}_{7}[k_{6},k_{15},k_{21},k_{23},k_{24},k_{26},k_{28},k_{32},k_{33}],\;\mathbb{K}_{8}[k_{8},k_{15},k_{19},k_{24},k_{28},k_{30},k_{32},k_{34},k_{36}],$
$\displaystyle\mathbb{K}_{9}[k_{9},k_{16},k_{20},k_{25},k_{29},k_{31},k_{33},k_{35},k_{36}],\;\mathbb{K}_{10}[k_{15},k_{16},k_{20},k_{25},k_{29},k_{31},k_{33},k_{35},k_{36}]\Big{\\}},$
$\displaystyle\Big{\\{}$ $\displaystyle
P_{1}[k_{15},k_{21},k_{23},k_{24},k_{25},k_{33}],\;P_{2}[k_{11},k_{13},k_{15},k_{16},k_{21},k_{25}],\;P_{3}[k_{12},k_{14},k_{15},k_{16},k_{21},k_{29}],$
$\displaystyle
P_{4}[k_{15},k_{21},k_{26},k_{28},k_{29},k_{33}],\;P_{5}[k_{15},k_{24},k_{28},k_{32},k_{33},k_{36}],\;P_{6}[k_{14},k_{15},k_{27},k_{28},k_{29},k_{35}],$
$\displaystyle
P_{7}[k_{15},k_{19},k_{28},k_{34},k_{35},k_{36}],\;P_{8}[k_{13},k_{15},k_{17},k_{19},k_{20},k_{31}],\;P_{9}[k_{15},k_{19},k_{24},k_{30},k_{31},k_{36}],$
$\displaystyle
P_{10}[k_{14},k_{15},k_{18},k_{19},k_{20},k_{35}],\;P_{11}[k_{13},k_{15},k_{22},k_{24},k_{25},k_{31}],\;P_{12}[k_{10},k_{13},k_{14},k_{15},k_{16},k_{20}],$
$\displaystyle
P_{13}[k_{15},k_{19},k_{20},k_{31},k_{35},k_{36}],\;P_{14}[k_{15},k_{28},k_{29},k_{33},k_{35},k_{36}],\;P_{15}[k_{15},k_{16},k_{21},k_{25},k_{29},k_{33}],$
$\displaystyle
P_{16}[k_{13},k_{15},k_{16},k_{20},k_{25},k_{31}],\;P_{17}[k_{15},k_{24},k_{25},k_{31},k_{33},k_{36}],\;P_{18}[k_{14},k_{15},k_{16},k_{20},k_{29},k_{35}]\Big{\\}}\bigg{\\}}.$

All elements $\mathbb{K}_{i}$, $i=1,2,\ldots,10$, are cubic pyramids and all
elements $P_{j}$, $j=1,2,\ldots,18$, are bipentatopes. The first nine vertices
(above) $k_{1},k_{2},\ldots,k_{9}$ are previously defined in the definition of
$K^{\ast}$. In addition, the remaining vertices $k_{10},\ldots,k_{36}$ can be
defined as follows

$\displaystyle k_{10}=\tfrac{1}{2}(k_{1}+k_{2}),\qquad$ $\displaystyle
k_{11}=\tfrac{1}{2}(k_{1}+k_{3}),\qquad$ $\displaystyle
k_{12}=\tfrac{1}{2}(k_{1}+k_{4}),\qquad$ $\displaystyle
k_{13}=\tfrac{1}{2}(k_{1}+k_{5}),$ $\displaystyle
k_{14}=\tfrac{1}{2}(k_{1}+k_{7}),\;$ $\displaystyle
k_{15}=\tfrac{1}{2}(k_{1}+k_{8}),\;$ $\displaystyle
k_{16}=\tfrac{1}{2}(k_{1}+k_{9}),\;$ $\displaystyle
k_{17}=\tfrac{1}{2}(k_{2}+k_{5}),$ $\displaystyle
k_{18}=\tfrac{1}{2}(k_{2}+k_{7}),\;$ $\displaystyle
k_{19}=\tfrac{1}{2}(k_{2}+k_{8}),\;$ $\displaystyle
k_{20}=\tfrac{1}{2}(k_{2}+k_{9}),\;$ $\displaystyle
k_{21}=\tfrac{1}{2}(k_{3}+k_{4}),$ $\displaystyle
k_{22}=\tfrac{1}{2}(k_{3}+k_{5}),\;$ $\displaystyle
k_{23}=\tfrac{1}{2}(k_{3}+k_{6}),\;$ $\displaystyle
k_{24}=\tfrac{1}{2}(k_{3}+k_{8}),\;$ $\displaystyle
k_{25}=\tfrac{1}{2}(k_{3}+k_{9}),$ $\displaystyle
k_{26}=\tfrac{1}{2}(k_{4}+k_{6}),\;$ $\displaystyle
k_{27}=\tfrac{1}{2}(k_{4}+k_{7}),\;$ $\displaystyle
k_{28}=\tfrac{1}{2}(k_{4}+k_{8}),\;$ $\displaystyle
k_{29}=\tfrac{1}{2}(k_{4}+k_{9}),$ $\displaystyle
k_{30}=\tfrac{1}{2}(k_{5}+k_{8}),\;$ $\displaystyle
k_{31}=\tfrac{1}{2}(k_{5}+k_{9}),\;$ $\displaystyle
k_{32}=\tfrac{1}{2}(k_{6}+k_{8}),\;$ $\displaystyle
k_{33}=\tfrac{1}{2}(k_{6}+k_{9}),$ $\displaystyle
k_{34}=\tfrac{1}{2}(k_{7}+k_{8}),\;$ $\displaystyle
k_{35}=\tfrac{1}{2}(k_{7}+k_{9}),\;$ $\displaystyle
k_{36}=\tfrac{1}{2}(k_{8}+k_{9}).$

The full set of vertices is illustrated in Figure 7. In addition, a few of the
individual elements of $\mathcal{D}K^{\ast}$ are shown in Figure 8. Finally,
several couplings between the cubic pyramid and bipentatope elements are shown
in Figure 9.

Figure 7: The cubic pyramid $K^{\ast}$ and the full set of vertices for the
subdivision operator $\mathcal{D}K^{\ast}$.

Figure 8: The cubic pyramid element $\mathbb{K}_{1}$ and the bipentatope
element $P_{13}$ of the decomposition $\mathcal{D}K^{\ast}$.

Figure 9: Couplings between cubic pyramid and bipentatope elements:
$\mathbb{K}_{10}$ and $P_{14}$ (top) and $\mathbb{K}_{2}$ and $P_{10}$
(bottom).

Next, the operator $\mathcal{M}$ can be defined as follows

$\displaystyle\mathcal{M}P^{*}=\bigg{\\{}$
$\displaystyle\mathbb{P}_{1}[p_{5},p_{10},p_{13},p_{16},p_{18},p_{20}],\;\mathbb{P}_{2}[p_{6},p_{11},p_{14},p_{17},p_{19},p_{20}],\;\mathbb{P}_{3}[p_{8},p_{10},p_{13},p_{16},p_{18},p_{20}],$
$\displaystyle\mathbb{P}_{4}[p_{8},p_{11},p_{14},p_{17},p_{19},p_{20}],\;\mathbb{P}_{5}[p_{8},p_{16},p_{17},p_{18},p_{19},p_{20}],\;\mathbb{P}_{6}[p_{8},p_{13},p_{14},p_{16},p_{17},p_{20}],$
$\displaystyle\mathbb{P}_{7}[p_{8},p_{10},p_{11},p_{18},p_{19},p_{20}],\;\mathbb{P}_{8}[p_{8},p_{10},p_{11},p_{13},p_{14},p_{20}],\;\mathbb{P}_{9}[p_{3},p_{8},p_{12},p_{15},p_{16},p_{17}],$
$\displaystyle\mathbb{P}_{10}[p_{4},p_{8},p_{9},p_{15},p_{18},p_{19}],\;\mathbb{P}_{11}[p_{1},p_{7},p_{8},p_{9},p_{10},p_{11}],\;\mathbb{P}_{12}[p_{2},p_{7},p_{8},p_{12},p_{13},p_{14}],$
$\displaystyle\mathbb{P}_{13}[p_{8},p_{15},p_{16},p_{17},p_{18},p_{19}],\;\mathbb{P}_{14}[p_{8},p_{12},p_{13},p_{14},p_{16},p_{17}],\;\mathbb{P}_{15}[p_{7},p_{8},p_{10},p_{11},p_{13},p_{14}],$
$\displaystyle\mathbb{P}_{16}[p_{8},p_{9},p_{10},p_{11},p_{18},p_{19}]\bigg{\\}}.$

All the elements $\mathbb{P}_{j}$, $j=1,2,\ldots,16$, are bipentatopes. The
first six vertices (above) $p_{1},p_{2},\ldots,p_{6}$ are previously defined
in the definition of $P^{\ast}$. In addition, the remaining vertices
$p_{7},\ldots,p_{20}$ can be defined as follows

$\displaystyle p_{7}=\tfrac{1}{2}(p_{1}+p_{2}),\qquad$ $\displaystyle
p_{8}=\tfrac{1}{2}(p_{1}+p_{3}),\qquad$ $\displaystyle
p_{9}=\tfrac{1}{2}(p_{1}+p_{4}),\qquad$ $\displaystyle
p_{10}=\tfrac{1}{2}(p_{1}+p_{5}),$ $\displaystyle
p_{11}=\tfrac{1}{2}(p_{1}+p_{6}),\;$ $\displaystyle
p_{12}=\tfrac{1}{2}(p_{2}+p_{3}),\;$ $\displaystyle
p_{13}=\tfrac{1}{2}(p_{2}+p_{5}),\;$ $\displaystyle
p_{14}=\tfrac{1}{2}(p_{2}+p_{6}),$ $\displaystyle
p_{15}=\tfrac{1}{2}(p_{3}+p_{4}),\;$ $\displaystyle
p_{16}=\tfrac{1}{2}(p_{3}+p_{5}),\;$ $\displaystyle
p_{17}=\tfrac{1}{2}(p_{3}+p_{6}),\;$ $\displaystyle
p_{18}=\tfrac{1}{2}(p_{4}+p_{5}),$ $\displaystyle
p_{19}=\tfrac{1}{2}(p_{4}+p_{6}),$ $\displaystyle
p_{20}=\tfrac{1}{2}(p_{5}+p_{6}).$

The full set of vertices is illustrated in Figure 10. In addition, a few of
the individual elements of $\mathcal{M}P^{\ast}$ are shown in Figure 11.

Figure 10: The bipentatope $P^{\ast}$ and the full set of vertices for the
subdivision operator $\mathcal{M}P^{\ast}$.

Figure 11: The elements $\mathbb{P}_{2}$, $\mathbb{P}_{4}$, and
$\mathbb{P}_{9}$ of the decomposition $\mathcal{M}P^{\ast}$.

We conclude this section by summarizing the refinement procedure for obtaining
hierarchical triangulations below:

$\displaystyle\tau_{0}=\mathcal{R}\Omega,\quad\tau_{0}=\check{\tau}_{0}\cup\hat{\tau}_{0};$
$\displaystyle\hat{\tau}_{0}^{*}=\mathcal{B}\hat{\tau}_{0},\quad\tau_{0}^{*}=\check{\tau}_{0}\cup\hat{\tau}_{0}^{*};$
$\displaystyle\check{\tau}_{1}=\mathcal{E}\check{\tau}_{0},\;\hat{\tau}_{1}=\mathcal{L}\hat{\tau}_{0}^{*},\quad\tau_{1}=\check{\tau}_{1}\cup\hat{\tau}_{1};$
$\check{\tau}_{2}=\mathcal{E}\check{\tau}_{1},\;\hat{\tau}_{2}=\mathcal{L}\hat{\tau}_{1},\quad\tau_{2}=\check{\tau}_{2}\cup\hat{\tau}_{2};$
(3)
$\displaystyle\check{\tau}_{k+1}=\mathcal{E}^{k}\check{\tau}_{1},\;\hat{\tau}_{k+1}=\mathcal{L}^{k}\hat{\tau}_{1},\quad\tau_{k+1}=\check{\tau}_{k+1}\cup\hat{\tau}_{k+1}.$

In accordance with (3), we can explicitly define the partition operator
$\mathcal{H}:\tau_{1}\rightarrow\tau_{2}$ as follows

$\displaystyle\mathcal{H}X=\begin{cases}\mathcal{E}X,&\mbox{if }X\;\mbox{is a
tesseract}\\\ \mathcal{L}X,&\mbox{if }X\;\mbox{is a cubic pyramid or
bipentatope.}\end{cases}$

Then

$\displaystyle\tau_{k+1}=\mathcal{H}^{k}\tau_{1}.$

Thus, we have a conforming sequence of successive hybrid finite element
triangulations $\left\\{\tau_{k}\right\\}$.

##  4 Theoretical results

In this section, we summarize the important theoretical results that govern
the hybrid triangulations introduced in the previous section.

###### Property 1

The hypercube $T^{n}$ can be divided into $2n$ hypercubic pyramids $K^{n}_{i}$
so that all edges of the $i$-th pyramid have the same length.

###### Theorem 1

The four-dimensional space is the only Euclidean space that satisfies Property
1.

###### Proof

For the sake of simplicity, we suppose that

$\displaystyle T^{n}=\left[t_{1}^{n}(-a,-a,\ldots,-a),\ldots
t_{2^{n}}^{n}(a,a,\ldots,a)\right]$

has all edges parallel to the coordinate axes. Then, Property 1 is valid iff

$\displaystyle\frac{1}{2}h[T^{n}]=l[T^{n}],$

where $l[T^{n}]$ is the edge length of $T^{n}$. The latter equality leads to

$\displaystyle\frac{1}{2}\sqrt{\sum_{k=1}^{n}\left(2a\right)^{2}}=2a\Rightarrow\sqrt{na^{2}}=2a,$

which is possible iff $n=4$. ∎

The result of Theorem 1 is that the four-dimensional space is unique among all
Euclidean spaces.

###### Theorem 2

Consider the cubic pyramids in $\mathcal{D}K^{\ast}$. All such pyramids are
congruent to the reference pyramid, i.e. $\mathbb{K}_{i}\subset[K^{*}],$
$i=1,2,\ldots,10$.

###### Proof

For all elements $\mathbb{K}_{i}$ $i=1,2,\ldots 10,$ of $\mathcal{D}K^{*}$ it
is possible to verify through direct calculation that:

  1. (i)

all edges of $\mathbb{K}_{i}$ are equal to $1$;

  2. (ii)

the base of $\mathbb{K}_{i}$ is a cube with edges equal to $1$;

  3. (iii)

the degeneracy measure of $\mathbb{K}_{i}$ is
$\delta(\mathbb{K}_{i})=4.18154$;

from which congruency with the reference pyramid follows. ∎

###### Theorem 3

Consider the bipentatopes in $\mathcal{D}K^{\ast}$. All such bipentatopes are
congruent.

###### Proof

Without loss of generality, we consider demonstrating this property for two
bipentatopes chosen from $\mathcal{D}K^{*}$, with the claim following from
repeating this procedure for all pairwise combinations. With this in mind, we
take our two bipentatopes chosen from $\mathcal{D}K^{*}$ to be $P_{10}$ and
$P_{14}$, see Figure 9. Any bipentatope can be partitioned into two pentatopes
as shown in Figure 4. We hence begin by dividing the bipentatopes $P_{10}$ and
$P_{14}$ into pentatopes whence

$\displaystyle
P_{10}=\hat{S}_{10}\left[p_{15},p_{18},p_{19},p_{20},p_{35}\right]\cup\check{S}_{10}\left[p_{14},p_{15},p_{18},p_{20},p_{35}\right],$
$\displaystyle
P_{14}=\hat{S}_{14}\left[p_{15},p_{28},p_{29},p_{33},p_{36}\right]\cup\check{S}_{14}\left[p_{15},p_{28},p_{29},p_{35},p_{36}\right].$

Let the generic affine transformations of the pentatopes $\hat{S}_{10}$ and
$\hat{S}_{14}$ be

$\displaystyle F_{10}\;:\;\hat{S}_{10}=A_{10}\mathcal{S}^{*}+B_{10},\quad
F_{14}\;:\;\hat{S}_{14}=A_{14}\mathcal{S}^{*}+B_{14}.$

We map $\hat{S}_{14}$ onto $\hat{S}_{10}$ by

$\displaystyle\hat{S}_{10}=F_{10}\left(F_{14}^{-1}\left(\hat{S}_{14}\right)\right)\Leftrightarrow
F_{14,10}\;:\;\hat{S}_{10}=A_{14,10}\hat{S}_{14}+B_{10}-A_{14,10}B_{14},$

where $A_{14,10}$ is the transitional matrix from $\hat{S}_{14}$ to
$\hat{S}_{10}$ and $F_{14,10}=F_{10}\circ F_{14}^{-1}$. We emphasize that
$p_{14}=F_{14,10}(p_{35})$, which assures that
$P_{10}=F_{14,10}\left(P_{14}\right)$. The claim that $P_{10}\cong P_{14}$
follows from the fact that the transition matrix

$\displaystyle A_{14,10}=\left(\begin{array}[]{cccc}-1&0&0&0\\\ 0&0&1&0\\\
0&-1&0&0\\\ 0&0&0&-1\\\ \end{array}\right)$

is orthogonal. ∎

###### Remark 1

Let all bipentatopes from $\mathcal{D}K^{*}$ be partitioned into two
pentatopes as it is shown in Figure 4. Then all of them belong to
$[\hat{S}_{10}]$ and $\delta\left([\hat{S}_{10}]\right)=5.$ Moreover, the
elements of $[\hat{S}_{10}]$ are invariant. For more details concerning this
class of pentatopes, we refer the reader to petrov2020properties .

###### Theorem 4

Consider the bipentatopes in $\mathcal{M}P^{\ast}$. All such bipentatopes are
congruent to the reference bipentatope, i.e.
$\mathbb{P}_{j}\subset[P^{\ast}],$ $j=1,2,\ldots,16$.

###### Proof

Firstly, we prove that all bipentatopes in $\mathcal{M}P^{*}$ belong to the
same class. We randomly choose two bipentatopes $\mathbb{P}_{2}$ and
$\mathbb{P}_{4}$ (see Figure 11) from $\mathcal{M}P^{*}$. The bipentatope
$\mathbb{P}_{2}$ is obtained by the coupling of the pentatopes

$\displaystyle\hat{S}_{2}\left[p_{6},p_{11},p_{14},p_{17},p_{20}\right]\quad\mbox{and}\quad\check{S}_{2}\left[p_{6},p_{11},p_{17},p_{19},p_{20}\right].$

The second bipentatope is partitioned as follows

$\displaystyle\mathbb{P}_{4}=\hat{S}_{4}\left[p_{8},p_{11},p_{14},p_{17},p_{20}\right]\cup\check{S}_{4}\left[p_{8},p_{11},p_{17},p_{19},p_{20}\right].$

The pentatope $\hat{S}_{2}$ is mapped onto $\hat{S}_{4}$ by the affine
transformation

$\hat{S}_{4}=F_{4}\left(F_{2}^{-1}\left(\hat{S}_{2}\right)\right),$

where

$\displaystyle F_{2}\;:\;\hat{S}_{2}=A_{2}\mathcal{S}^{*}+B_{2},\quad
F_{4}\;:\;\hat{S}_{4}=A_{4}\mathcal{S}^{*}+B_{4}.$

We present the map $F_{2,4}=F_{4}\circ F_{2}^{-1}$ in a matrix form by

$\displaystyle F_{2,4}\;:\;\hat{S}_{4}=A_{2,4}\hat{S}_{2}+B_{4}-A_{2,4}B_{4}.$

The transitional matrix

$\displaystyle A_{2,4}=\left(\begin{array}[]{cccc}1&0&0&0\\\ 0&0&0&-1\\\
0&0&1&0\\\ 0&-1&0&0\\\ \end{array}\right)$

is orthogonal. Additionally, $F_{2,4}$ keeps the facet $\mathbb{P}_{2,-1,-6}$
invariant, which assures $\mathbb{P}_{4}=F_{2,4}\left(\mathbb{P}_{2}\right)$.
Therefore, $\mathbb{P}_{2}\cong\mathbb{P}_{4}$.

On the other hand, all edges of the bipentatopes have the same length equal to
$\frac{1}{2}$, i.e. they are identical.

It remains to prove that $\mathbb{P}_{2}$ belongs to $[P^{*}]$. This follows
directly from $A_{2,*}=2I$, where $A_{2,*}$ is the transitional matrix from
$\mathbb{P}_{2}$ to $P^{*}$ and $I$ is the identity matrix. ∎

##  5 Quadrature rules for the cubic pyramid

###  5.1 Methodology

In this section, we consider the problem of approximating the integral of a
function $f(\boldsymbol{x})$ on $K^{\ast}$, the reference cubic pyramid. This
can be effectively accomplished via a quadrature rule whence

$\int_{-1}^{0}\int_{x_{4}}^{-x_{4}}\int_{x_{4}}^{-x_{4}}\int_{x_{4}}^{-x_{4}}f(\boldsymbol{x})\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}\,\mathrm{d}x_{3}\,\mathrm{d}x_{4}\approx\sum_{j=1}^{N}\omega_{j}f(\boldsymbol{x}_{j}),$
(4)

where $\\{\omega_{j}\\}$ are a set of $N$ weights and
$\\{\boldsymbol{x}_{j}\\}$ are an associated set of $N$ points—known as
abscissa—which together specify a _quadrature rule_. If a rule is capable of
exactly computing the integrals of a generic constant and all polynomials in
the set

$\Xi(p)=\\{\psi_{ijkq}(\boldsymbol{x})\mid 0<i+j+k+q\leq p\text{ and
}i,j,k,q\geq 0\\},$

where $\psi_{ijkq}(\boldsymbol{x})$ is as given in Eq. (1), then the rule is
said to be of strength $p$. When defining the abscissa of a rule, it is
typical to constrain them to be strictly inside of our domain and be arranged
symmetrically. This ensures that we never evaluate $f(\boldsymbol{x})$ outside
of its domain. In addition, to reduce the likelihood of catastrophic
cancellation it is also customary to require the weights to be positive.
Evidently if the rule is to integrate a generic constant mode correctly it
must necessarily be the case that

$\sum_{j=1}^{N}\omega_{j}=2.$

Substituting the polynomials from our set into Eq. (4) and enforcing equality
we require for all $\psi\in\Xi(p)$ that

$\int_{-1}^{0}\int_{x_{4}}^{-x_{4}}\int_{x_{4}}^{-x_{4}}\int_{x_{4}}^{-x_{4}}\psi(\boldsymbol{x})\,\mathrm{d}x_{1}\,\mathrm{d}x_{2}\,\mathrm{d}x_{3}\,\mathrm{d}x_{4}=0=\sum_{j=1}^{N}\omega_{j}\psi(\boldsymbol{x}_{j}),$

where in evaluating the integral we have exploited the orthogonality
relationship associated with the polynomials in our basis. This can be
regarded as a non-linear least squares problem in which we have $|\Xi(p)|+1$
equations and $5N$ unknowns. However, unlike a typical least squares problem,
it is one where we seek a solution with zero residual. Assuming the abscissa
are known we remark that the associated weights may be determined through the
solution of a _linear_ least squares problem. This enables us to reduce the
number of degrees of freedom in our non-linear problem to $4N$.

In order to enforce symmetry, it is advantageous to decompose our $N$ abscissa
into symmetry orbits. The symmetry orbits of the cubic pyramid are similar to
those of the hexahedron, with an additional parameter that accounts for the
extent of the pyramid in the fourth dimension. Therefore, one may follow the
procedure for generating orbits for the hexahedron (see
witherden2015identification ), and enumerate all symmetry orbits. Thereafter,
upon adding a parameter to account for the fourth dimension, the symmetry
orbits of the cubic pyramid are written as follows

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

subject to the constraints $0<\alpha,\beta,\gamma\leq-\delta$ and
$-1\leq\delta\leq 0$. Here, the operator $\chi$ returns all possible (unique)
signed permutations of the input arguments.

Now, let us provide an example of how these symmetry orbits are used. If
$N=156$, then one possible arrangement of points is four $S_{1}$ orbits, two
$S_{2}$ orbits, one $S_{3}$ orbit, three $S_{4}$ orbits, two $S_{5}$ orbits,
and one $S_{7}$ orbit whence

$N=4|S_{1}|+2|S_{2}|+|S_{3}|+3|S_{4}|+2|S_{5}|+|S_{7}|=156,$

as required. Inspecting the number of arguments for each of the orbits and
summing we find that this decomposition has a total of $26$ degrees of
freedom. Together these fix the locations of all of our abscissa and hence
become the parameters in our non-linear least squares problem. We note here
that this represents a substantial decrease compared with an asymmetric rule
which has some $4N=624$ degrees of freedom. Of course, for a given $N$ there
are typically a very large number of distinct decompositions into symmetry
orbits. For our example of $N=156$ there are some $8518$ unique
decompositions.

Our overall goal when identifying quadrature rules is to find one with the
minimum value of $N$ for a particular strength $p$. Unfortunately, given a $p$
it is not possible to determine $N$ theoretically; instead it must be treated
as a hyper-parameter to the overall optimisation process. Similarly, once
given an $N$ it is, in general, not possible to know _a priori_ which
particular symmetric decomposition will yield a rule–if any. Thus, it is also
necessary to treat the specific decomposition as a hyper-parameter, too.

We have implemented the aforementioned system into the open source symmetric
quadrature rule package Polyquad witherden2015identification . Running this
modified version of Polyquad on a workstation, we have identified symmetric
quadrature rules on the cubic pyramid for $2\leq p\leq 12$. The total number
of points required for these rules can be seen in Table 1. The abscissa for
all of our rules are inside of our reference pyramid, and feature positive
weights. The rules themselves are provided in quadruple precision in the
electronic supplemental material. An example of the $p=5$ quadrature rule is
provided in Table 2.

Table 1: Number of points $N$ required for a rule of strength $p$ inside of the cubic pyramid. $p$ | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---|---|---|---|---|---|---
$N$ | 7 | 10 | 21 | 29 | 50 | 65 | 114 | 163 | 251 | 323 | 552
Table 2: The $p=5$ quadrature rule on the cubic pyramid $K^{\ast}$ with $N=29$ points. The remaining rules are provided in the electronic supplemental material. Orbit | Abscissa | Weight
---|---|---
$S_{1}$ | -0.32846339581882957902277574789024 | 0.084184481885656624083314191434215
$S_{2}$ | 0.59942266549153142521578075246447 | 0.11340384654119720346152792514806
| -0.74912930534502345701037962143888 |
$S_{2}$ | 0.76070275297226284532398968584027 | 0.089004419809134534555110633350815
| -0.95754362470717092670167055158225 |
$S_{3}$ | 0.43384835705078179649267535015814 | 0.033983660519520996993103370708523
| -0.5766200661792940506119682920933 |
$S_{3}$ | 0.69719510426558290268820473081072 | 0.05368707948202312148400343648805
| -0.91852788682445245487884982935091 |

###  5.2 Numerical experiments

The objective of this section is to numerically validate a subset of the
quadrature rules in Table 1. Towards this end, we performed a series of
numerical experiments on the rules with strengths $p=2,3,\ldots,9$. These
rules are capable of exactly integrating certain monomial functions by
construction, as well as approximately integrating certain transcendental
functions. In what follows, we limit our focus to: a) weighted combinations of
monomial functions, and b) analytic transcendental functions. In addition, we
execute all of the numerical experiments in quadruple precision arithmetic in
order to demonstrate the full precision of the rules.

####  5.2.1 Polynomial integration

In this section, we evaluate the ability of our quadrature rules to integrate
a weighted combination of monomial functions on the reference cubic pyramid
$K^{\ast}$. More specifically, we consider the following polynomial function
of degree $m$

$\displaystyle
f_{\text{poly}}\left(\boldsymbol{x};m\right)=\sum_{r=0}^{m}\;\sum_{s=0}^{m-r}\;\sum_{t=0}^{m-r-s}\;\sum_{v=0}^{m-r-
s-t}C_{rstv}\,x_{1}^{r}x_{2}^{s}x_{3}^{t}x_{4}^{v},$

where the constants $C_{rstv}$ are given by the following formula

$\displaystyle
C_{rstv}=24\frac{\left(r+1\right)\left(s+1\right)\left(t+1\right)\left(v+1\right)}{\left(m+1\right)\left(m+2\right)\left(m+3\right)\left(m+4\right)}.$

One may observe that this function contains all possible distinct monomials of
degree $m$. It turns out that the integral of this function on the domain
$K^{\ast}$ can be computed exactly as follows

$\displaystyle
J_{\infty}\left(\boldsymbol{x};m\right)=\int_{K^{\ast}}f_{\text{poly}}\left(\boldsymbol{x}\right)d\boldsymbol{x}=\sum_{r=0}^{m}\;\sum_{s=0}^{m-r}\;\sum_{t=0}^{m-r-s}\;\sum_{v=0}^{m-r-
s-t}C_{rstv}I_{rstv},$

where

$\displaystyle
I_{rstv}=\int_{K^{\ast}}x_{1}^{r}x_{2}^{s}x_{3}^{t}x_{4}^{v}\,d\boldsymbol{x}=\frac{\left(1+\left(-1\right)^{r}\right)\left(1+\left(-1\right)^{s}\right)(1+\left(-1\right)^{t})\left(-1\right)^{v}}{\left(1+r\right)\left(1+s\right)\left(1+t\right)\left(r+s+t+v+4\right)},$

is the exact integral of each monomial.

The approximate integral of $f_{\text{poly}}$ can be computed by evaluating
the following quadrature formula

$\displaystyle
J_{p}=\sum_{j=1}^{N}\omega_{j}f_{\text{poly}}\left(\boldsymbol{x}_{j}\right),$

where the number of points $N=N\left(p\right)$, weights
$\omega_{j}=\omega_{j}\left(p\right)$, and abscissa
$\boldsymbol{x}_{j}=\boldsymbol{x}_{j}\left(p\right)$ are functions of the
quadrature rule strength $p$. Evidently, we obtain exact integration
$J_{p}=J_{\infty}$ when $p\geq m$, and approximate integration $J_{p}\approx
J_{\infty}$ when $p<m$.

Figure 12 shows the quadrature error produced by integrating $f_{\text{poly}}$
for different values of $m$ and $p$. In accordance with expectations, we
obtain exact integration whenever $p\geq m$, to within machine precision
($10^{-32}$).

Figure 12: Absolute error in the numerical integration of $f_{\text{poly}}$
over the reference cubic pyramid $K^{\ast}$. The dashed line marks the
threshold of machine precision.

####  5.2.2 Transcendental integration

In this section, we evaluate the ability of our quadrature rules to integrate
a set of transcendental functions on a family of successively refined meshes.
For this purpose, we consider the following functions

$\displaystyle f_{1}\left(\boldsymbol{x}\right)$ $\displaystyle=\sin\left(\pi
x_{1}^{2}\right)\sin\left(\pi x_{2}^{2}\right)\sin\left(\pi
x_{3}^{2}\right)\sin\left(\pi x_{4}^{2}\right),$ $\displaystyle
f_{2}\left(\boldsymbol{x}\right)$
$\displaystyle=\exp\left(x_{1}^{2}\right)\exp\left(x_{2}^{2}\right)\exp\left(x_{3}^{2}\right)\exp\left(x_{4}^{2}\right),$
$\displaystyle f_{3}\left(\boldsymbol{x}\right)$
$\displaystyle=\exp\left(x_{1}\right)\exp\left(\tfrac{1}{2}x_{2}\right)\exp\left(\tfrac{1}{3}x_{3}\right)\exp\left(\tfrac{1}{4}x_{4}\right).$

Note that the first two functions are symmetric with respect to the
coordinates $x_{1},\ldots,x_{4}$, whereas the last function is asymmetric.

The integral of each transcendental function was computed on the domain
$\Omega=\left[0,1\right]^{4}$ as follows

$\displaystyle
J_{\infty}=\int_{\Omega}f_{\text{trans}}\left(\boldsymbol{x}\right)d\boldsymbol{x}.$

The analytical solution of this integral is generally unknown. As a result,
the ‘exact’ integration was carried out using a vectorized adaptive
integration routine in Matlab. This routine approximates the integral of a
smooth function to an arbitrary level of precision by successively subdividing
the domain of integration, and leveraging nested Gauss-Kronrod quadrature
rules to estimate the integration error on each subinterval. One may consult
the work of Shampine shampine2008matlab ; shampine2008vectorized  for details
of the Matlab implementation, and Notaris notaris2016gauss  for a general
review of Gauss-Kronrod quadrature rules. In our case, we used the adaptive
Matlab routine to approximate the integrals to within machine precision.

The transcendental functions above were also integrated using the
aforementioned quadrature rules on the cubic pyramid. In order to facilitate
this process, the domain $\Omega$ was covered with a uniform mesh of $M^{4}$
tesseract elements where $M$ is a positive integer. Thereafter, each tesseract
was subdivided into 8 cubic pyramid elements in accordance with the partition
operator $\mathcal{B}$, yielding a total of $N_{e}=8M^{4}$ cubic pyramid
elements. The quadrature rules were transformed (mapped) from the reference
cubic pyramid $K^{\ast}$ to the individual cubic pyramids in the mesh via a
linear mapping procedure. The approximate integral of each function was then
computed in the following fashion

$\displaystyle
J_{p}=\sum_{i=1}^{N_{e}}\sum_{j=1}^{N}\omega_{j}^{K_{i}}f_{\text{trans}}\left(\boldsymbol{x}_{j}^{K_{i}}\right),$

where $\omega_{j}^{K_{i}}$ are the quadrature weights and
$\boldsymbol{x}_{j}^{K_{i}}$ are the point locations in each element $K_{i}$.

Figures 13, 14, and 15 illustrate the quadrature errors for integrating
$f_{1}$, $f_{2}$, and $f_{3}$ on $\Omega$ for mesh parameters
$M=1,2,\ldots,10$. In each case, the error converges at a rate of $h^{p+1}$ or
higher. In addition, the pairs of even and odd quadrature rules converge at
the same rates, e.g. rules with $p=2$ and $p=3$ converge at a rate of $h^{4}$.
Upon combining these insights together, we estimate that the rules will
typically converge at a rate of $h^{p+2}$ for even values of $p$, and
$h^{p+1}$ for odd values of $p$. The precise reason for this behavior is
unknown, although it is likely due to fortunate cancellations of truncation
error terms for quadrature rules with even values of $p$.

In addition, we note that increasing the number of quadrature points (by
increasing $p$) does not always yield a more accurate result when integrating
transcendental functions. For example, the 50 point rule with $p=6$ sometimes
outperforms the 65 point rule with $p=7$ (see Figure 15). This behavior is
likely due to a convenient alignment between the symmetry orbits of the $p=6$
rule and the local topology of the transcendental function $f_{3}$. We note
that this trend is not general, as it does not hold for $f_{1}$ and $f_{2}$.
In fact, in most cases the higher degree rules outperform the lower degree
rules, as expected.

Figure 13: Absolute value of the error in the numerical integration of $f_{1}$
over the domain $\Omega$ for mesh parameters $M=1,2,\ldots,10$. Figure 14:
Absolute value of the error in the numerical integration of $f_{2}$ over the
domain $\Omega$ for mesh parameters $M=1,2,\ldots,10$.

Figure 15: Absolute value of the error in the numerical integration of $f_{3}$
over the domain $\Omega$ for mesh parameters $M=1,2,\ldots,10$.

##  6 Conclusion

In this work, we have presented a novel refinement strategy for four-
dimensional cubic pyramids. Given a cubic pyramid our method subdivides it
into a conforming set of smaller cubic pyramids and bipentatopes. Moreover,
all of the edges of each cubic pyramid have the same length, as do the edges
of each bipentatope, which corresponds to an optimal configuration.
Furthermore, we have also developed and evaluated a new set of polynomial
quadrature rules inside of the cubic pyramid. Together, these results open up
new pathways for four-dimensional meshing for space-time finite element
methods.

For the sake of completeness, we also provided a comprehensive example of a
four-dimensional conformal hybrid mesh on a canonical domain. The initial mesh
contains tesseract and cubic pyramid elements. The mesh is conformal because
the couplings between the tesseracts and the cubic pyramids are consistent,
without the need for an interface subdomain. Following the formation of the
initial mesh, it is refined by subdividing the tesseract elements in the
standard fashion (splitting along planes that are orthogonal to the coordinate
axes), and by subdividing the cubic pyramids in the aforementioned novel
fashion. We believe that this strategy can be generalized to a much broader
class of hybrid meshes, including curved meshes, without significant
modifications.

## Declarations

#### Funding

The authors did not receive support from any organization for the submitted
work.

#### Conflict of interest/Competing interests

The authors have no conflicts of interest to declare that are relevant to the
content of this article.

#### Availability of data and material

The quadrature rules are available as electronic supplemental material.

#### Code availability

The Polyquad code is available at https://github.com/PyFR/Polyquad.

## References

  * (1) Ainsworth, M., Fu, G.: A lowest-order composite finite element exact sequence on pyramids.  Computer Methods in Applied Mechanics and Engineering 324, 110–127 (2017)
  * (2) Arnold, D., Falk, R., Winther, R.: Finite element exterior calculus: from Hodge theory to numerical stability.  Bulletin of the American Mathematical Society 47(2), 281–354 (2010)
  * (3) Arnold, D.N., Falk, R.S., Winther, R.: Finite element exterior calculus, homological techniques, and applications.  Acta Numerica 15, 1 (2006)
  * (4) Bedrosian, G.: Shape functions and integration formulas for three-dimensional finite element analysis.  International Journal for Numerical Methods in Engineering 35(1), 95–108 (1992)
  * (5) Bergot, M., Cohen, G., Duruflé, M.: Higher-order finite elements for hybrid meshes using new nodal pyramidal elements.  Journal of Scientific Computing 42(3), 345–381 (2010)
  * (6) Brandts, J., Korotov, S., Křížek, M.: Simplicial finite elements in higher dimensions.  Applications of Mathematics 52(3), 251–265 (2007)
  * (7) Cao, C., Qin, Q.H., Yu, A.: A new hybrid finite element approach for three-dimensional elastic problems.  Archives of Mechanics 64(3), 261–292 (2012)
  * (8) Caplan, P.C.: Four-dimensional anisotropic mesh adaptation for spacetime numerical simulations.  Ph.D. thesis, Massachusetts Institute of Technology (2019)
  * (9) Chan, J., Warburton, T.: A comparison of high order interpolation nodes for the pyramid.  SIAM Journal on Scientific Computing 37(5), A2151–A2170 (2015)
  * (10) Chan, J., Warburton, T.: hp-finite element trace inequalities for the pyramid.  Computers & Mathematics with Applications 69(6), 510–517 (2015)
  * (11) Chan, J., Warburton, T.: Orthogonal bases for vertex-mapped pyramids.  SIAM Journal on Scientific Computing 38(2), A1146–A1170 (2016)
  * (12) Chan, J., Warburton, T.: A short note on a Bernstein–Bezier basis for the pyramid.  SIAM Journal on Scientific Computing 38(4), A2162–A2172 (2016)
  * (13) Coulomb, J.L., Zgainski, F.X., Marechal, Y.: A pyramidal element to link hexahedral, prismatic and tetrahedral edge finite elements.  IEEE Transactions on Magnetics 33(2), 1362–1365 (1997)
  * (14) Coxeter, H.: Regular and semi-regular polytopes. I.  Mathematische Zeitschrift 46(1), 380–407 (1940)
  * (15) Coxeter, H.S.M.: Regular Polytopes.  Courier Corporation (1973)
  * (16) Devloo, P.R., Duran, O., Gomes, S.M., Ainsworth, M.: High-order composite finite element exact sequences based on tetrahedral–hexahedral–prismatic–pyramidal partitions.  Computer Methods in Applied Mechanics and Engineering 355, 952–975 (2019)
  * (17) Dumont, S., Jourdan, F.: A space-time finite element method for elastodynamics problems with a moving loading zone.  In: Eleventh International Conference on Computational Structures Technology (2012)
  * (18) Frontin, C.V., Walters, G.S., Witherden, F.D., Lee, C.W., Williams, D.M., Darmofal, D.L.: Foundations of space-time finite element methods: polytopes, interpolation, and integration.  arXiv preprint arXiv:2012.08701 (2020)
  * (19) Gillette, A.: Serendipity and tensor product affine pyramid finite elements.  The SMAI Journal of Computational Mathematics 2, 215–228 (2016)
  * (20) Korotov, S., Křížek, M.: Red refinements of simplices into congruent subsimplices.  Computers & Mathematics with Applications 67(12), 2199–2204 (2014)
  * (21) Langer, U., Steinbach, O.: Space-Time Methods: Applications to Partial Differential Equations, vol. 25.  Walter de Gruyter GmbH & Company (2019)
  * (22) Maday, Y., Mavriplis, C., Patera, A.: Nonconforming mortar element methods: Application to spectral discretizations.  In: Domain decomposition methods: Proceedings of the Second International Symposium, pp. 392–418 (1988)
  * (23) McLaren-Young-Sommerville, D.: An Introduction to the Geometry of n Dimensions.  Dover (1958)
  * (24) Meshkat, S., Talmor, D.: Generating a mixed mesh of hexahedra, pentahedra and tetrahedra from an underlying tetrahedral mesh.  International Journal for Numerical Methods in Engineering 49(1-2), 17–30 (2000)
  * (25) Neumüller, M., Karabelas, E.: Generating admissable space-time meshes for moving domains in (d+1) dimensions.  In: Space-Time Methods: Applications to Partial Differential Equations, pp. 185–206 (2019)
  * (26) Nigam, N., Phillips, J.: Numerical integration for high order pyramidal finite elements.  ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 46(2), 239–263 (2012)
  * (27) Notaris, S.E.: Gauss-Kronrod quadrature formulae–A survey of fifty years of research.  Electronic Transactions on Numerical Analysis 45, 371–404 (2016)
  * (28) Pantano, A., Averill, R.C.: A penalty-based interface technology for coupling independently modeled 3D finite element meshes.  Finite Elements in Analysis and Design 43(4), 271–286 (2007)
  * (29) Petrov, M.S., Todorov, T.D.: Stable subdivision of 4D polytopes.  Numerical Algorithms 79(2), 633–656 (2018)
  * (30) Petrov, M.S., Todorov, T.D.: Properties of the multidimensional finite elements.  Applied Mathematics and Computation (2020).  In Press
  * (31) Reberol, M., Lévy, B.: Low-order continuous finite element spaces on hybrid non-conforming hexahedral-tetrahedral meshes.  arXiv preprint arXiv:1605.02626 (2016)
  * (32) Seshaiyer, P., Suri, M.: hp submeshing via non-conforming finite element methods.  Computer Methods in Applied Mechanics and Engineering 189(3), 1011–1030 (2000)
  * (33) Shampine, L.F.: Matlab program for quadrature in 2D.  Applied Mathematics and Computation 202(1), 266–274 (2008)
  * (34) Shampine, L.F.: Vectorized adaptive quadrature in Matlab.  Journal of Computational and Applied Mathematics 211(2), 131–140 (2008)
  * (35) Steinbach, O.: Space-time finite element methods for parabolic problems.  Computational Methods in Applied Mathematics 15(4), 551–566 (2015)
  * (36) Steinbach, O., Yang, H.: An algebraic multigrid method for an adaptive space–time finite element discretization.  In: International Conference on Large-Scale Scientific Computing, pp. 66–73. Springer (2017)
  * (37) Todorov, T.D.: The optimal refinement strategy for 3-D simplicial meshes.  Computers & Mathematics with Applications 66(7), 1272–1283 (2013)
  * (38) Toulopoulos, I.: Space-time finite element methods stabilized using bubble function spaces.  Applicable Analysis 99(7), 1153–1170 (2020)
  * (39) Williams, D.M., Frontin, C.V., Miller, E.A., Darmofal, D.L.: A family of symmetric, optimized quadrature rules for pentatopes.  Computers & Mathematics with Applications 80(5), 1405–1420 (2020)
  * (40) Witherden, F.D., Vincent, P.E.: On the identification of symmetric quadrature rules for finite element methods.  Computers & Mathematics with Applications 69(10), 1232–1241 (2015)
  * (41) Yamakawa, S., Gentilini, I., Shimada, K.: Subdivision templates for converting a non-conformal hex-dominant mesh to a conformal hex-dominant mesh without pyramid elements.  Engineering with Computers 27(1), 51–65 (2011)
  * (42) Yamakawa, S., Shimada, K.: Increasing the number and volume of hexahedral and prism elements in a hex-dominant mesh by topological transformations.  In: Proceedings of the 12th International Meshing Roundtable, pp. 403–413 (2003)
  * (43) Yamakawa, S., Shimada, K.: Converting a tetrahedral mesh to a prism–tetrahedral hybrid mesh for FEM accuracy and efficiency.  International Journal for Numerical Methods in Engineering 80(1), 74–102 (2009)
  * (44) Yamakawa, S., Shimada, K.: 88-Element solution to Schneiders’ pyramid hex-meshing problem.  International Journal for Numerical Methods in Biomedical Engineering 26(12), 1700–1712 (2010)
  * (45) Yamakawa, S., Shimada, K.: Automatic all-hex mesh generation of thin-walled solids via a conformal pyramid-less hex, prism, and tet mixed mesh.  In: Proceedings of the 20th International Meshing Roundtable, pp. 125–141. Springer (2011)
  * (46) Zamboj, M.: Sections and shadows of four-dimensional objects.  Nexus Network Journal 20(2), 475–487 (2018)