# Design and performance of a space–time virtual element method for the heat
equation on prismatic meshes ††thanks: The authors have been funded by the
Austrian Science Fund (FWF) through the projects F 65 and P 33477 (I.
Perugia), and the Italian Ministry of University and Research through the PRIN
project “NA-FROM-PDEs” and PNRR-M4C2-I1.4-NC-HPC-Spoke6 (S. Gómez). S. Gómez
acknowledges the kind hospitality of the Erwin Schrödinger International
Institute for Mathematics and Physics (ESI), where part of this research was
developed.

Sergio Gómez , Lorenzo Mascotto , Ilaria Perugia55footnotemark: 5 Department
of Mathematics, University of Pavia, 27100 Pavia, Italy
(sergio.gomez01@universitadipavia.it)Erwin Schrödinger Institute for
Mathematics and Physics, University of Vienna, AustriaDepartment of
Mathematics and Applications, University of Milano-Bicocca, 20125 Milan, Italy
(lorenzo.mascotto@unimib.it)Faculty of Mathematics, University of Vienna, 1090
Vienna, Austria (ilaria.perugia@univie.ac.at)IMATI-CNR, Pavia, Italy

###### Abstract

We present a space–time virtual element method for the discretization of the
heat equation, which is defined on general prismatic meshes and variable
degrees of accuracy. Strategies to handle efficiently the space–time mesh
structure are discussed. We perform convergence tests for the $h$\- and
$hp$-versions of the method in case of smooth and singular solutions, and test
space–time adaptive mesh refinements driven by a residual-type error
indicator.

AMS subject classification: 35K05; 65N12; 65N30.

Keywords: virtual element methods; heat equation; space–time methods;
polytopic meshes.

##  1 Introduction

Space–time Galerkin methods aim at approximating solutions to time-dependent
partial differential equations treating the time variable as an additional
space variable. Even though the foundation of space–time Galerkin methods
traces back to the 70ies of the last century [15] and a few contributions were
developed in the twenty years to follow [2, 3, 8, 14, 9], only in the last two
decades there has been a growing attention on this topic, mainly due to the
improved performance of computers.

Compared to the time-stepping approach, space–time Galerkin methods have some
important upsides and features: the discrete solution can be evaluated on the
whole space–time domain and not only at a finite number of times without
additional post-processing of the discrete solution; such methods allow for
space–time adaptivity; it is possible to design space–time parallel solvers.

In this paper, we focus on the approximation of solutions to the heat
equation. Several space–time methods have been designed to this aim and they
can be classified into two main groups. The first one is based on the
discretization of a standard Petrov-Galerkin formulation [6]; see [1, 21] for
continuous finite element methods, [20] for a wavelet method, [26, 4] for
discontinuous Galerkin methods, [16] for an isogeometric method, [24] for a
coercive method based on a Hilbert transformation of the test space, and [25]
for a mixed finite element method. Residual-type error indicators for the
method of [21] were considered in [22, 23]. Possible drawbacks of employing
continuous finite elements are that suboptimal convergence rates are obtained
for some singular solutions, and incompatible boundary and initial conditions
cannot be naturally handled.

The second group is based on first order system least squares discretizations
(FOSLS); see [11, 12, 27] for finite element methods and [17] for isogeometric
methods. Space–time FOSLS finite elements naturally provide reliable and
efficient error indicators; see [11, 18, 12]. However, they require the
computation of an additional vector-valued flux variable.

In this paper, we extend the nonconforming space–time virtual element method
of [13] to the case of general prismatic meshes and nonuniform degrees of
accuracy. This method allows for the use of space–time meshes with hanging
facets (nodes, edges, faces) and is based on discontinuous in time test and
trial functions that are solutions to local space–time problems with
polynomial data. This approach provides a natural framework for space–time
adaptivity without need of re-meshing neighbouring elements. If the space–time
mesh is decomposed into separate time-slabs, the global linear system can be
split into much smaller systems that can be solved sequentially. Due to the
nonconformity across time-like facets, the design and implementation of the
method are independent of the spatial dimension. No artificial compatibility
of initial and boundary conditions is enforced on the trial virtual element
space, and optimal convergence rates have been achieved also for singular
solutions. However, the stability analysis of the method relies on a discrete
inf-sup condition and the method requires a stabilization term for the spatial
Laplacian that needs to be carefully designed; see [13].

The main goals of this manuscript are:

  * •

to construct a space–time virtual element method on general prismatic meshes,
possibly with hanging facets and variable degrees of accuracy;

  * •

to discuss suitable strategies to handle efficiently the space–time mesh
structure;

  * •

to investigate $h$\- and $hp$\- refinements numerically;

  * •

to test space–time adaptive mesh refinements driven by a residual-type error
indicator.

In the remainder of this section, we introduce some notation and the model
problem. Finally, we outline the structure of the paper.

###  1.1 Notation

We denote the first partial derivative with respect to the time variable $t$
by $\partial_{t}{}$, and the spatial gradient and Laplacian operators by
$\nabla_{\mathbf{x}}{},\ \Delta_{\mathbf{x}}{}$, respectively.

Standard notation for Sobolev spaces is employed. For a given domain $D$ in
$\mathbb{R}^{d}$, $d\in\mathbb{N}$, $H^{s}(D)$ represents the standard Sobolev
space of order $s$ in $\mathbb{N}$ endowed with the standard inner product
$(\cdot,\cdot)_{s,D}$, the seminorm ${\left|{\cdot}\right|}_{s,D}$, and the
norm ${\left\|{\cdot}\right\|}_{s,D}$. In particular, we let $H^{0}(D)$ be the
space $L^{2}(D)$ of Lebesgue square integrable functions over $D$ and
$H_{0}^{1}(D)$ be the subspace of functions in $H^{1}(D)$ with zero trace on
$\partial D$. If $s$ is a fractional or negative number, then the Sobolev
space $H^{s}(D)$ is defined by means of interpolation and duality,
respectively. We denote the duality product between $H^{-1}(D)$ and
$H^{1}_{0}(D)$ by $\langle\cdot,\cdot\rangle$. The Sobolev spaces on $\partial
D$ are defined analogously and denoted by $H^{s}(\partial D)$.

Given $s$ in $\mathbb{R}$, a time interval $(a,b)$, and a Banach space
$(Z,{\left\|{\cdot}\right\|}_{Z})$, we introduce the Bochner space
$H^{s}(a,b;Z)$. In particular, we define

$Y:=L^{2}(0,T;H_{0}^{1}(\Omega)),\quad X:=\left\\{v\in Y\cap
H^{1}(0,T;H^{-1}(\Omega))\mid
v(\mathbf{x},0)=0\;\forall\mathbf{x}\in\Omega\right\\}.$ (1)

We endow $Y$ and $X$ with the norms

${\left\|{v}\right\|}_{Y}^{2}:=\nu\int_{0}^{T}{\left|{v(\cdot,t)}\right|}_{1,\Omega}^{2}\mbox{dt},\qquad\qquad{\left\|{v}\right\|}_{X}^{2}:=c_{H}{\left\|{\partial_{t}{v}}\right\|}_{L^{2}(0,T;H^{-1}(\Omega))}^{2}+{\left\|{v}\right\|}_{Y}^{2},$
(2)

respectively, where $\nu$ and $c_{H}$ are the problem coefficients defined at
the beginning of Section 1.2 below, and we have set

${\left\|{\phi}\right\|}_{L^{2}(0,T;H^{-1}(\Omega))}:=\sup_{0\neq v\in
Y}\frac{\int_{0}^{T}\langle\phi,v\rangle\mbox{dt}}{{\left\|{v}\right\|}_{Y}}.$

We denote the space of polynomials in $d$ variables of degree at most $p$ on a
domain $D\subset\mathbb{R}^{d}$ by $\mathbb{P}_{p}(D)$.

###  1.2 Model problem

We consider the heat equation on the space–time domain
$Q_{T}:=\Omega\times(0,T)$, where $\Omega\subset\mathbb{R}^{d}$, $d=1,\ 2,\
3$, and $T>0$ are the (bounded) spatial domain and the final time,
respectively.

Associated with $Q_{T}$, we introduce the space-like surface
$\Sigma_{0}:=\Omega\times\left\\{0\right\\}$ and the time-like surface
$\Sigma_{D}:=\partial\Omega\times(0,T)$.

Let $f:Q_{T}\to\mathbb{R}$, $u_{0}:\Sigma_{0}\to\mathbb{R}$, and
$g:\Sigma_{D}\to\mathbb{R}$ denote the source term, initial condition, and
Dirichlet boundary condition, respectively. We assume $f\in
L^{2}(0,T;H^{-1}(\Omega))$. We further consider $c_{H}>0$ and $\nu>0$ given
positive constant volumetric heat capacity and thermal conductivity,
respectively.

The heat equation in strong formulation reads: find a function
$u:Q_{T}\to\mathbb{R}$ such that

$\begin{cases}c_{H}\partial_{t}{u}-\nu\Delta_{\mathbf{x}}u=f&\text{in
}Q_{T},\\\ u=u_{0}&\text{on }\Sigma_{0},\\\ u=g&\text{on
}\Sigma_{D}.\end{cases}$ (3)

To simplify the presentation, we henceforth assume that $u_{0}=0$ and $g=0$.
Nonetheless, in Sections 4 and 5, we present numerical results for solutions
with inhomogeneous initial and boundary conditions. Inhomogeneous boundary
conditions can be treated as in the elliptic case and inhomogeneous initial
conditions can be dealt with by standard lifting arguments.

Define the space–time bilinear form $b(\cdot,\cdot):X\times Y\to\mathbb{R}$ as

$b(u,v):=\int_{0}^{T}\langle
c_{H}\partial_{t}{u},v\rangle\mbox{dt}+\int_{0}^{T}\int_{\Omega}\nu\nabla_{\mathbf{x}}u\cdot\nabla_{\mathbf{x}}v\,\mbox{d}\mathbf{x}\mbox{dt}.$
(4)

We consider the Petrov-Galerkin weak formulation of (3) as in [6]:

$\begin{cases}\text{find }u\in X\text{ such that}\\\
b(u,v)=\int_{0}^{T}\langle f,v\rangle\mbox{dt}\qquad\forall v\in
Y.\end{cases}$ (5)

Problem (5) is well posed; see, e.g., [21, Corollary $2.3$].

#### Structure of the paper.

In Section 2, we design a space–time virtual element method for the
approximation of solutions to (5) on general prismatic meshes; in Section 3,
strategies to handle efficiently the space–time mesh structure are discussed,
including flagging strategies that might be used for other polytopic methods
as well; in Section 4, we assess the convergence of the $h$\- and
$hp$-versions of the method on some singular solutions; in Section 5, a
computable residual-type error indicator is defined and used to lead adaptive
mesh refinements.

##  2 The space–time virtual element method

We describe a nonconforming space–time virtual element method for the
approximation of solutions to (5). We proceed in several steps: we introduce
general prismatic meshes in Section 2.1; design local virtual element spaces
and describe their degrees of freedom (DoFs) in Section 2.2; show that the
choice of the DoFs allows for the computation of orthogonal projectors onto
polynomial spaces in Section 2.3; define global nonconforming space–time
virtual element spaces in Section 2.4; detail the discrete bilinear forms in
Section 2.5; present the method in Section 2.6.

###  2.1 General prismatic meshes

We consider sequences of meshes $\\{\mathcal{T}_{h}\\}$ consisting of
nonoverlapping prismatic-type polytopes covering $Q_{T}$ in the following
sense: each element $\mathcal{K}$ in $\mathcal{T}_{h}$ can be written as
$\mathcal{K}_{\mathbf{x}}\times\mathcal{K}_{t}$ for some open $d$-dimensional
polytope $\mathcal{K}_{\mathbf{x}}$ contained in $\Omega$ with boundary
$\partial\mathcal{K}_{\mathbf{x}}$, and some time subinterval
$\mathcal{K}_{t}:=(a_{t},b_{t})$ of $(0,T)$. We denote the diameter of
$\mathcal{K}_{\mathbf{x}}$ and the length of $\mathcal{K}_{t}$ by
$h_{\mathcal{K}_{\mathbf{x}}}$ and $h_{\mathcal{K}_{t}}$, respectively.

We call “mesh facet” any intersection
$\partial\mathcal{K}^{1}\cap\partial\mathcal{K}^{2}$,
$\partial\mathcal{K}^{1}\cap(\partial\Omega\times(0,T))$,
$\partial\mathcal{K}^{1}\cap(\Omega\times\\{0\\})$, or
$\partial\mathcal{K}^{1}\cap(\Omega\times\\{T\\})$, for given
$\mathcal{K}^{1},\mathcal{K}^{2}$ in $\mathcal{T}_{h}$, that has positive
$d$-dimensional measure and is contained in a $d$-dimensional hyperplane. For
each element $\mathcal{K}$ in $\mathcal{T}_{h}$, we identify two types of
nonoverlapping mesh facets: space-like facets $K_{\mathbf{x}}$, whose union
gives $\mathcal{K}_{\mathbf{x}}\times\\{a_{t}\\}$, and time-like facets
$F:=F_{\mathbf{x}}\times F_{t}$, where
$F_{\mathbf{x}}\subset\partial\mathcal{K}_{\mathbf{x}}$ is a facet of
$\mathcal{K}_{\mathbf{x}}$ and $F_{t}\subset\mathcal{K}_{t}$. 111For our
purposes, we do not need to classify the facets of $\mathcal{K}$ on
$\mathcal{K}_{\mathbf{x}}\times\\{b_{t}\\}$. One may consider that the facet
$\mathcal{K}_{\mathbf{x}}\times\\{b_{t}\\}$ is a unique facet of
$\mathcal{K}$. We collect the space-like and time-like facets of $\mathcal{K}$
into the sets $\mathcal{F}_{\mathcal{K}}^{\text{space}}$ and
$\mathcal{F}_{\mathcal{K}}^{\text{time}}$.

Associated with each element $\mathcal{K}$ in $\mathcal{T}_{h}$, we define $K$
as the $(d+1)$-dimensional closed polytope $K$ whose interior is $\mathcal{K}$
and whose boundary facets are those in the set
$\mathcal{F}_{\mathcal{K}}^{\text{space}}\cup\mathcal{F}_{\mathcal{K}}^{\text{time}}\cup\\{\mathcal{K}_{\mathbf{x}}\times\\{b_{t}\\}\\}$.

For any time-like facet $F=F_{\mathbf{x}}\times F_{t}$ in
$\mathcal{F}_{\mathcal{K}}^{\text{time}}$, we define

$h_{F_{\mathbf{x}}}:=\begin{cases}\min\\{h_{\mathcal{K}_{\mathbf{x}}},h_{\widetilde{\mathcal{K}}_{\mathbf{x}}}\\}&\text{if
}F=\mathcal{K}\cap\widetilde{\mathcal{K}}\text{ for some
}\widetilde{\mathcal{K}}=\widetilde{\mathcal{K}}_{\mathbf{x}}\times\widetilde{\mathcal{K}}_{t}\in\mathcal{T}_{h},\\\
h_{K_{\mathbf{x}}}&\text{if }F_{\mathbf{x}}\subset\partial\Omega.\end{cases}$

We define the sets of _space-like_ and _time-like_ facets of $\mathcal{T}_{h}$
as

$\begin{split}\mathcal{F}_{h}^{\text{space}}:=\bigcup\limits_{K\in\mathcal{T}_{h}}\mathcal{F}_{\mathcal{K}}^{\text{space}},\qquad\qquad\mathcal{F}_{h}^{\text{time}}:=\bigcup\limits_{K\in\mathcal{T}_{h}}\mathcal{F}_{\mathcal{K}}^{\text{time}}.\end{split}$

Hanging nodes ($1+1$ dimensional case), edges ($2+1$ dimensional case), and
faces ($3+1$ dimensional case) are included within this structure of the mesh.
In the absence of hanging facets, i.e., for tensor-product-in-time meshes, it
is not necessary to distinguish between $\mathcal{K}$ and $K$, and all the
definitions in this work reduce to those in [13]. For each $F$ in
$\mathcal{F}_{h}^{\text{time}}$, we fix a normal unit vector $\mathbf{n}_{F}$
in $\mathbb{R}^{d+1}$ of the form $(\mathbf{n}_{F_{\mathbf{x}}},0)$,
$\mathbf{n}_{F_{\mathbf{x}}}\in\mathbb{R}^{d}$. In Figure 1, we illustrate the
definitions in this section.

$\mathcal{K}$$\widetilde{\mathcal{K}}^{1}$$\widetilde{\mathcal{K}}^{2}$

$\mathcal{K}$$\mathcal{K}_{\mathbf{x}}$$\mathcal{K}_{t}$

$K$$K_{\mathbf{x}}^{1}$$K_{\mathbf{x}}^{2}$$F^{2}$$F^{1}$$F^{3}$

Figure 1: Example of a prismatic space–time mesh. Left panel: The prismatic
partition $\mathcal{T}_{h}$ of the space–time domain. Central panel: Zoom of
the element $\mathcal{K}\in\mathcal{T}_{h}$. Right panel: The associated
closed hexagon $K$ with two space-like facets $K_{\mathbf{x}}^{1}$ and
$K_{\mathbf{x}}^{2}$, and three time-like facets $F^{1},F^{2}$, and $F^{3}$.
By definition,
$h_{F_{\mathbf{x}}^{1}}=\min\\{h_{\mathcal{K}_{\mathbf{x}}},h_{\widetilde{\mathcal{K}}^{1}_{\mathbf{x}}}\\}$,
$h_{F_{\mathbf{x}}^{2}}=\min\\{h_{\mathcal{K}_{\mathbf{x}}},h_{\widetilde{\mathcal{K}}^{2}_{\mathbf{x}}}\\}$,
and $h_{F_{\mathbf{x}}^{3}}=h_{\mathcal{K}_{\mathbf{x}}}$.

For a given mesh $\mathcal{T}_{h}$, we define the broken Sobolev space of
order $s\in\mathbb{R}^{+}$ as

$H^{s}(\mathcal{T}_{h}):=\left\\{v\in L^{2}(\Omega)\ \middle|\ v{}_{|K}\in
H^{s}(K)\;\forall K\in\mathcal{T}_{h}\right\\},$

and endow it with the standard broken norm
${\left\|{\cdot}\right\|}_{s,\mathcal{T}_{h}}$ and seminorm
${\left|{\cdot}\right|}_{s,\mathcal{T}_{h}}$. We denote the space of piecewise
polynomials of maximum degree $p$ in $\mathbb{N}$ over $\mathcal{T}_{h}$ by
$\mathcal{S}_{p}(\mathcal{T}_{h})$.

###  2.2 Local space–time virtual element spaces

We introduce local space–time virtual element spaces, extending the
construction in [13] to more general meshes as in Section 2.1.

Introduce the scaling factors

$\widetilde{c}_{H}^{\mathcal{K}}:=h_{\mathcal{K}_{t}},\qquad\qquad\widetilde{\nu}^{\mathcal{K}}:=h_{\mathcal{K}_{\mathbf{x}}}^{2}.$
(6)

Given $p$ in $\mathbb{N}$ and $\mathcal{K}$ an element in $\mathcal{T}_{h}$,
we define the local space–time virtual element space as

$\begin{split}V_{h}(K):=\Big{\\{}v\in L^{2}(\mathcal{K})\mid\
&\widetilde{c}_{H}^{\mathcal{K}}\partial_{t}{v_{h}}-\widetilde{\nu}^{\mathcal{K}}\Delta_{\mathbf{x}}v_{h}\in\mathbb{P}_{p-1}\left(\mathcal{K}\right);\\\\[1.99997pt]
&v_{h}{}_{|\mathcal{K}_{\mathbf{x}}}\in\mathbb{P}_{p}\left(\mathcal{K}_{\mathbf{x}}\right);\;\mathbf{n}_{F_{\mathbf{x}}}\cdot\nabla_{\mathbf{x}}v_{h}{}_{|F}\in\mathbb{P}_{p}\left(F\right)\,\forall
F\in\mathcal{F}_{\mathcal{K}}^{\text{time}}\Big{\\}}.\end{split}$ (7)

Functions in $V_{h}(K)$ are not known in closed form. Yet,
$\mathbb{P}_{p}(\mathcal{K})$ is contained in $V_{h}(K)$.

Given an element $\mathcal{K}$ of $\mathcal{T}_{h}$ and any of its time-like
facets $F$, let
$\\{m_{\alpha}^{\mathcal{K}}\\}_{\alpha=1}^{\dim(\mathbb{P}_{p-1}(\mathcal{K}))}$,
$\\{m_{\beta}^{F}{}\\}_{\beta=1}^{\dim(\mathbb{P}_{p}(F))}$, and
$\\{m_{\gamma}^{\mathcal{K}_{\mathbf{x}}}\\}_{\gamma=1}^{\dim(\mathbb{P}_{p}(\mathcal{K}_{\mathbf{x}}))}$
be given bases of $\mathbb{P}_{p-1}(\mathcal{K})$, $\mathbb{P}_{p}(F)$, and
$\mathbb{P}_{p}(\mathcal{K}_{\mathbf{x}})$, respectively. Those basis elements
are assumed to be invariant with respect to translations and dilations.

Introduce the following set of linear functionals:

  * •

the _bulk_ moments

$\frac{1}{|\mathcal{K}|}\int_{K}v_{h}\
m_{\alpha}^{\mathcal{K}}\mbox{d}\mathbf{x}\
\mbox{dt},\qquad\alpha=1,\dots,\dim(\mathbb{P}_{p-1}(\mathcal{K}));$ (8)

  * •

for all _time-like_ facets $F\in\mathcal{F}_{\mathcal{K}}^{\text{time}}$, the
_time-like_ moments

$\frac{1}{|F|}\int_{F}v_{h}\ m_{\beta}^{F}\mbox{d}S\
\mbox{dt},\qquad\beta=1,\dots,\dim(\mathbb{P}_{p}(F));$ (9)

  * •

the _space-like_ moments

$\frac{1}{|\mathcal{K}_{\mathbf{x}}|}\int_{\mathcal{K}_{\mathbf{x}}}v_{h}(\cdot,t_{n-1})m_{\gamma}^{\mathcal{K}_{\mathbf{x}}}\mbox{d}\mathbf{x},\qquad\gamma=1,\dots,\dim(\mathbb{P}_{p}(\mathcal{K}_{\mathbf{x}})).$
(10)

These linear functionals constitute a set of unisolvent degrees of freedom for
the space $V_{h}(K)$. To see this, it suffices to extend [13, Lemma $2.1$] to
the case of the general prismatic-type elements introduced in Section 2.1.

###  2.3 Polynomial projections

We define four orthogonal projectors onto polynomial spaces as in [13]: for
all elements $\mathcal{K}$ in $\mathcal{T}_{h}$ and all $\varepsilon>0$,

  1. 1.

we introduce
$\Pi^{N}_{p}:H^{\frac{1}{2}+\varepsilon}(\mathcal{K}_{t};L^{2}(\mathcal{K}_{\mathbf{x}}))\cap
L^{2}(\mathcal{K}_{t};H^{1}(\mathcal{K}_{\mathbf{x}}))\to\mathbb{P}_{p}(\mathcal{K})$
as

$\displaystyle\int_{\mathcal{K}_{t}}\int_{\mathcal{K}_{\mathbf{x}}}\nabla_{\mathbf{x}}q_{p}^{\mathcal{K}}\cdot\nabla_{\mathbf{x}}\left(\Pi^{N}_{p}{v}-v\right)\
\mbox{d}\mathbf{x}\mbox{dt}=0\quad\forall
q_{p}^{\mathcal{K}}\in\mathbb{P}_{p}\left(\mathcal{K}\right)\backslash\mathbb{P}_{p}\left(\mathcal{K}_{t}\right),$
(11a)
$\displaystyle\int_{\mathcal{K}_{t}}\int_{\mathcal{K}_{\mathbf{x}}}q_{p-1}(t)\left(\Pi^{N}_{p}{v}-v\right)\mbox{d}\mathbf{x}\mbox{dt}=0\quad\forall
q_{p-1}\in\mathbb{P}_{p-1}\left(\mathcal{K}_{t}\right),$ (11b)
$\displaystyle\int_{\mathcal{K}_{\mathbf{x}}}\left(\Pi^{N}_{p}{v}(\mathbf{x},t_{n-1})-v(\mathbf{x},t_{n-1})\right)\mbox{d}\mathbf{x}=0;$
(11c)

  2. 2.

we introduce
$\Pi^{\star}_{p}:H^{\frac{1}{2}+\varepsilon}(\mathcal{K}_{t};L^{2}(\mathcal{K}_{\mathbf{x}}))\to\mathbb{P}_{p}(\mathcal{K})$
as

$\displaystyle\int_{\mathcal{K}_{t}}\int_{\mathcal{K}_{\mathbf{x}}}q_{p-1}^{\mathcal{K}}\left(\Pi^{\star}_{p}{v}-v\right)\mbox{d}\mathbf{x}\mbox{dt}=0\quad\forall
q_{p-1}^{\mathcal{K}}\in\mathbb{P}_{p-1}\left(\mathcal{K}\right),$ (12a)
$\displaystyle\int_{\mathcal{K}_{\mathbf{x}}}q_{p}^{\mathcal{K}_{\mathbf{x}}}\left(\Pi^{\star}_{p}{v}(\mathbf{x},t_{n-1})-v(\mathbf{x},t_{n-1})\right)\mbox{d}\mathbf{x}=0\quad\forall
q_{p}^{\mathcal{K}_{\mathbf{x}}}\in\mathbb{P}_{p}\left(\mathcal{K}_{\mathbf{x}}\right);$
(12b)

  3. 3.

we introduce
$\Pi^{0,\mathcal{K}}_{p-1}:L^{2}(\mathcal{K})\to\mathbb{P}_{p-1}(\mathcal{K})$
as

$(q_{p-1}^{\mathcal{K}},v-\Pi^{0,\mathcal{K}}_{p-1}v)_{0,\mathcal{K}}=0\qquad\forall
q_{p-1}^{\mathcal{K}}\in\mathbb{P}_{p-1}(\mathcal{K});$ (13)

  4. 4.

for all _time-like_ facets $F$ in $\mathcal{F}_{\mathcal{K}}^{\text{time}}$,
we introduce $\Pi^{0,F}_{p}:L^{2}(F)\to\mathbb{P}_{p}(F)$ as

$(q_{p}^{F},v-\Pi^{0,F}_{p}v)_{0,F}=0\qquad\forall
q_{p}^{F}\in\mathbb{P}_{p}(F).$ (14)

Let $v$ be a function in $V_{h}(K)$ with given degrees of freedom (8)–(10).
Then, the proof of the well posedness and computability of the above four
projectors follows by extending the results in [13, Section $2.2$] to the case
of the general prismatic-type elements in Section 2.1.

###  2.4 Global nonconforming space–time virtual element spaces

We design a global virtual element space $Y_{h}$, consisting of functions that
are discontinuous in time and nonconforming in space.

To this aim, we introduce the jump operator on each time-like facet $F$ as the
functional
$\left[\\!\left[\cdot\right]\\!\right]_{F}:H^{\frac{1}{2}+\varepsilon}(\mathcal{T}_{h})\to[L^{2}(F)]^{d+1}$,
$\varepsilon>0$, given by

$\left[\\!\left[v\right]\\!\right]_{F}:=\begin{cases}v_{|\mathcal{K}_{1}}\mathbf{n}_{K_{1}}^{F}+v_{|\mathcal{K}_{2}}\mathbf{n}_{K_{2}}^{F}&\text{if
}F\subset\partial\mathcal{K}_{1}\cap\partial\mathcal{K}_{2}\text{ is an
internal facet, for }\mathcal{K}_{1},\mathcal{K}_{2}\in\mathcal{T}_{h}\\\
v_{|\mathcal{K}_{3}}\mathbf{n}_{K_{3}}^{F}&\text{if
}F\subset\partial\mathcal{K}_{3}\text{ is a boundary facet, for
}\mathcal{K}_{3}\in\mathcal{T}_{h}.\end{cases}$

We define the Sobolev nonconforming space of order $p$ associated with the
mesh $\mathcal{T}_{h}$ as

$\begin{split}\mathcal{NC}_{p}(\mathcal{T}_{h}):=\Big{\\{}v\in
L^{2}(Q_{T})\Big{|}&\,v_{|\mathcal{K}}\in
L^{2}(\mathcal{K}_{t},H^{1}(\mathcal{K}_{\mathbf{x}}))\;\forall\mathcal{K}\in\mathcal{T}_{h};\\\
&\int_{F}q_{p}^{F}\left[\\!\left[v\right]\\!\right]_{F}\cdot\mathbf{n}_{F}\mbox{d}S\
=0\quad\forall q_{p}^{F}\in\mathbb{P}_{p}\left(F\right),\;\forall
F\in\mathcal{F}_{h}^{\text{time}}\Big{\\}}.\end{split}$ (15)

We define the global virtual element space as follows:

$\begin{split}Y_{h}:=\Big{\\{}v_{h}\in L^{2}(Q_{T})\ \Big{|}\
&v_{h}{}_{|\mathcal{K}}\in V_{h}(K)\;\;\forall K\in\mathcal{T}_{h};\
v_{h}\in\mathcal{NC}_{p}(\mathcal{T}_{h})\Big{\\}}.\end{split}$

###  2.5 The discrete bilinear forms

Since functions in the space–time virtual element space $Y_{h}$ are not
available in closed form, we discretize the bilinear forms by computable
counterparts. On each element $\mathcal{K}$, define the local continuous
bilinear form and the induced norm

$a^{\mathcal{K}}(u_{h},v_{h}):=\nu(\nabla_{\mathbf{x}}u_{h},\nabla_{\mathbf{x}}v_{h})_{0,\mathcal{K}},\qquad\qquad{\left|{v_{h}}\right|}_{Y(\mathcal{K})}^{2}:=a^{\mathcal{K}}(v_{h},v_{h})=\nu{\left\|{\nabla_{\mathbf{x}}v_{h}}\right\|}_{0,\mathcal{K}}^{2}.$

Let $S^{K}:V_{h}(K)\times V_{h}(K)\to\mathbb{R}$ be any symmetric bilinear
form, which is computable via the degrees of freedom and satisfies the
following property: there exist $0<c_{*}<c^{*}$ independent of $K$ such that

$c_{*}{\left|{v_{h}}\right|}_{Y(\mathcal{K})}^{2}\leq\nu
S^{K}(v_{h},v_{h})\leq
c^{*}{\left|{v_{h}}\right|}_{Y(\mathcal{K})}^{2}\qquad\qquad\forall v_{h}\in
V_{h}(K)\cap\ker(\Pi^{N}_{p}).$ (16)

Define

$a_{h}^{K}(u_{h},v_{h}):=a^{\mathcal{K}}(\Pi^{N}_{p}u_{h},\Pi^{N}_{p}v_{h})+\nu
S^{K}((I-\Pi^{N}_{p})u_{h},(I-\Pi^{N}_{p})v_{h}).$ (17)

Extending [13, Lemma 2.7] to the case of general prismatic-type elements as in
Section 2.1, it is possible to prove that there exist
$0<\alpha_{*}<\alpha^{*}$ independent of $K$ such that the following local
stability bounds are valid:

$\alpha_{*}{\left|{v_{h}}\right|}_{Y(\mathcal{K})}^{2}\leq
a_{h}^{K}(v_{h},v_{h})\leq\alpha^{*}{\left|{v_{h}}\right|}_{Y(\mathcal{K})}^{2}\qquad\qquad\forall
v_{h}\in V_{h}(K).$ (18)

Besides, we introduce upwind-type terms, which allows for imposing weakly the
continuity in time of the trial functions: for all space-like facets
$K_{\mathbf{x}}\subset\Omega\times\\{t^{*}\\}$, $t^{*}\in[0,T]$,

$\mathcal{U}^{K_{\mathbf{x}}}(u_{h}):=\begin{cases}c_{H}\Pi^{\star}_{p}u_{h}{}_{|\mathcal{K}}(\cdot,t^{*})&\text{if
}t^{*}=0\\\
c_{H}\left(\Pi^{\star}_{p}u_{h}{}_{|\mathcal{K}^{+}}(\cdot,t^{*})-\Pi^{\star}_{p}u_{h}{}_{|\mathcal{K}^{-}}(\cdot,t^{*})\right)&\text{if
}t^{*}>0,\end{cases}$ (19)

where $K_{\mathbf{x}}\in\mathcal{F}_{\mathcal{K}}^{\text{space}}$ if
$t^{*}=0$, and $K_{\mathbf{x}}\in F^{\text{space}}_{\mathcal{K}^{+}}$ is a
subset of $\partial\mathcal{K}^{+}\cap\partial\mathcal{K}^{-}$ if $t^{*}>0$.

######  Remark 1.

The case of nonzero initial conditions can be dealt with by modifying the
definition of the upwind functional at time $t=0$ as follows: for all space-
like facets $K_{\mathbf{x}}\subset\Omega\times\\{0\\}$,

$\mathcal{U}^{K_{\mathbf{x}}}(u_{h}):=c_{H}\left(\Pi^{\star}_{p}u_{h}{}_{|\mathcal{K}}(\cdot,0)-u_{0}\right).$

The computability of the upwind terms follows from the definition of the
space-like moment degrees of freedom (10) and the computability of the
projector $\Pi^{\star}_{p}$.

We define a discrete counterpart of the bilinear form $b(\cdot,\cdot)$ in (4)
as follows:

$b_{h}(u_{h},v_{h}):=\sum_{\mathcal{K}\in\mathcal{T}_{h}}\left[c_{H}(\partial_{t}{\Pi^{\star}_{p}}u_{h},v_{h})_{0,\mathcal{K}}+a_{h}^{K}(u_{h},v_{h})\right]+\sum_{K_{\mathbf{x}}\in\mathcal{F}_{h}^{\text{space}}}\left(\mathcal{U}^{K_{\mathbf{x}}}(u_{h}),v_{h}{}_{|\mathcal{K}^{+}}\right)_{0,K_{\mathbf{x}}}.$
(20)

A computable stabilization satisfying (18) is given by

$\begin{split}S^{K}(u_{h},v_{h})&:=p^{2}h_{\mathcal{K}_{\mathbf{x}}}^{-2}(\Pi^{0,\mathcal{K}}_{p-1}u_{h},\Pi^{0,\mathcal{K}}_{p-1}v_{h})_{0,\mathcal{K}}+\sum_{F\in\mathcal{F}_{\mathcal{K}}^{\text{time}}}ph_{F_{\mathbf{x}}}^{-1}(\Pi^{0,F}_{p}u_{h},\Pi^{0,F}_{p}v_{h})_{0,F}\\\
&\quad+ph_{\mathcal{K}_{t}}h_{\mathcal{K}_{\mathbf{x}}}^{-2}(u_{h},v_{h})_{0,\mathcal{K}_{\mathbf{x}}}.\end{split}$
(21)

###  2.6 The method

Henceforth, we assume that $f$ belongs to $L^{2}(Q_{T})$. The space–time
virtual element method for the approximation of solutions to (5) reads

$\begin{cases}\text{find }u_{h}\in Y_{h}\text{ such that}\\\
b_{h}(u_{h},v_{h})=(f,\Pi^{0}_{p-1}v_{h})_{0,\Omega}\qquad\forall v_{h}\in
Y_{h}.\end{cases}$ (22)

For Cartesian-product-in-time meshes, the well posedness of the method and _a
priori_ error estimates were proven in [13, Sections 3 and 4] based on a
discrete inf-sup argument. That analysis can be extended to the meshes
introduced in Section 2.1, using the norms

$\displaystyle{\left\|{v}\right\|}_{Y(\mathcal{T}_{h})}^{2}$
$\displaystyle:=\sum_{\mathcal{K}\in\mathcal{T}_{h}}{\left\|{v}\right\|}_{Y(\mathcal{K})}^{2},$
$\displaystyle{\left\|{v}\right\|}_{X(\mathcal{T}_{h})}^{2}$
$\displaystyle:={\left\|{v}\right\|}_{Y(\mathcal{T}_{h})}^{2}+{\left\|{\mathfrak{N}_{h}\Pi^{\star}_{p}v}\right\|}_{Y(\mathcal{T}_{h})}^{2}$
$\displaystyle\quad+\frac{c_{H}}{2}\Big{(}{\left\|{\Pi^{\star}_{p}v(\cdot,0)}\right\|}_{0,\Omega}^{2}+\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\scriptsize\begin{tabular}[]{c}$K_{\mathbf{x}}\in\mathcal{F}_{\mathcal{K}}^{\text{space}}$\\\
$K_{\mathbf{x}}\not\subset\Omega\times\\{0\\}$\end{tabular}}{\left\|{\mathcal{U}^{K_{\mathbf{x}}}(v)}\right\|}_{0,K_{\mathbf{x}}}^{2}+{\left\|{\Pi^{\star}_{p}v(\cdot,T)}\right\|}_{0,\Omega}^{2}\Big{)},$

where the discrete Newton potential
$\mathfrak{N}_{h}:\mathcal{S}_{p}(\mathcal{T}_{h})\rightarrow Y_{h}$ is
defined as follows: for all $v_{h}\in Y_{h}$,

$\begin{split}a_{h}(\mathfrak{N}_{h}\phi_{h},v_{h})=c_{H}\Big{(}(\partial_{t}{\phi_{h}},v_{h})_{0,Q_{T}}+\left(\phi_{h}(\cdot,0),v_{h}(\cdot,0)\right)_{0,\Omega}+\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\scriptsize\begin{tabular}[]{c}$K_{\mathbf{x}}\in\mathcal{F}_{\mathcal{K}}^{\text{space}}$\\\
$K_{\mathbf{x}}\not\subset\Omega\times\\{0\\}$\end{tabular}}\Big{(}\mathcal{U}^{K_{\mathbf{x}}}(\phi_{h}),v_{h}{}_{|_{\mathcal{K}^{+}}}\Big{)}_{0,K_{\mathbf{x}}}\Big{)}.\end{split}$

In particular, we have the following result.

######  Theorem 2.1.

There exists $\gamma_{I}>0$ independent of $\mathcal{T}_{h}$ such that

$\sup_{0\neq v_{h}\in
Y_{h}}\frac{b_{h}(u_{h},v_{h})}{{\left\|{v_{h}}\right\|}_{Y(\mathcal{T}_{h})}}\geq\gamma_{I}{\left\|{u_{h}}\right\|}_{X(\mathcal{T}_{h})}\qquad\qquad\forall
u_{h}\in Y_{h}.$

Therefore, method (22) is well posed. An inspection of [13, Proposition 3.2]
reveals that the constant $\gamma_{I}$ only depends on the stability constants
$\alpha_{*}$ and $\alpha^{*}$ in (18).

##  3 Handling space–time meshes and variable degrees

We focus on the handling of the general prismatic mesh structure. More
precisely, we describe refinements procedures in Section 3.1; a time-slab
flagging strategy to split method (22) into smaller linear systems in Section
3.2; an element-topology flagging strategy in Section 3.3. Further, we discuss
space–time virtual elements with variable degrees of accuracy in Section 3.4.

###  3.1 Mesh refinements

We describe a procedure to refine general prismatic space–time elements. For
the sake of presentation, we consider the $(1+1)$ dimensional case only; the
extension to any spatial dimensions follows with a minor effort.

As discussed in Section 2.1, each mesh consists of rectangular elements
$\mathcal{K}$, with boundary given by the union of four (two space-like, two
time-like) straight segments. Each straight segment may be the union of
aligned edges of the element, which constitute the boundary of an associated
closed polytope $K$; see Section 2.1. Similar to [5, Figure 7], regardless of
the number of existing hanging nodes from previous refinements, a given $K$ to
be refined is split into four siblings by connecting the centroid of $K$ with
the midpoints of each straight segment of the boundary.

In Figure 2, we show an example of an element refinement. By this procedure,
at most new five nodes are generated, fewer in presence of previously
generated hanging nodes.

Figure 2: Element refinement strategy. We connect the centroid of the element
with the midpoints of each straight segment of the boundary, regardless of the
presence of previously generated hanging nodes. The red dots denote the nodes
of $K$; the blue dots denote the newly created nodes.

###  3.2 Time-slab flagging strategy

We present a flagging strategy that allows for the decomposition of the
algebraic linear system stemming from (22) into smaller linear systems. To
this aim, we assume that the first mesh of a sequence $\\{\mathcal{T}_{h}\\}$
is a “tensor-product-in-time” mesh, which can be arranged into time-slabs. On
this mesh, method (22) can be assembled and solved sequentially with respect
to the time-slabs, i.e., it can be interpreted as a time-stepping scheme.

Assume now that we are given a mesh at the refinement step $n-1$, which is
split into time-slabs. We explain here how to identify a time-slab partition
of the refined mesh at step $n$. The time interval $[0,T]$ can be partitioned
as $0=:t_{0}<\dots<t_{\ell}:=T$, so that there exists at least one element
$\mathcal{K}=\mathcal{K}_{\mathbf{x}}\times\mathcal{K}_{t}$ of the mesh at
step $n$ with $\mathcal{K}_{\mathbf{x}}$ contained in
$\Omega\times\\{t_{j}\\}$ for a $j=0,\dots,\ell-1$, and with
$\mathcal{K}_{t}=(t_{i},t_{j})$ for given $i,j=0,\dots,\ell$, $i<j$. The time-
slabs may be generated cycling on the $t_{j}$, $j=1,\dots,\ell-1$: if all the
elements in the mesh at step $n$ are either below or above $t_{j}$, then a new
time-slab is created by flagging all the elements that lie below the given
$t_{j}$ but have not been already allocated in a previously identified time-
slab.

Based on the new flagging, at the $n$-th refinement step, we can assemble and
solve method (22) sequentially on the newly created time-slabs. In Figure 3,
we illustrate this flagging procedure with an example.

xt$i$$i$$i$$i$$i$$ii$$ii$

xt$i$$i$$ii$$ii$$i$$i$$ii$$ii$$iii$$iii$

Figure 3:  Left panel: we start with a given space–time mesh with prescribed
time-slab structure and flagging. Right panel: time-slab flagging based on the
proposed strategy after one refinement step.

###  3.3 Element-topology flagging strategy

A common issue in the implementation of polytopic methods is the lack of a
reference element; this renders the computation of the local matrices more
expensive than for methods based on simplicial or Cartesian meshes.

Assume that the elements of a given mesh $\mathcal{T}_{h}$ can be grouped into
a uniformly finite number of equivalence classes (up to dilations and
translations), i.e., there exists a set $\mathcal{C}$ of “reference elements”,
with $\operatorname{card}(\mathcal{C})$ bounded uniformly for all meshes, such
that each element $K$ of $\mathcal{T}_{h}$ is equivalent to a “reference
element” in $\mathcal{C}$. In this case, local matrices need to be computed
only for the “reference elements”, thus hastening the assembling of the final
system. In adaptive mesh refinements, different “reference elements” may
appear. For this reason, we introduce a flag associated with the element
topology that identifies the corresponding “reference element”. Similar ideas
were used in [10] in a space-only VEM context.

For instance, in Figure 4 (left panel), we only have one “reference element”
(a space–time square); in Figure 4 (central panel), we have two types of
elements (a space–time square and a space–time pentagon with two time-like
facets on the left; in Figure 4 (right panel), we have two types of elements
(a space–time square and a space–time pentagon with two space-like facets at
the bottom. Inside each element, we denote the associated “topology-flag” with
a natural number.

xt$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$

xt$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$2$$2$

xt$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$3$$3$

Figure 4: Left panel: a mesh consisting of equivalent elements. Central panel:
a mesh with nonmatching time-like facets. Right panel: a mesh with nonmatching
space-like facets.

In Figure 5, we present three situations that are more elaborated than those
in Figure 4.

The above definition of equivalence classes refers to the geometry of the mesh
elements. This is enough to define a flagging strategy in case of a uniform
degree of accuracy. For the case of variable degrees of accuracy in Section
3.4 below, this definition can be extended taking into account the
distribution of degrees of accuracy. Namely, we flag two elements in the same
way if the following conditions are satisfied: they have the same geometric
flag; they have the same degree of accuracy; the degrees of accuracy assigned
to their time-like facets are the same. The existence of a maximum degree of
accuracy and a maximum number of time-like facets for all elements would allow
for a uniformly bounded number of reference elements in the extended sense.
Without a fixed maximum degree of accuracy (e.g., in the $p$-version of the
method, where $p\to\infty$), one cannot expect to have a finite number of
reference elements; this is also the case for the $p$-version of any standard
finite element discretization.

xt$1$$1$$1$$1$$1$$1$$1$$1$$2$$4$$1$$1$$5$$1$$1$$3$$3$$3$$3$

xt$1$$1$$1$$1$$1$$1$$1$$1$$2$$4$$1$$5$$5$$1$$3$$3$$3$$3$$3$$3$$3$$3$

xt$1$$1$$1$$1$$1$$2$$4$$2$$4$$1$$1$$5$$1$$1$$3$$3$$3$$3$$3$$3$$3$$3$

Figure 5: Three different meshes with element-topology flagging.

###  3.4 Space–time virtual element spaces with variable degrees of accuracy

Let $\mathcal{T}_{h}$ be a given space–time polytopic mesh consisting of
$N_{\mathcal{T}_{h}}$ elements and
$\mathbf{p}\in\mathbb{N}^{N_{\mathcal{T}_{h}}}$ be a given distribution of
degrees of accuracy. More precisely, we sort the elements of $\mathcal{T}_{h}$
as $\\{\mathcal{K}_{j}\\}_{j=1}^{N_{\mathcal{T}_{h}}}$, and denote the degree
of accuracy in each element $\mathcal{K}_{j}$ by $p_{j}$,
$j=1,\dots,N_{\mathcal{T}_{h}}$.

Given the vector $\mathbf{p}$, we fix the degrees of freedom associated with
each element according to the following maximum strategy:

  * •

in each element $\mathcal{K}_{j}$, $j=1,\dots,N_{\mathcal{T}_{h}}$, we take
bulk moments (8) up to degree $p_{j}-1$;

  * •

on each internal _time-like_ facet $F$ shared by two different elements
$\mathcal{K}_{j}$ and $\mathcal{K}_{\ell}$ for given
$j,\ell=1,\dots,N_{\mathcal{T}_{h}}$, we take _time-like_ moments (9) up to
degree $\max(p_{j},p_{\ell})$;

  * •

on each boundary _time-like_ facet $F$ on the boundary of the element
$\mathcal{K}_{j}$ for a given $j=1,\dots,N_{\mathcal{T}_{h}}$, we take _time-
like_ moments (9) up to degree $p_{j}$;

  * •

if
$\mathcal{K}_{j}=\mathcal{K}_{\mathbf{x}}{}_{,j}\times\mathcal{K}_{t}{}_{,j}$
for a given $j=1,\dots,N_{\mathcal{T}_{h}}$, we take the _space-like_ moments
(10) on $\mathcal{K}_{\mathbf{x}}{}_{,j}$ up to degree $p_{j}$.

We collect the _time-like_ polynomial degrees in the vector
$\mathbf{p}^{\text{time}}\in\mathbb{N}^{\operatorname{card}(\mathcal{F}_{h}^{\text{time}})}$
and the _space-like_ polynomial degrees in the vector
$\mathbf{p}^{\text{space}}\in\mathbb{N}^{N_{\mathcal{T}_{h}}}$. We order the
_time-like_ facets in $\mathcal{F}_{h}$ as
$\\{F_{j}\\}_{j=1}^{\operatorname{card}(\mathcal{F}_{h}^{\text{time}})}$.

Given $\widetilde{c}_{H}^{\mathcal{K}}$ and $\widetilde{\nu}^{\mathcal{K}}$ as
in (6), the corresponding local space on $K_{j}$ reads

$\begin{split}V_{h}(K_{j}):=\Big{\\{}v\in L^{2}(\mathcal{K}_{j})\mid\
&\widetilde{c}_{H}^{\mathcal{K}}\partial_{t}{v_{h}}-\widetilde{\nu}^{\mathcal{K}}\Delta_{\mathbf{x}}v_{h}\in\mathbb{P}_{p_{j}-1}\left(\mathcal{K}_{j}\right);\;v_{h}{}_{|\mathcal{K}_{\mathbf{x}}{}_{,j}}\in\mathbb{P}_{p_{j}^{\text{space}}}\left(\mathcal{K}_{\mathbf{x}}{}_{,j}\right);\\\\[1.99997pt]
&\mathbf{n}_{F_{k},\mathbf{x}}\cdot\nabla_{\mathbf{x}}v_{h}{}_{|F_{k}}\in\mathbb{P}_{p_{k}^{\text{time}}}\left(F_{k}\right)\,\forall
F_{k}\in\mathcal{F}_{K_{j}}^{\text{time}}\Big{\\}}.\end{split}$

The global space $Y_{h}$ is constructed by the nonconforming coupling of the
_time-like_ degrees of freedom (9). An immediate consequence of this maximum
strategy is that $\mathbb{P}_{p_{j}}(\mathcal{K}_{j})$ is contained in
$V_{h}(K_{j})$.

To illustrate the maximum strategy, we provide an example in Figure 6, where
we consider a uniform Cartesian mesh of $4$ elements with different degrees of
accuracy.

xt$1$$2$$3$$4$

xt$1$$2$$3$$4$$1$$2$$2$$3$$4$$4$

xt

Figure 6: Left panel: initial distribution of degrees of accuracy over the
elements. Central panel: “polynomial degrees” on _space-like_ (blue) and
_time-like_ (red) facets. Right panel: corresponding degrees of freedom; the
orange dots denote the bulk moments; the blue dots denote the _space-like_
moments; the red dots denote the _time-like_ moments.

##  4 Numerical investigation: convergence tests

We assess the convergence of the $h$\- and $hp$-versions of the virtual
element method (VEM) in (22). Since the virtual element solution $u_{h}$ is
not known in closed form and the error in the $X(\mathcal{T}_{h})$ norm is not
computable, we report the following error quantities: given $u_{h}$ the
solution to (22),

$\begin{split}&\mathcal{E}^{Y}:={\left\|{u-\Pi^{N}_{p}u_{h}}\right\|}_{Y(\mathcal{T}_{h})},\qquad\qquad\mathcal{E}^{N}:={\left\|{\Pi^{N}_{p}(\mathfrak{N}_{h}\Pi^{\star}_{p}(u-u_{h}))}\right\|}_{Y(\mathcal{T}_{h})},\\\
&(\mathcal{E}^{U})^{2}:=\frac{c_{H}}{2}\Bigg{(}{\left\|{\Pi^{\star}_{p}(u-u_{h})(\cdot,0)}\right\|}_{0,\Omega}^{2}+\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\scriptsize\begin{tabular}[]{c}$K_{\mathbf{x}}\in\mathcal{F}_{\mathcal{K}}^{\text{space}}$\\\
$K_{\mathbf{x}}\not\subset\Omega\times\\{0\\}$\end{tabular}}\\!\\!\\!\\!\\!\\!\\!\\!{\left\|{\mathcal{U}^{K_{\mathbf{x}}}(\Pi^{\star}_{p}(u-u_{h}))}\right\|}_{0,K_{\mathbf{x}}}^{2}+{\left\|{\Pi^{\star}_{p}(u-u_{h})(\cdot,T)}\right\|}_{0,\Omega}^{2}\Bigg{)},\\\
&\big{(}\mathcal{E}^{X}\big{)}^{2}:=\big{(}\mathcal{E}^{Y}\big{)}^{2}+\big{(}\mathcal{E}^{N}\big{)}^{2}+\big{(}\mathcal{E}^{U}\big{)}^{2}.\end{split}$
(23)

###  4.1 Test cases

We consider test cases with coefficients $\nu=1$ and $c_{H}=1$. The right-hand
side $f$, and the boundary and initial conditions are computed accordingly to
the exact solutions below.

#### Test case 1.

We define the analytic function

$\displaystyle u_{1}(x,t):=\exp(-t)\sin(\pi x)$ $\displaystyle\forall(x,t)\in
Q_{T}:=(0,1)\times(0,1).$ (24)

#### Test case 2.

For $\alpha>\frac{1}{2}$, we define the function

$\displaystyle u_{2}(x,t):=\sin(\pi x)\ t^{\alpha}$
$\displaystyle\forall(x,t)\in Q_{T}:=(0,1)\times(0,0.1),$ (25)

which belongs to
$H^{\alpha+1/2-\varepsilon}(0,1;\mathcal{C}^{\infty}(\Omega))$,
$\varepsilon>0$.

#### Test case 3.

We define the function

$\displaystyle
u_{3}(x,t):=\sum_{n=0}^{\infty}\frac{4}{(2n+1)\pi}\sin((2n+1)\pi\
x)\exp(-(2n+1)^{2}\pi^{2}t)$ $\displaystyle\forall(x,t)\in
Q_{T}:=(0,1)\times(0,1),$ (26)

which is the Fourier series of the solution to (3) with zero source term $f$,
initial condition $u_{0}=1$, and homogeneous Dirichlet boundary conditions
$g$. In the numerical experiments, the series in (25) is truncated at $n=250$.
The function $u_{3}$ belongs to $H^{s}(0,1;H^{1}_{0}(0,1))$ for any $s<1/4$
and to $H^{1}(0,1;H^{-1}(\Omega))\cap L^{2}(0,1;H_{0}^{1}(\Omega))\backslash
H^{1}(Q_{T})$; see [19]. In particular, $u_{3}$ is singular at the interface
of the (incompatible) initial and boundary conditions.

###  4.2 $h$\- and $hp$-versions for singular solutions

The performance of the $h$-version of the method on smooth solutions was
investigated in [13] and is therefore omitted here. The test case 1 is used in
Section 5 below.

We focus on the convergence of the $h$\- and $hp$-versions of the method for
singular solutions. To that aim, we consider the test cases 2 and 3. Notably,
we want to assess exponential convergence in terms of the cubic root of the
number of degrees of freedom for the $hp$-version on certain geometrically
refined space–time meshes.

First, we consider the test case 2. For the $h$-version of the method, we
consider uniform degree of accuracy $p=1$ and a sequence of uniform Cartesian
space–time meshes with $h_{t}=2h_{x}=0.2\times 2^{-i},\ i=1,\ldots,6$. For the
$hp$-version of the method, we proceed similarly as in [4, Example 2]: we fix
a partition of the spatial domain with $h_{x}=0.05$ and consider a sequence of
temporal meshes geometrically graded towards $t=0$ with grading factor
$\sigma_{t}=0.1$. In addition, the degree of accuracy $p$ is increased by $1$
from one time slab to the next one. In Figure 7, we depict the first three
meshes with varying degrees of accuracy.

Figure 7: First three meshes employed in the $hp$ refinements for the test
case 2 with exact solution $u_{2}$ in (25). The space–time domain is
$Q_{T}=(0,1)\times(0,0.1)$. For a better understanding of the figure, we scale
the $t$-coordinates by 10. In colours, we represent the local degrees of
accuracy: red: $p=1$; blue: $p=2$; green: $p=3$.

In Figures 8 and 9, we show the errors in (23) in _semilogy_ scale for
$\alpha=0.55$ and $\alpha=0.75$, respectively. Exponential convergence in
terms of the cubic root of the number of degrees of freedom is observed for
the $hp$-version of the method and both values of $\alpha$. In all cases, the
$h$-version is outperformed and displays only an algebraic decay of the error.

Figure 8: Convergence of the errors in (23) for the $h$\- and $hp$-versions of
the method. We consider the test case 2 with exact solution $u_{2}$ in (25),
$\alpha=0.55$.

Figure 9: Convergence of the errors in (23) for the $h$\- and $hp$-versions of
the method. We consider the test case 2 with exact solution $u_{2}$ in (25),
$\alpha=0.75$.

Next, we focus on the test case $3$. For the $h$-version of the method, we
consider uniform degree of accuracy $p=1$, and a sequence of uniform Cartesian
space–time meshes with $h_{x}=0.5h_{t}=2^{-i},\ i=1,\ldots,8$. For the
$hp$-version of the method, we proceed similarly as in [19, Section 7.4]: we
consider a sequence of space–time meshes geometrically graded towards $x=0,\
x=1$ and $t=0$ with grading factors $\sigma_{x}=\sigma_{t}=0.25$. In Figure
10, we depict the first three meshes with varying degrees of accuracy.

Figure 10: First three meshes employed in the $hp$ refinements for the test
case 2 with exact solution $u_{2}$ in (25). The space–time domain is
$Q_{T}=(0,1)\times(0,1)$. In colours, we represent the local degrees of
accuracy: red: $p=1$; blue: $p=2$; green: $p=3$.

Figure 11: Convergence of the errors in (23) for the $h$\- and $hp$-versions
of the method. We consider the test case 3 with exact solution $u_{3}$ in
(26).

In Figure 11, we show the errors in (23) in _semilogy_ scale. Exponential
convergence in terms of the cubic root of the number of degrees of freedom is
observed for the $hp$-version of the method; only algebraic convergence is
observed for the $h$-version.

##  5 Numerical investigation: an adaptive procedure

We consider a standard adaptive algorithm of the form

$\textbf{SOLVE}\quad\Longrightarrow\quad\textbf{ESTIMATE}\quad\Longrightarrow\quad\textbf{MARK}\quad\Longrightarrow\quad\textbf{REFINE}.$
(27)

We base the MARK step on a Dörfler marking strategy [7] with parameter
$\theta$ that we shall specify at each occurrence:

$\text{mark a set }\mathcal{A}\qquad\text{ such that
}\qquad\sum_{K\in\mathcal{A}}\eta_{K}^{2}\geq\theta\eta^{2}.$

The REFINE step involves only space–time mesh element refinements as described
in Section 3.1. For the ESTIMATE step, we use the following local, computable
residual-type error indicator: given an element
$\mathcal{K}\in\mathcal{T}_{h}$,

$\eta_{K}^{2}:=\sum_{i=1}^{5}\eta_{K,i}^{2},$

where

$\begin{split}\eta_{K,1}^{2}&:=\nu^{-1}\frac{h_{\mathcal{K}_{\mathbf{x}}}^{2}}{p^{2}}{\left\|{f+\nu\Delta_{\mathbf{x}}\Pi^{N}_{p}u_{h}-c_{H}\partial_{t}{\Pi^{\star}_{p}}u_{h}}\right\|}_{0,\mathcal{K}}^{2},\\\
\eta_{K,2}^{2}&:=\frac{1}{2}\nu^{-1}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\scriptsize\begin{tabular}[]{c}$F\in\mathcal{F}_{\mathcal{K}}^{\text{time}},$\\\
$F\not\subset\partial\Omega\times(0,T)$\end{tabular}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{h_{F_{\mathbf{x}}}}{p}{\left\|{\nu\left[\\!\left[\nabla_{\mathbf{x}}\Pi^{N}_{p}u_{h}\right]\\!\right]}\right\|}_{0,F}^{2},\\\
\eta_{K,3}^{2}&:=\nu\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\scriptsize\begin{tabular}[]{c}$F\in\mathcal{F}_{\mathcal{K}}^{\text{time}},$\\\
$F\subset\partial\Omega\times(0,T)$\end{tabular}}\\!\\!\\!\\!\\!\\!\\!\\!ph_{F_{\mathbf{x}}}^{-1}{\left\|{\Pi^{N}_{p}u_{h}}\right\|}_{0,F}^{2}+\frac{1}{2}\nu\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\scriptsize\begin{tabular}[]{c}$F\in\mathcal{F}_{\mathcal{K}}^{\text{time}},$\\\
$F\not\subset\partial\Omega\times(0,T)$\end{tabular}}\\!\\!\\!\\!\\!\\!\\!\\!ph_{F_{\mathbf{x}}}^{-1}{\left\|{\left[\\!\left[\Pi^{N}_{p}u_{h}\right]\\!\right]}\right\|}_{0,F}^{2},\\\
\eta_{K,4}^{2}&:=\sum_{K_{\mathbf{x}}\in\mathcal{F}_{\mathcal{K}}^{\text{space}}}c_{H}^{-1}{\left\|{\mathcal{U}^{K_{\mathbf{x}}}(u_{h})}\right\|}^{2}_{0,K_{\mathbf{x}}},\\\
\eta_{K,5}^{2}&:=\nu
S^{K}((I-\Pi^{N}_{p})u_{h},(I-\Pi^{N}_{p})u_{h}).\end{split}$ (28)

The above local residual-type error indicator consists of five terms:
$\eta_{K,1}$ is the internal residual of the projected discrete solution;
$\eta_{K,2}$ is the boundary residual involving the normal trace of the
gradient of the projected discrete solution; $\eta_{K,3}$ is a term due to the
nonconformity in space involving the jump of traces on time-like facets;
$\eta_{K,4}$ is related to the upwind terms in the scheme; $\eta_{K,5}$ is a
correction term due to the virtual element stabilization of the method.

The global error indicator and its parts read

$\eta^{2}:=\sum_{i=1}^{5}\eta_{i}^{2},\qquad\qquad\eta_{i}^{2}:=\sum_{\mathcal{K}\in\mathcal{T}_{h}}\eta_{K,i}^{2}.$
(29)

###  5.1 Assessment of the reliability and efficiency of the error indicator

We test the performance of the proposed error indicator $\eta$ in (29). To
this aim, we introduce the effectivity index related to the computable error
$\mathcal{E}^{Y}$:

$\text{effectivity index}:=\frac{\eta}{\mathcal{E}^{Y}}.$ (30)

In Figure 12 (left panel), we assess numerically the reliability and
efficiency of $\eta$ for the test case $1$ with exact solution $u_{1}$ in (24)
under uniform mesh refinements, starting from a uniform space–time Cartesian
mesh with $h_{x}=h_{t}=0.1$. The effectivity indices for $p=1,2,3$, tend to
constant values, whence the error indicator $\eta$ appears to be efficient and
reliable with respect to the error $\mathcal{E}^{Y}$.

In Figure 12 (right panel), we show all the errors in (23) and the five terms
appearing in the error indicator (29) for degree of accuracy $p=1$. The error
indicator $\eta$ decays with order $\mathcal{O}(N_{DoFs}^{-\frac{1}{2}})$),
i.e., slower than that for the error $\mathcal{E}^{N}$ in (23), which decays
with order $\mathcal{O}(N_{DoFs}^{-1})$. This suggests that the error
indicator $\eta$ is not efficient with respect to the error $\mathcal{E}^{N}$.

Figure 12: The test case 1 with exact solution $u_{1}$ in (24). Left panel:
Effectivity index. Right panel: Comparison of the errors in (23) with the
terms appearing in the error indicator (29) for $p=1$.

Next, we focus on the test cases 2 and 3 with singular solutions $u_{2}$
($\alpha=0.55$ and $0.75$) in (25) and $u_{3}$ in (26) on the space–time
domains $Q_{T}=(0,1)\times(0,1)$ and $Q_{T}=(0,1)\times(0,0.1)$, respectively.
In Figures 13–15, we show the decay of the errors in (23) and of the five
terms appearing in the error indicator under uniform mesh refinements,
starting with Cartesian meshes with $h_{x}=10h_{t}=0.1$ (for the test case 2)
and $h_{x}=h_{t}=0.1$ (for the test case 3). In all cases, the effectivity
indices shown on the left panels tend to constant values, which suggests that
the error indicator (29) is efficient and reliable with respect to the error
$\mathcal{E}^{Y}$ also for singular solutions.

For the test case 2, the errors $\mathcal{E}^{Y}$ and $\mathcal{E}^{N}$ decay
with orders $\mathcal{O}(N_{DoFs}^{-\frac{1}{2}(\alpha+\frac{1}{2})})$ and
$\mathcal{O}(N_{DoFs}^{-\frac{1}{2}(\alpha-\frac{1}{2})})$, respectively; see
Figures 13 and 14. The error indicator $\eta$ decays with the same order as
that of the error $\mathcal{E}^{Y}$, while it it not reliable with respect to
the error $\mathcal{E}^{N}$.

For the test case 3, the error indicator $\eta$ and all the errors in (23)
decay with the same order, namely $\mathcal{O}(N_{DoFs}^{-\frac{s}{2}})$ with
$s=\frac{1}{4}$; see Figure 15.

In summary, the above experiments seem to indicate that the error indicator
$\eta$ is reliable and efficient for the error $\mathcal{E}^{Y}$ but not for
the error $\mathcal{E}^{N}$.

Figure 13:  The test case 2 with exact solution $u_{2}$ in (25) and
$\alpha=0.55$. Left panel: Effectivity index. Right panel: Comparison of the
errors in (23) with the terms appearing in the error indicator (29) for $p=2$

.

Figure 14: The test case 2 with exact solution $u_{2}$ in (25) and
$\alpha=0.75$. Left panel: Effectivity index. Right panel: Comparison of the
errors in (23) with the terms appearing in the error indicator (29) for $p=2$

.

Figure 15: The test case 3 with exact solution $u_{3}$ in (26). Left panel:
Effectivity index. Right panel: Comparison of the errors in (23) with the
terms appearing in the error indicator (29) for $p=2$.

###  5.2 Adaptive mesh refinements

We test the performance of method (22) under adaptive mesh refinements as
described in (27). We consider the test cases 2 and 3 with exact solutions
$u_{2}$ ($\alpha=0.55$) in (25) and $u_{3}$ in (26), respectively. The marking
step is dictated by the error indicator in (29).

We are also interested in comparing the results with those obtained with an
adaptive procedure for the continuous finite element method (FEM) of (5) in
[21]:

$\begin{cases}\text{find }\widetilde{u}_{h}\in\widetilde{X}_{h}\text{ such
that}\\\
b(\widetilde{u}_{h},\widetilde{v}_{h})=(f,\widetilde{v}_{h})_{0,Q_{T}}\quad\quad\forall\widetilde{v}_{h}\in\widetilde{Y}_{h}.\end{cases}$
(31)

Above, $\widetilde{Y}_{h}$ is a space of continuous piecewise polynomials over
a space–time simplicial tessellation of $Q_{T}$ and $\widetilde{X}_{h}$ is the
subspace of $\widetilde{Y}_{h}$ of functions with zero initial condition. We
recall the residual-type error indicator introduced in [22, 23]:

$\begin{split}\widetilde{\eta}^{2}:=\sum_{i=1}^{2}\widetilde{\eta}_{i}^{2},&\qquad\widetilde{\eta}_{i}^{2}=\sum_{K\in\mathcal{T}_{h}}\widetilde{\eta}_{K,i}^{2},\\\
\widetilde{\eta}_{K,1}^{2}:=\frac{h_{K}^{2}}{p^{2}}{\left\|{f+\nu\Delta_{\mathbf{x}}\widetilde{u}_{h}-c_{H}\partial_{t}{\widetilde{u}_{h}}}\right\|}_{0,K}^{2},&\qquad\widetilde{\eta}_{K,2}^{2}:=\frac{1}{2}\sum_{F\in\mathcal{F}_{\mathcal{K}}^{\text{time}}}\frac{h_{F}}{p}{\left\|{\nu\left[\\!\left[\nabla_{\mathbf{x}}\widetilde{u}_{h}\right]\\!\right]}\right\|}_{0,F}^{2},\end{split}$
(32)

and the error quantity

$E^{Y}={\left\|{u-\widetilde{u}_{h}}\right\|}_{Y}.$ (33)

For the VEM, we start with a mesh with $1$ element; for the continuous FEM, we
start with a structured simplicial mesh with $2$ elements.

In Figure 16, we show the errors of both methods under uniform and adaptive
mesh refinements for the test case $2$ with $\alpha=0.55$. Uniform and
adaptive mesh refinements for the VEM in (22) lead to higher convergence rates
than those for the continuous FEM in (31). For the Dörfler marking strategy,
we set $\theta=0.99$ for the VEM and $\theta=0.9$ for the FEM.

Figure 16: The test case 2 with exact solution $u_{2}$ in (25) and
$\alpha=0.55$. Left panel: Error $E^{Y}$ for the continuous FEM with $p=2$
under uniform and adaptive mesh refinements, with convergence rates of
approximately $\mathcal{O}(N_{DoFs}^{-0.28})$ and
$\mathcal{O}(N_{DoFs}^{-0.40})$, respectively. Central panel: Error
$\mathcal{E}^{Y}$ for the VEM with $p=2$ under uniform and adaptive mesh
refinements, with convergence rates of approximately
$\mathcal{O}(N_{DoFs}^{-0.52})$ and $\mathcal{O}(N_{DoFs}^{-1})$,
respectively. Right panel: The effectivity index for the adaptive VEM.

In Figure 17, we show the $Y$-type errors for both methods under uniform and
adaptive mesh refinements with Dörfler marking parameter $\theta=0.9$ for the
test case $3$ and $p=1$. The adaptive procedure for the VEM (22) leads to
higher convergence rates compared to those obtained for uniform refinements.
For the continuous FEM (31), although the adaptive procedure produces meshes
that are refined towards the bottom corners as in Figure 19 (right panels),
the error $E^{Y}$ does not converge to zero. This experiment seems to suggest
that the presented nonconforming approach is able to capture correctly
singularities arising from incompatible initial and boundary data.

Figure 17: The test case 3 with exact solution $u_{3}$ in (26). Left panel:
Error $E^{Y}$ for the continuous finite element with $p=1$ under uniform and
adaptive refinements. Central panel: Error $\mathcal{E}^{Y}$ for the VEM with
$p=1$ under uniform and adaptive mesh refinements, with convergence rates of
approximately $\mathcal{O}(N_{DoFs}^{-0.13})$ and
$\mathcal{O}(N_{DoFs}^{-33})$, respectively. Right panel: The effectivity
index for the adaptive VEM.

In Figures 18 and 19, we plot some meshes for the test cases with exact
solutions $u_{2}$ ($\alpha=0.55$) and $u_{3}$, respectively, produced by the
adaptive procedure driven by the VEM error indicator in (29) and the
continuous FEM error indicator in (32).

(a)

(b)

(c)

(d)

Figure 18: Meshes generated by the adaptive schemes driven by the VEM error
indicator $\eta$ in (29) (left panel) and the continuous FEM error indicator
$\widetilde{\eta}$ (right panel) for the test case with exact solutions
$u_{2}$ ($\alpha=0.55$) in (25).

(a)

(b)

(c)

(d)

Figure 19: Meshes generated by the adaptive schemes driven by the VEM error
indicator $\eta$ in (29) (left panel) and the continuous FEM error indicator
$\widetilde{\eta}$ in (32) (right panel) for the test case with exact
solutions $u_{3}$ in (26).

Next, in Tables 1 and 2, we focus on the adaptive mesh refinements driven by
the VEM error indicator in (29) and report the number of time-slabs and
“reference elements” as in Sections 3.2 and 3.3 for some adaptively generated
meshes;. The time-slab and the element-topology strategies allow us to hasten
considerably the assembling and the solving time of method (22).

Mesh | Number of reference elements | Total number of elements | Number of time-slabs
---|---|---|---
$m_{2}$ | 1 | 4 | 2
$m_{4}$ | 1 | 52 | 5
$m_{6}$ | 16 | 265 | 10
$m_{8}$ | 19 | 1189 | 21
$m_{10}$ | 24 | 5110 | 45
$m_{12}$ | 24 | 21883 | 103
Table 1: Performance of flagging strategies for the test case $2$ with $\alpha=0.55$, $p=2$ and adaptive refinements. Mesh | Number of reference elements | Total number of elements | Number of time-slabs
---|---|---|---
$m_{2}$ | 1 | 4 | 2
$m_{4}$ | 3 | 22 | 4
$m_{6}$ | 9 | 91 | 5
$m_{8}$ | 21 | 361 | 9
$m_{10}$ | 33 | 1204 | 13
$m_{12}$ | 35 | 3493 | 21
Table 2: Performance of flagging strategies for the test case $3$ with $p=1$
and adaptive refinements.

##  6 Conclusions

We extended the virtual element framework of [13] to the case of general
prismatic space–time meshes with possible hanging nodes and nonuniform degrees
of accuracy, useful in $hp$-adaptive procedures. We discussed flagging
strategies to handle the space–time meshes data structure; this improves the
performance of the method. We investigated numerically the $hp$-version of the
method and demonstrated the expected exponential convergence in terms of
suitable roots of the number of degrees of freedom for singular solutions. A
residual-type error indicator was introduced and an $h$-adaptive refinement
procedure was tested based on it. This error indicator appears to be reliable
and efficient for the error $\mathcal{E}^{Y}$ but not for the error
$\mathcal{E}^{N}$ in (23). We investigated numerically the overall virtual
element adaptive procedure. The results obtained for certain singular
solutions show that it outperforms a corresponding one for a continuous finite
element method.

## References

  * [1] A. Aziz and P. Monk.  Continuous finite elements in space and time for the heat equation.  Math. Comp., 52(186):255–274, 1989.
  * [2] I. Babuška and T. Janik.  The $h$-$p$ version of the finite element method for parabolic equations. Part I. The $p$-version in time.  Numer. Methods Partial Differential Equations, 5(4):363–399, 1989\.
  * [3] I. Babuška and T. Janik.  The $h$-$p$ version of the finite element method for parabolic equations. Part II. The $h$-$p$ version in time.  Numer. Methods Partial Differential Equations, 6(4):343–369, 1990\.
  * [4] A. Cangiani, Z. Dong, and E.H. Georgoulis.  $hp$-version space-time discontinuous Galerkin methods for parabolic problems on prismatic meshes.  SIAM J. Sci. Comput., 39(4):A1251–A1279, 2017.
  * [5] A. Cangiani, E. H. Georgoulis, T. Pryer, and O. J. Sutton.  A posteriori error estimates for the virtual element method.  Numer. Math., 137(4):857–893, 2017.
  * [6] R. Dautray and J.-L. Lions.  Mathematical Analysis and Numerical Methods for Science and Technology, volume 5, Evolution Problems I. Springer-Verlag, 1992.
  * [7] W. Dörfler.  A convergent adaptive algorithm for Poisson’s equation.  SIAM J. Numer. Anal., 33(3):1106–1124, 1996.
  * [8] K. Eriksson, C. Johnson, and V. Thomée.  Time discretization of parabolic problems by the discontinuous Galerkin method.  ESAIM Math. Model. Numer. Anal., 19(4):611–643, 1985.
  * [9] D. French and T. Peterson.  A continuous space-time finite element method for the wave equation.  Math. Comp., 65(214):491–506, 1996.
  * [10] M. Frittelli, A. Madzvamuse, and I. Sgura.  Virtual element method for elliptic bulk-surface PDEs in three space dimensions.  Numer. Methods Partial Differential Equations, 2023.
  * [11] T. Führer and M. Karkulik.  Space–time least-squares finite elements for parabolic equations.  Comput. Math. Appl., 92:27–36, 2021.
  * [12] G. Gantner and R. Stevenson.  Further results on a space-time FOSLS formulation of parabolic PDEs.  ESAIM Math. Model. Numer. Anal., 55(1):283–299, 2021.
  * [13] S. Gómez, L. Mascotto, A. Moiola, and I. Perugia.  Space-time virtual elements for the heat equation.  https://arxiv.org/abs/2212.05343, 2022.
  * [14] T. J. R. Hughes and J. R. Stewart.  A space-time formulation for multiscale phenomena.  J. Comput. Appl. Math., 74(1-2):217–229, 1996.
  * [15] P. Jamet.  Galerkin-type approximations which are discontinuous in time for parabolic equations in a variable domain.  SIAM J. Numer. Anal., 15(5):912–928, 1978.
  * [16] U. Langer, S. E. Moore, and M. Neumüller.  Space–time isogeometric analysis of parabolic evolution problems.  Comput. Methods Appl. Mech. Engrg., 306:342–363, 2016.
  * [17] M. Montardini, M. Negri, G. Sangalli, and M. Tani.  Space–time least–squares isogeometric method and efficient solver for parabolic problems.  Math. Comp., 89(323):1193–1227, 2020.
  * [18] A. Schafelner and P. S. Vassilevski.  Numerical results for adaptive (negative norm) constrained first order system least squares formulations.  Comput. Math. Appl., 95:256–270, 2021.
  * [19] D. Schötzau and C. Schwab.  Time discretization of parabolic problems by the $hp$-version of the discontinuous Galerkin finite element method.  SIAM J. Numer. Anal., 38(3):837–875, 2000.
  * [20] Ch. Schwab and R. Stevenson.  Space-time adaptive wavelet methods for parabolic evolution problems.  Math. Comp., 78(267):1293–1318, 2009.
  * [21] O. Steinbach.  Space-time finite element methods for parabolic problems.  Comput. Methods Appl. Math., 15(4):551–566, 2015.
  * [22] O. Steinbach and H. Yang.  Comparison of algebraic multigrid methods for an adaptive space–time finite-element discretization of the heat equation in 3D and 4D.  Numer. Linear Algebra Appl., 25(3):e2143, 2018.
  * [23] O. Steinbach and H. Yang.  Space-time finite element methods for parabolic evolution equations: discretization, a posteriori error estimation, adaptivity and solution.  Space-Time Methods: Application to Partial Differential Equations, 25:207–248, 2019.
  * [24] O. Steinbach and M. Zank.  Coercive space-time finite element methods for initial boundary value problems.  Electron. Trans. Numer. Anal., 52:154–194, 2020.
  * [25] R. Stevenson and J. Westerdiep.  Stability of Galerkin discretizations of a mixed space–time variational formulation of parabolic evolution equations.  IMA J. Numer. Anal., 41(1):28–47, 2021.
  * [26] J.J. Sudirham, J.J.W. Van der Vegt, and R.M.J. Van Damme.  Space-time discontinuous Galerkin method for advection-diffusion problems on time-dependent domains.  Appl. Numer. Math., 56(12):1491–1518, 2006.
  * [27] K. Voronin, C. S. Lee, M. Neumüller, P. Sepulveda, and P. S. Vassilevski.  Space-time discretizations using constrained first-order system least squares (CFOSLS).  J. Comput. Phys., 373:863–876, 2018.