# Complexity of direct and iterative solvers on space-time formulations versus
time–marching schemes for h-refined grids towards singularities

Marcin Skotniczny(1), Anna Paszyńska(2), Sergio Rojas(3), Maciej Paszyński(1)
(1)AGH University of Science and Technology, Kraków, Poland
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(2) Jagiellonian University, Kraków, Poland
e-mail<EMAIL_ADDRESS>
(3) Instituto de Matemáticas. Pontificia Universidad Católica de Valparaíso,
Chile
e-mail<EMAIL_ADDRESS>
###### Abstract

We study computational complexity aspects for Finite Element formulations
considering hypercubic space–time full and time–marching discretization
schemes for $h$–refined grids towards singularities. We perform a relatively
comprehensive study of comparing the computational time via time complexities
of direct and iterative solvers. We focus on the space-time formulation with
refined computational grids and on the corresponding time slabs, namely,
computational grids obtained by taking the “cross-sections” of the refined
space-time mesh. We compare the computational complexity of the space-time
formulation and the corresponding time–marching scheme. Our consideration
concerns the computational complexity of the multi-frontal solvers, the
iterative solvers, as well as the static condensation. Numerical experiments
with Octave confirm our theoretical findings.

###### keywords:

space-time formulation , time marching scheme , $h$ refinements ,
computational complexity , multi-frontal solver , iterative solver

##  1 Introduction

Several approaches exist in the Finite Element (FE) community to deal with
time-dependent Partial Differential Equations (PDEs). Nowadays, two common
approaches are classical time–marching FE schemes (see e.g. [2]) and space-
time formulations, where time is treated as another spatial variable (see e.g.
[1]). Both methodologies enjoy several advantages, implying in general that
one can run with some benefit concerning the other depending on the result of
interest. For instance, if only the final state is of interest, time–marching
schemes are typically preferred due to their flexibility in the election for
schemes to be considered in the time discretization procedure and their
optimal storage requirements as a consequence of their sequential nature.
However, when considering problems where the continuous history of the
evolution problem is of interest, typically, a space-time formulation will be
preferred, despite the natural increment of the complexity associated with its
resolution due to the dimensional increment in the discrete formulation.
Moreover, space-time formulations run with an advantage in the presence of
space–time singularities since they allow for the construction of robust a
posteriori error estimates for local space-time mesh refinements.
Time–marching schemes must remain with the sequential nature of the
refinement, implying a possible increment in the storage requirements
depending on the nature of the singularity. Finally, space-time formulations
allow for parallel static condensation, thus possibly being faster than
time–marching schemes when a sufficient amount of computational power is
available. Even though there is still extensive research on the space-time
formulations and the time–marching schemes, to the authors’ knowledge, there
are no works related to theoretical aspects of the computational complexity of
both methodologies available in the literature.
In this work, we study complexity aspects of time–marching and space-time
formulations, being the first step to further work in complexity aspects. We
explore the complexity of iterative and direct solvers when considering
hypercubic $h$–refined grids towards singularities. For this purpose, we
propose a general method-independent strategy simulating the best-scenario
possibility for each case. We proceed as follows: To simulate an $h$–adaptive
space–time procedure, we start from a uniform $d$–dimensional mesh, where
$d-1$ dimensions correspond to the spatial discretization and the last
dimension to the time discretization, and we consider regular hypercubic
geometrical refinements towards singularities. To simulate a time–marching
scheme, we consider a sequence of $d-1$ dimensional meshes obtained from
considering time cuts of the space–time grids, with a fixed time step of
length equal to the smallest element dimension in the time axis direction. To
motivate these assumptions, the reader can consider a parabolic problem, for
instance, the heat equation, as the PDE of interest. Therefore, the assumption
for the space–time formulation is natural, assuming that information of the
singularity is available, for instance, an a posteriori error estimator. In
contrast, the assumption for the time–marching scheme will be the ideal
scenario for the explicit-in-time Euler scheme ensuring that the
Courant–Friedrichs–Lewy condition (CFL condition) [3] is satisfied, and also
information on the singularity is provided.
To derive the estimations, we start by considering refinements toward a space-
time “edge”, resulting from a point traveling through space and time. This
space-time refinement pattern corresponds to a sequence of spatial meshes
refined towards “points” located on the space-time edge at particular time
moments. Next, we consider the refinement towards a space-time “face”,
resulting from an edge traveling through space and time. Again, this space-
time refinement pattern corresponds to a sequence of spatial meshes refined
towards “edges” located on the space-time face at a particular moment.
Finally, we consider the refinement towards a space-time “hyperface” resulting
from a face traveling through space and time. This space-time refinement
pattern corresponds to a sequence of spatial meshes refined towards “faces”
located on the space-time hyperface at a particular time instant. We compare
the computational complexity of the iterative solver executed over the $d$
dimensional space-time domain to the computational complexity of the iterative
solver performed multiple times over the $d$-1 dimensional meshes during the
time–marching scheme. We compare the computational complexity of the direct
solver run over the $d$ dimensional space-time domain to the computational
complexity of the direct solver executed multiple times over $d$-1 dimensional
meshes within the time–marching scheme.
We summarize the findings of this work in the following two points:

  * 1.

This paper derives a general formula for the number of unknowns for the mesh
of arbitrary dimension $d$, refined towards a singularity of dimension $q$.
This paper generalizes the results discussed in [22], for the case of point
singularity ($q=0$), edge singularity ($q=1$), face singularity ($q=2$),
hyperface singularity ($q=3$), focusing on the computational complexities of
space-time formulations and time–marching schemes. In our general estimates,
we do not consider the polynomial order of approximation $p$, and we assume
that this is a constant in our formulas. Additionally, the singularities in
our case may have an arbitrary shape, as is presented in Figure 13. Therefore,
we only fix the dimension of the space $d$ and the dimension of the
singularity $q$.

  * 2.

This paper estimates the time complexity of finite element method simulations
performed on adaptive $d$–dimensional space–time meshes. We also assess the
time complexity of finite element method simulations performed on a sequence
of adaptive $d-1$–dimensional meshes resulting from the corresponding
time–marching scheme. In particular, we estimate $N$ = the number of unknowns
for the whole adaptive space–time formulation, $n$= the number of unknowns
from a mesh from a sequence of adaptive meshes resulting from the
time–marching scheme. We estimate the time complexity of the iterative and
direct solvers. Thus, we have the lower bounds (expressed by the number of
unknowns) and upper bounds (described by the time complexity of the sparse
direct solver).

The structure of the paper is the following. We start Section 2 with some
preliminaries, including a general algorithm refining the mesh towards a given
singularity, followed by an overview of the sparse Gaussian elimination,
usually implemented using a multi-frontal solver approach, and our notion of
element partition tree that allows constructing an ordering for sparse matrix
permutation in the space–time formulations. We also discuss the recursive
formula for estimating the time complexity of our sparse direct solvers based
on arrangements constructed from the element partition trees. Section 3
summarizes our findings in the context of space–time formulations and
time–marching schemes. We also present numerical verification using Octave
codes in Section 4. The paper is concluded in the last section. In Appendix A,
we derive time complexity for a multi-dimensional grid refined toward point
singularity. In Appendix B, we derive time complexity for multi-dimensional
grids refined toward arbitrary singularities. In this case, we also include
the dependency in the polynomial approximation orders.

##  2 Preliminaries

###  2.1 Complexity literature review and assumptions

Computational complexity, especially time complexity, is one of the most
fundamental concepts of theoretical computer science. It was first defined in
1965 by Hartmanis and Stearns [4]. In this paper we estimate the computational
complexity of solving systems of linear equations for multidimensional meshes
with arbitrary dimensional similarities encountered in space –time formulation
as well as time–marching schemes. The systems of equations generated by non-
stationary problems can be solved by either direct [5, 6] or iterative solvers
[21]. The time complexity of iterative solvers, in general, can be estimated
as $n_{iter}\times N$, where $n_{iter}$ is the number of iterations, and in
general, it depends on the spectral properties of the matrix, and it grows
with the problem size $N$. The time complexity of direct solvers [5, 6] for
certain classes of meshes, especially regular meshes, is well known. In
particular, for three-dimensional uniform grids, the computational complexity
is of the order of ${\cal O}(N^{2})$. For the two-dimensional grids, the
complexity is of the order of ${\cal O}(N^{1.5})$ [7, 8]. The sparse direct
solvers rely on the smart permutation of the matrix, resulting in its banded
structure and efficient sparse factorization, avoiding zero entries. The
problem of finding an optimal order of elimination of unknowns for the direct
solver, in general, is indeed NP-complete [10]. There are several heuristic
algorithms analyzing the sparsity pattern of the resulting matrix [12, 13, 14,
15, 16]. We focus on adaptive grids, refined in space–time and space domains.
For three-dimensional grids adapted towards the point, edge, and face, the
computational complexities are ${\cal O}(N),{\cal O}(N)$, and ${\cal
O}(N^{1.5})$, respectively [9]. These estimates assume a prescribed order of
eliminating variables [17]. Similarly, for two-dimensional grids refined
towards a point or edge, it is ${\cal O}(N)$ [18]. Again, these estimates
assume a prescribed order of elimination of variables [19]. We generalized the
results of [9] into multidimensional adaptive grids with arbitrary dimensional
singularities. There are also solvers based on hierarchical matrices
decomposition [11]. Their time complexity generally depends on the number of
unknowns in the mesh. However, the constant in front of the complexity also
grows with the mesh dimension, sparsity of the matrix and rank of sub-blocks.
We leave their analysis in the space–time set up for our future work.

###  2.2 Hierarchical basis functions and construction of an $h$-adaptive mesh

For this work, we employ $d$-dimensional hypercube elements (rectangles in 2D,
hexahedrons in 3D, and octachorons in 4D), and we consider hierarchical
polynomial basis functions of order $p$. Hierarchical basis functions are
constructed by the tensor product of one-dimensional hierarchical shape
functions glued together to obtain globally continuous piecewise polynomial
functions with compact support. We refer to [20] for a more general
description.
We consider shape functions over vertices, edges, faces, and interiors in
three dimensions, while in four dimensions over vertices, edges, faces,
hyperfaces, and interiors.
We identify nodes of the mesh with basis functions. We also consider nodes’
support, defined as support of basis functions associated with the node. We
have basis functions assigned to vertices, edges, faces, hyperfaces (in higher
dimensions), and interiors. In general, the support of the nodal function
spans over all the adjacent elements sharing the node. For example, in the
case of two-dimensional regular mesh, the support of the vertex node spans
into four elements sharing the node, the support of the edge node spans into
two elements sharing the node, and the support of the interior node is equal
to the single element.
We focus on $h$-adaptive meshes, where we employ the 1-irregularity rule,
stating that an element can be broken only once without breaking its
neighbors. To construct an $h$-adaptive mesh around a singularity, we start
with the one initial element and iteratively refine all elements that overlap
with the singularity, ensuring that the $1$-irregularity rule is being
followed. For example, Algorithm 1 can be used to build $h$-adaptive mesh
around a singularity with shape $S$ with refinement level $R$, where the
symbol “ $[]$ ” stands for initiation as an empty array. We assume that
information on the singularity location is available. An exemplary singularity
mesh construction is shown in Figure 1.

Figure 1: Mesh refined over a point singularity resulting from a gradient of a
two-dimensional function. The shape of the singularity $S$ is a point in the
corner.

procedure Constructmesh($RootElement,R,S$)

$root\leftarrow RootElement$

$G\leftarrow[]$ $\triangleright$ $G[r]$ contains all elements of refinement
level $r$.

$G[0]\leftarrow[root]$ $\triangleright$ Initialize the mesh with single
element.

for $r\leftarrow 1$ to $R$ do $\triangleright$ Step 1: Refine all elements as
necessary.

$H\leftarrow G[r-1]$ $\triangleright$ Only elements of refinement $r-1$ can be
refined further.

$G[r-1]\leftarrow[]$

$G[r]\leftarrow[]$

for all $e\in H$ do

if $e\text{ overlaps }S$ then

if $i>1$ then

$K\leftarrow G[r-2]$

$G[r-2]\leftarrow[]$

for all $f\in K$ do

if $f\text{ shares a vector with }e$ then

$G[r-1]\leftarrow\textsc{concat}(G[r-1],\textsc{refine}(f))$

else

$G[r-2]\leftarrow\textsc{append}(G[r-2],f)$

end if

end for

end if

$G[r]\leftarrow\textsc{concat}(G[R],\textsc{refine}(e))$ $\triangleright$ If
analyzed element overlaps with the singularity, refine it into smaller
elements

else

$G[r-1]\leftarrow\textsc{append}(G[r-1],e)$ $\triangleright$ …otherwise, leave
it unrefined.

end if

end for

end for

$M\leftarrow[]$

for $r\leftarrow 0$ to  $R$ do

$M\leftarrow\textsc{concat}(M,G[r])$

end for

return $M$

end procedure

Algorithm 1 mesh construction over a singularity

###  2.3 Element partition tree

The element partition tree is used for the construction of ordering for sparse
matrix permutation in order to speed up the multi-frontal solver. Using some
partitioning strategy, an element partition tree is created by recursive
partitioning the mesh elements into two parts [17], [18],[19]. An example of
an element partition is presented in Figure 2.
An element partition tree for a mesh consisting of a set of elements $E$ is a
binary tree defined as $T=(E,V,c_{1},c_{2},e)$ with the following properties:

  1. 1.

$V$ is a set of tree nodes,

  2. 2.

$c_{1}:V\rightarrow V$ and $c_{2}:V\rightarrow V$ are the functions assigning
left and right child to a node, respectively,

  3. 3.

$e:V\rightarrow P(A)$ is a function assigning subsets of a set of all elements
of the mesh $E$ to a node,

  4. 4.

the root called node $ROOT$ contains all elements of the mesh, in other words,
$e(ROOT)=E$,

  5. 5.

each node $node$ for which $|e(node)|=1$ is a leaf in the tree,

  6. 6.

each node $node$ for which $|e(node)|>1$ has exactly two children; $c_{1}(n)$
and $c_{2}(n)$,

  7. 7.

for each node $node$ $e(node)=e(c_{1}(node))\cup e(c_{2}(node))$ .

Figure 2: An example for an element partition tree for a small adaptive mesh.
The post-order traversal of the tree results in the ordering of elements
$(10,9,7,8,6,4,5,3,1,2)$.

###  2.4 Ordering generation and solving using generated ordering

Given an element partition tree, we can generate a row elimination order for
the matrix using a post-order traversal of the element partition tree. At each
traversed partition tree node, we list all basis functions with support
entirely contained by the tree node elements which have not been listed
already. This produces a permutation (or ordering) of all nodes. Figure 3
presents an exemplary element partition tree with denoted basis functions for
each node, which generates the ordering
$(1,4,5,7,8,9,10,12,15,11,16,17,13,14,2,3,6)$.

Figure 3: An example of element partition tree ordering generation. Nodes
marked in green are eliminated at the given tree node. This example results in
the permutation of $(1,4,5,7,8,9,10,12,15,11,16,17,13,14,2,3,6)$.

s

###  2.5 Sparse Gaussian elimination

We focus on the Gaussian elimination adapted to work over sparse matrices. The
algorithm solves a linear equation system represented by $Mx=A$ ($M$ is a
sparse matrix and $A$ is a vector) by multiplying both sides by some ordering
$P$ and solving $PMx=PA$ ($P$ being a permutation matrix representing some
ordering). In our case, the permutation of the matrix is based on the post-
order traversal of the element partition tree. The practical implementation of
the algorithm is the multi-frontal solver [5, 6], avoiding zeros in matrices
by constructing the elimination tree internally based on the proposed ordering
and the sparsity pattern of the matrix. In our numerical experiments, we
employ the multi-frontal solver from Octave.

###  2.6 Time complexity of element partition tree based solvers

This section shows that an element partition tree-based ordering gives a
recursive formula for computation complexity that is easy to calculate. First,
let us analyze the Gaussian elimination algorithm. It comprises two main
steps–a first step corresponding to a forward elimination and a second step
given by a backward substitution. The second step requires several operations
proportional to the number of non-zero elements in the row form matrix.
However, each non-zero element has to be non-zero at the beginning of the
algorithm or originates from an operation performed in the first step.
Therefore, the second step does not add anything to the computational
complexity of the whole algorithm. For practical reasons, most sparse matrix
algorithms keep an element in the memory even if it has been modified to be
$0$. For the sake of brevity, we call a non-zero element to any element that
is or has previously been set to a non-zero value, without regard to whether
it is equal to $0$ at a given time.
The computational complexity of the first step can be analyzed as a sum of the
complexities of eliminating rows for each element partition tree node. Let us
make the following set of observations:

  1. 1.

A non-zero element in the initial matrix happens when the two basis functions
corresponding to that row and column have overlapping supports. Let us call
the graph created by considering the initial matrix to be an adjacency matrix
of a graph as an overlap graph. Two graph nodes cannot be neighbors in an
overlap graph unless the supports of their corresponding basis functions
overlap.

  2. 2.

When a row is eliminated, the new non-zero elements are created on the
intersection of columns and rows that has non-zero values in the eliminated
row or corresponding column. If we analyze the matrix as a graph, then
elimination of the row corresponding to a graph node produce edges between all
pairs of nodes that were neighbors of the node being removed.

  3. 3.

If at any given time during the forward elimination step a non-zero element
exists on the intersection of a row and a column corresponding to two basis
functions, then either those two basis functions have corresponding graph
nodes that are neighbors in the overlap graph, or that there exists a path
between those two nodes in the overlap graph that traverses only elements that
have been eliminated already.

  4. 4.

All variables corresponding to the neighboring nodes of the graph node of a
variable $x_{j}$ in the overlap graph are either:

    1. (a)

listed in one of the element partition tree nodes that are descendants of the
element partition tree node listing the variable $x_{j}$ – and those variables
are eliminated already by the time this variable is eliminated, or

    2. (b)

listed in the same element partition tree node as the variable $x_{j}$, or

    3. (c)

having the support of the corresponding basis function intersected by the
boundary of the submesh of the element partition tree node containing the
variable $x_{j}$ – those graph nodes are listed in one of the ancestors of the
element partition tree node listing the variable $x_{j}$.

Thus, in the overlap graph, there are no edges between nodes that belong to
two different element partition tree nodes that are not in an ancestor-
descendant relationship. At the same time, any path that connects a pair of
non-neighboring nodes in the overlap graph has to go through at least one
graph node corresponding to a variable that is listed in a common ancestor of
the element partition tree nodes containing the variables from that pair of
nodes.

These observations lead to a conclusion that a removal of a single row
requires no more than $(z-1)z$ subtractions, where $z$ is the amount of
variables related with graph nodes that were not removed earlier, have
overlapping support and are listed at the same tree node. Last leads to the
following recursive formula for the complexity of removing all mesh nodes at a
single tree node:

$\begin{array}[]{r @{\;} l}T(node)&=\displaystyle
T(child_{1}(node))+T(child_{2}(node))+\sum_{i=1}^{a}(b-i)(b-i+1)\vskip 3.0pt
plus 1.0pt minus 1.0pt\\\ &=\displaystyle
T(child_{1}(node))+T(child_{2}(node))+\frac{1}{3}a(a^{2}+6b+3b^{2}+3ab+3a+2)\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ &=\displaystyle
T(child_{1}(node))+T(child_{2}(node))+\mathcal{O}(ab^{2}),\end{array}$ (1)

where $a$ denotes the number of variables removed at a given tree node, and
$b$ is equal to $a$ plus the number of variables on the interface of that tree
node.
In the following sections, we denote the complexity of the removal of $a$
variables belonging to a tree node from $b$ variables having overlapping
support with the corresponding node as:

$C_{r}(a,b)=\frac{1}{3}a(a^{2}+6b+3b^{2}+3ab+3a+2)=\mathcal{O}(ab^{2}).$ (2)

##  3 Computational complexities of space–time formulations and time–marching
schemes

This section considers some possible space–time mesh refinement patterns
resulting from different space–time singularities (cf. A and B). We look at
the structure of the $d$-dimensional space–time refined mesh. Namely, we focus
on

  * 1.

space–time “edge”, resulting from a point traveling through space and time,

  * 2.

space–time “face”, resulting from an edge traveling through space and time,

  * 3.

space–time “hyperface”, resulting from a face traveling through space and
time.

We consider these space–time refined meshed in three-dimensions ($d=3$) and in
four-dimensions ($d=4$). These space–time refined $d$-dimensional meshes
correspond to the following sequences of $d-1$ dimensional refined meshes
employed by the time–marching scheme:

  * 1.

sequence of spatial meshes refined towards “points” located on the space–time
edge at particular time moments,

  * 2.

sequence of spatial meshes refined towards “edges” located on the space–time
face at a particular time moment,

  * 3.

sequence of spatial meshes refined towards “faces” located on the space–time
hyperface at a particular time moment.

This correspondence is illustrated in Figure 4 for the three-dimensional
space–time mesh with edge singularity, and the resulting sequence of two-
dimensional meshes refined to the corresponding points.

Figure 4: The refinement towards a space–time “edge” resulting from a point
traveling through space and time. This space–time refinement pattern
corresponds to a sequence of spatial meshes refined towards “points” located
on the space–time edge at particular time moments.

We assume that the time step size in the time–marching scheme is equal to the
diameter along the time axis of the smallest element in the space–time mesh.

###  3.1 Computational complexities of direct solvers for space–time
formulations and time–marching schemes

We compare the computational complexity of the direct solver executed over the
$d$ dimensional space–time domain to the computational complexity of the
direct solver executed multiple times over $d$-1 dimensional meshes within the
time–marching scheme. In Table 1 the mesh dimensions and direct solver time
complexity for space–time adaptive meshes are presented. We include the 4D
space–time uniform mesh, 4D space–time mesh refined towards the hyperface
singularity, 4D space–time mesh refined towards the face singularity, 4D
space–time mesh refined towards edge singularity, as well as 3D space–time
uniform mesh, 3D space–time mesh refined towards face singularity, and 3D
space–time mesh refined towards edge singularity. We compute sizes $N$ of all
these grids and express it by the number of refinement levels $r$.
Additionally, we provide the estimates for the direct solver execution time.
The estimations presented in Table 1 are based on theoretical results
presented in Appendix A and Appendix B, see Table 8 and 9. An interesting
observation is that the computational complexity of the direct solver executed
for the space–time mesh refined towards an edge is linear ${\cal O}(N)$.
Next, in Table 2 we evaluate the computational complexity for time–marching
scheme. For each $d$-dimensional space–time refined toward the
$q$-singularity, we construct a sequence of $2^{r}$ meshes. Each of the grids
from the sequence is $d-1$ dimensional, refined towards $q-1$ singularity.
Such a sequence of grids provides the solution with similar accuracy to the
space–time grid. We estimate the dimensions $n$ of the spatial meshes from the
sequence and the computational complexity of the direct solver executing
$2^{r}$ times, once for each grid from the sequence. The estimations presented
in Table 2 are based on theoretical results presented in Appendix A and
Appendix B, see Table 8 and 9. An interesting observation is that the
computational complexity of the time–marching scheme corresponding to the
space–time mesh refined towards an edge is $r$ times higher than the
space–time complexity. However, the time marching scheme is cheaper for all
other space–time grids.

space–time mesh | space–time mesh size | space–time mesh
---|---|---
| $N$ | direct solver complexity
4D uniform | ${\cal O}(2^{4r})$ | ${\cal O}(N^{\frac{9}{4}})={\cal O}(2^{9r})$
4D hyperface | ${\cal O}(8^{r})={\cal O}(2^{3r})$ | ${\cal O}(N^{2})={\cal O}(2^{6r})$
4D face | ${\cal O}(4^{r})={\cal O}(2^{2r})$ | ${\cal O}(N^{\frac{3}{2}})={\cal O}(2^{3r})$
4D edge | ${\cal O}(2^{r})$ | ${\cal O}(N)={\cal O}(2^{r})$
3D uniform | ${\cal O}(2^{3r})$ | ${\cal O}(N^{2})={\cal O}(2^{6r})$
3D face | ${\cal O}(4^{r})={\cal O}(2^{2r})$ | ${\cal O}(N^{\frac{3}{2}})={\cal O}(2^{3r})$
3D edge | ${\cal O}(2^{r})$ | ${\cal O}(N)={\cal O}(2^{r})$
Table 1: space–time formulations, mesh dimensions and direct solver complexities space–time mesh | Sequence of | Single spatial | Total direct solver
---|---|---|---
| spatial meshes | mesh size | complexity for a sequence
4D uniform |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D uniform | $n={\cal O}(2^{3r})$ | $2^{r}\times{\cal O}(n^{2})=$
$N={\cal O}(2^{4r})$ |  |  | ${\cal O}(2^{7r})$
4D hyperface |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D face | $n={\cal O}(2^{2r})$ | $2^{r}\times{\cal O}(n^{\frac{3}{2}})=$
$N={\cal O}(2^{3r})$ |  |  | ${\cal O}(2^{4r})$
4D face |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D edge | $n={\cal O}(2^{r})$ | $2^{r}\times{\cal O}(n)=$
$N={\cal O}(2^{2r})$ |  |  | ${\cal O}(2^{2r})$
4D edge |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D point | $n={\cal O}(r)$ | $2^{r}\times{\cal O}(r)=$
$N={\cal O}(2^{r})$ |  |  | ${\cal O}(r2^{r})$
3D uniform |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 2D uniform | $n={\cal O}(2^{2r})$ | $2^{r}\times{\cal O}(n^{\frac{3}{2}})=$
$N={\cal O}(2^{3r})$ |  |  | ${\cal O}(2^{4r})$
3D face |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 2D edge | $n={\cal O}(2^{r})$ | $2^{r}\times{\cal O}(2^{r})=$
$N={\cal O}(2^{2r})$ |  |  | ${\cal O}(2^{2r})$
3D edge |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 2D point | $n={\cal O}(r)$ | $2^{r}\times{\cal O}(r)=$
$N={\cal O}(2^{r})$ |  |  | ${\cal O}(r2^{r})$
Table 2: time–marching schemes corresponding to space–time formulations

###  3.2 Computational complexities of iterative solvers for space–time
formulations and time–marching schemes

Finally, we focus on the iterative solver. The computational complexities of
the iterative solver for both space–time grids and the time–marching scheme
are estimated as in Table 3. The computational complexity for the space–time
grid is equal to ${\cal O}(N_{iter}\times N)$ where $N$ is the size of the
space–time grid, and $N_{iter}$ is the number of iterations of the iterative
solver for the space–time grid. The computational complexity for the
space–time grid is equal to ${\cal O}(2^{r}\times n_{iter}\times n)$ where
$2^{r}$ is the number of steps of the time–marching scheme, $n$ is the size of
the spatial grid from the sequence, and $n_{iter}$ is the number of iterations
on the spatial grid. It is obvious to assume that $N_{iter}>>n_{iter}$, but
the exact numbers of iterations are problem-dependent. The computational
complexity of the iterative solver for space–time grids is higher than the
computational complexity of the iterative solver for the time–marching scheme.
An interesting case is the space–time grid refined to the space–time edge,
following the trajectory of the point object. In this case, the complexity of
iterative solver for space–time grid is ${\cal O}(N_{iter}\times N)={\cal
O}(N_{iter}\times 2^{r})$ while the complexity of the iterative solver for
time marching grids is ${\cal O}(2^{r}\times n_{iter}\times n)={\cal
O}(n_{iter}\times r2^{r})$. The iterative solver for space–time grid is
cheaper than for the time marching scheme, if $N_{iter}<r\times n_{iter}$.

space–time | Sequence of | space–time | time–marching scheme
---|---|---|---
mesh | spatial meshes | iterative | iterative
|  | solver complexity | solver complexity
4D uniform |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D uniform |  | ${\cal O}(2^{r}\times n_{iter}\times 2^{3r})=$
$N={\cal O}(2^{4r})$ |  | ${\cal O}(N_{iter}\times 2^{4r})$ | ${\cal O}(n_{iter}\times 2^{4r})$
4D hyperface |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D face |  | ${\cal O}(2^{r}\times n_{iter}\times 2^{2r})=$
$N={\cal O}(2^{3r})$ |  | ${\cal O}(N_{iter}\times 2^{3r})$ | ${\cal O}(n_{iter}\times 2^{3r})$
4D face |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D edge |  | ${\cal O}(2^{r}\times n_{iter}\times 2^{r})=$
$N={\cal O}(2^{2r})$ |  | ${\cal O}(N_{iter}\times 2^{2r})$ | ${\cal O}(n_{iter}\times 2^{2r})$
4D edge |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 3D point |  | ${\cal O}(2^{r}\times n_{iter}\times r)=$
$N={\cal O}(2^{r})$ |  | ${\cal O}(N_{iter}\times 2^{r})$ | ${\cal O}(n_{iter}\times r2^{r})$
3D uniform |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 2D uniform |  | ${\cal O}(2^{r}\times n_{iter}\times 2^{2r})=$
$N={\cal O}(2^{3r})$ |  | ${\cal O}(N_{iter}\times 2^{3r})$ | ${\cal O}(n_{iter}\times 2^{3r})$
3D face |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 2D edge |  | ${\cal O}(2^{r}\times n_{iter}\times 2^{r})=$
$N={\cal O}(2^{2r})$ |  | ${\cal O}(N_{iter}\times 2^{2r})$ | ${\cal O}(n_{iter}\times 2^{2r})$
3D edge |  ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}\times}$ 2D point |  | ${\cal O}(2^{r}\times n_{iter}\times r)=$
$N={\cal O}(2^{r})$ |  | ${\cal O}(N_{iter}\times 2^{r})$ | ${\cal O}(n_{iter}\times r2^{r})$
Table 3: Comparison of complexities of iterative solvers for space–time
formulations and time–marching schemes.

Additionally, we compare the computational complexity of the iterative solver
executed over the $d$ dimensional space–time domain to the computational
complexity of the iterative solver executed multiple times over the $d$-1
dimensional meshes during the time–marching scheme.

###  3.3 Impact of polynomial order of approximation

Notice that we have not included the $p$-factor in the computational
complexity estimates for the arbitrary shape of singularity. However, we can
easily estimate the computational complexity of the static condensation
performed at the beginning of the computations with a higher-order finite
element method with hierarchical basis functions.
The computational complexity of the static condensation over a single element
is equal to ${\cal O}(p^{3d})$. We are eliminating the degrees of freedom from
the interior of the element; we have $(p-1)^{d}$ degrees of freedom there; the
matrix is dense, so the complexity of elimination is the cube of the number of
degrees of freedom. In other words, the complexity of static condensation for
a grid with $N$ elements is ${\cal O}(Np^{3d})$. With this observation in
mind, we can estimate the complexities of static condensations for space–time
and time–marching grids. It is illustrated in Table 4. The complexity of
static condensations is higher for the space–time mesh. However, for the
space–time formulation, we can always execute all the elemental computations
at the same time, in parallel, while for the time–marching scheme, we can only
run the static condensation over a single time-step mesh in parallel. In this
sense, the space–time method can outperform the time–marching scheme.

space–time | Sequence of | space–time | time–marching scheme
---|---|---|---
mesh | spatial meshes | complexity of static | complexity of static
|  | condensation | condensation
4D uniform |  $2^{r}$ steps $\times$ |  | ($2^{r}$ steps)
$N={\cal O}(2^{4r})$ | 3D uniform | ${\cal O}(p^{12}2^{4r})$ | ${\cal O}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}}\times p^{9}2^{3r})$
4D hyperface |  $2^{r}$ steps $\times$ |  | ($2^{r}$ steps)
$N={\cal O}(2^{3r})$ | 3D face | ${\cal O}(p^{12}2^{3r})$ | ${\cal O}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}}\times p^{9}2^{2r})$
4D face |  $2^{r}$ steps $\times$ |  | ($2^{r}$ steps)
$N={\cal O}(2^{2r})$ | 3D edge | ${\cal O}(p^{12}2^{2r})$ | ${\cal O}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}}\times p^{9}2^{r})$
4D edge |  $2^{r}$ steps $\times$ |  | ($2^{r}$ steps)
$N={\cal O}(2^{r})$ | 3D point | ${\cal O}(p^{12}2^{r})$ | ${\cal O}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}}\times p^{9}r)$
3D uniform |  $2^{r}$ steps $\times$ |  | ($2^{r}$ steps)
$N={\cal O}(2^{3r})$ | 2D uniform | ${\cal O}(p^{9}2^{3r})$ | ${\cal O}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}}\times p^{6}2^{2r})$
3D face |  $2^{r}$ steps $\times$ |  | ($2^{r}$ steps)
$N={\cal O}(2^{2r})$ | 2D edge | ${\cal O}(p^{9}2^{2r})$ | ${\cal O}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}}\times p^{6}2^{r})$
3D edge |  $2^{r}$ steps $\times$ |  | ($2^{r}$ steps)
$N={\cal O}(2^{r})$ | 2D point | ${\cal O}(p^{9}2^{r})$ | ${\cal O}({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}2^{r}}\times p^{6}r)$
Table 4: Comparison of complexities of static condensations for space–time
formulations and time–marching schemes.

##  4 Numerical results

In this section we provide numerical experiments for verification of the
theoretical findings summarized in Tables 1-3. We employ our Octave codes to
generate the structure of $d$-dimensional computational grids with
$q$-dimensional singularity, generated using $r$ refinement levels. For
simplicity of implementation, we generate our matrices by looking at relations
between finite elements. Rows and columns in matrices correspond to finite
elements. Non-zero entries in a row mean that two elements, one related to the
row, and one related to the column, are adjacent through a $d-1$ dimensional
face. This way of generating matrices influences the computational complexity
constant, ignoring the polynomial order of approximation, but the dependence
on $N$ equal here the number of elements is of the same order as if we include
all the relations of the basis functions.
We run experiments using Octave on a Linux cluster node equipped with 2.4GHz
processor with 64 GB of RAM. We cannot factorize more than eight refinements
for the face singularity in four dimensions and nine refinement levels for the
face singularity in three dimensions because of a lack of memory during the
factorization process. The comparisons of execution times for four-dimensional
face singularity versus a sequence of three-dimensional edge singularities are
presented in Table 5. The comparisons of execution times for three-dimensional
face singularity versus a sequence of two-dimensional edge singularities are
presented in Table 6. We employ AMD ordering and multi-frontal solver as
implemented in the Octave, e.g., for the face singularity with nine refinement
levels in three dimensions, we run our matrix generation script:

F9_3=Face(9,3);

N = 87381

We compute the AMD permutation

p=amd(F9_3);

and we plug it into the LU factorization, measuring the execution time

tic; lu(F9_3(p,p)); toc

Elapsed time is 46.2016 seconds.

We know this ordering is different from the one proposed in our paper.
Nevertheless, the results show up to one order of magnitude times faster
execution times of time–marching schemes with edge singularities than one call
for the space–time domain with the face singularity.
The comparisons of execution times for four-dimensional edge singularity
versus a sequence of three-dimensional point singularities are presented in
Table 7. We can perform 15 refinements over the space–time mesh this time with
the Octave solver. The space–time mesh with a multi-frontal solver is, in this
case, faster than the time marching scheme, up to the 13 refinement level.
With 14 or 15 refinements, the complexity of processing the space–time mesh is
higher than the complexity of processing the time–marching solver.

space–time | r | N | space–time | Sequence of | time–marching scheme
---|---|---|---|---|---
mesh |  |  | solver time [s] | spatial meshes | solver time [s]
4D face | 6 | 4095 | 0.56 |  $64\times$ 3D edge | $64\times 0.0037=0.0236$
4D face | 7 | 16383 | 13.79 |  $128\times$ 3D edge | $128\times 0.024=3.072$
4D face | 8 | 65535 | 234 |  $256\times$ 3D edge | $256\times 0.063=16.128$

Table 5: Execution times for four dimensional mesh refined towards face
singularity and the corresponding sequence of three-dimensional meshes refined
towards edge singularities.

.

space–time | r | N | space–time | Sequence of | time–marching scheme
---|---|---|---|---|---
mesh |  |  | solver time [s] | spatial meshes | solver time [s]
3D face | 6 | 1365 | 0.044 |  $64\times$ 2D edge | $64\times 0.0008=0.05$
3D face | 7 | 5461 | 0.28 |  $128\times$ 2D edge | $128\times 0.001=0.128$
3D face | 8 | 21845 | 2.55 |  $256\times$ 2D edge | $256\times 0.003=0.76$
3D face | 9 | 87381 | 46.20 |  $512\times$ 2D edge | $512\times 0.0074=3.78$
Table 6: Execution times for three dimensional mesh refined towards face
singularity and the corresponding sequence of two-dimensional meshes refined
towards edge singularities.

space–time | r | N | space–time | Sequence of | time–marching scheme
---|---|---|---|---|---
mesh |  |  | solver time [s] | spatial meshes | solver time [s]
4D edge | 6 | 434 | 0.008 |  $64\times$ 3D point | $64\times 0.0004=0.0256$
4D edge | 7 | 882 | 0.026 |  $128\times$ 3D point | $128\times 0004=0.0512$
4D edge | 8 | 1778 | 0.069 |  $256\times$ 3D point | $256\times 0.0005=0.128$
4D edge | 9 | 3570 | 0.15 |  $512\times$ 3D point | $512\times 0.0006=0.3$
4D edge | 10 | 7154 | 0.38 |  $1024\times$ 3D point | $1024\times 0.0006=0.61$
4D edge | 11 | 14332 | 1.01 |  $2048\times$ 3D point | $2048\times 0.0007=1.43$
4D edge | 12 | 28665 | 2.45 |  $4096\times$ 3D point | $4096\times 0.0008=3.27$
4D edge | 13 | 57330 | 5.52 |  $8192\times$ 3D point | $8192\times 0.0007=5.73$
4D edge | 14 | 114674 | 13.59 |  $16384\times$ 3D point | $16384\times 0.0007=11.46$
4D edge | 15 | 229362 | 38.96 |  $32768\times$ 3D point | $32768\times 0.0009=29.49$

Table 7: Execution times for four dimensional mesh refined towards edge
singularity and the corresponding sequence of three-dimensional meshes refined
towards point singularities.

.

##  5 Conclusions

To estimate the computational complexity of full space–time formulations and
time–marching schemes for hypercubic elements, we simulate several possible
scenarios for the resulting matrices when considering adaptivity toward
singularities. In particular, we consider refinements towards the point, edge,
face, and hyperface singularities over space–time mesh. In our idealized case,
we refine all the elements that contain the prescribed point, edge, face or
hyperface singularity. Thus, we obtain several representative refined
$d$-dimensional computational meshes, where we assume that we perform
refinements towards $q$-dimensional manifold ($q<d$) representing the
singularities. For each of these representative meshes, we estimate the number
of degrees of freedom (second column in Table 1), the computational complexity
of the multi-frontal solver (third column in Table 1), the computational
complexity of the iterative solver (third column in Table 3), and the
computational complexity of the static condensation (third column in Table 4).
On the other hand, we generated a sequence of refined $d-1$ dimensional
computational meshes, representing the “cross-sections” of the $d$ dimensional
space–time mesh. In this case, we performed refinements towards
$q-1$-dimensional manifold representing the cross-section of the
$q$-dimensional singularity. Our theoretical estimations and the numerical
experiments imply that the time–marching scheme is competitive only in this
idealized case when the number of time steps is equal to the element size in
the time dimension. This concerns the computational complexity of the multi-
frontal solver (third column in Table 1 versus the fourth column in Table 2),
the iterative solver (third column in Table 3 and fourth column in Table 3),
and the static condensation (third and fourth column in Table 4). We also
present numerical experiments, confirming the predicted theoretical behaviors.

We understand that our assumptions are the best possible idealistic scenarios.
In the real life applications, the computational complexity of space–time
formulation is more competitive due to

  * 1.

Increased number of time steps in higher-order and accuracy time marching
schemes, where the time-step size can be actually smaller than the temporal
dimension of the smallest elements in the space–time mesh.

  * 2.

Extensive parallelization of the computational process, where for example the
static condensation for the space–time formulation can be performed fully in
parallel, and the static condensation of the time–marching scheme has the
limitation of the single size of the time-step mesh.

  * 3.

Sequential nature of the time–marching scheme, where the iterative solver has
to be executed in a sequence for each time-step mesh, and it cannot be
parallelized once for the entire computational space–time mesh.

  * 4.

The cost of generation of the refined computational meshes is in general
ignored in our estimations (it is assumed to be linear), while in general, it
is an iterative procedure that requires several solves, and in the space–time
setup it can be performed once for the entire mesh, but in the time–marching
scheme it has to be performed for each time-step mesh.

Nevertheless, our estimate constitutes the lower bounding case of the
computational complexities for both space–time and time–marching schemes.

In our future work, we plan to extend our computational complexity estimates
into parallel distributed memory [24] and shared-memory machines [25],
considering direct and iterative solvers. Our future work will also involve
the computational complexity analysis of the hierarchical matrices solvers
(H-matrices) in the space–time setup. In general, the complexity of processing
these matrices is proportional to the problem size $N$, multiplied by an
additional factor that grows with the dimension size.

## Acknowledgments

This work from the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreement No 777778
(MATHROCKS). The work of SR has also been supported by the Chilean grant ANID
Fondecyt No 3210009.

##  Appendix A Computational complexity for h-adapted meshes towards a point
singularity

This Appendix analyzes the time complexity of the direct solver being run on
hierarchical meshes adapted toward a point singularity.

###  A.1 h-adapted mesh towards a point singularity

We define as a point singularity mesh a mesh refined hierarchically towards a
single point. For example, let us consider a point singularity in a
$d$-dimensional space refined towards a point $Q$ until some arbitrary
refinement level $r$. We denote such mesh as $\mathcal{G}_{Q}^{d}(r)$. The
singularity point can be placed either inside an element or on an element’s
boundary. In particular, the singularity point can be placed on the whole mesh
boundary. Figure 5(b) illustrates examples of meshes with point singularities.

Figure 5: Examples of point singularities.

(a) Examples of $h$-adaptive singularity meshes in 2D. Green dots denote the
singularity position.

(b) Examples of $h$-adaptive singularity meshes in 3D. Elements containing
singularities marked in red.

Figure 6 illustrates a process of building such a mesh using Algorithm 1.

Figure 6: A $G_{q}^{2}(4)$ mesh construction for a point $q$ close to the
middle of the mesh.

###  A.2 Analysis: point singularity placed on the boundary of the mesh

As an example of a point singularity mesh, we analyze a relatively simple case
with a singularity point $Q$ placed in one of the corners of the mesh. We show
a way to generate an ordering that results in the linear complexity of the
solver. Figure 7 shows examples of $h$-adaptive meshes of this type in two and
three dimensions.

Figure 7: Meshes in 2D and 3D with corner point singularity marked in red.

For this kind of meshes, it is not required to perform an extra step to ensure
that the $1$-irregularity rule is met. Therefore, the number of mesh elements
$N_{e}$ is described by the following formula:

$N_{e}(G_{Q}^{d}(r))=\begin{cases}1&\text{if }r=0,\\\
N_{e}(G_{Q}^{d}(r-1))+2^{d}-1&\text{if }r>0,\end{cases}$ (3)

where $r$ is the refinement level of the mesh, and $d$ is the dimensionality
of the mesh. Formula (3) is equivalent to:

$\displaystyle N_{e}(G_{Q}^{d}(r))$
$\displaystyle=r(2^{d}-1)+1=\mathcal{O}(2^{d}r).$ (4)

In other words, the number of elements in a corner point singularity mesh
grows linearly with the refinement level $r$.
After the first refinement, the mesh has $(2p+1)^{d}$ variables. Further
refinement level adds a layer of $(2^{d}-1)$ elements, and for each new
element, $p^{d}$ variabes are created. The following formula describes the
number of the basis functions on such a mesh for $r\geq 1$:

$\displaystyle
N_{v}(G_{Q}^{d}(r))={(2p+1)}^{d}+(r-1)(2^{d}-1)p^{d}=\mathcal{O}(2^{d}rp^{d}).$
(5)

If $p$ is constant, the number of variables is linearly proportional to the
refinement level $r$. Figure 8 visually explains how the number of basis
functions grows.

Figure 8: Each layer of refinement adds $2^{d}-1$ elements and each element
adds $p^{d}$ variables.

####  A.2.1 Time complexity

We use the element partition tree method of ordering generation to calculate
the time complexity. The element partition tree can be built by recursively
removing the layer of $2^{d}-1$ least refined (or largest) elements. Figure 9
illustrates the structure of such an element partition tree.

Figure 9: Corner point singularity element partition tree in 2D. Green nodes
are removed at given partition tree node, gray nodes are removed deeper in the
tree, and red nodes are removed in the ancestor of the tree node.

Let us denote each layer by $s_{i}$, where $i$ is the number of the shell,
counting from the most refined one, $s_{0}$ being the particular case of the
single element in the corner. It is easy to see that the whole tree is
recursive and that each layer $s_{i}$, with $i>0$, contains $2^{d}-1$
elements.
During the solution process at each layer $s_{i}$ ($i\geq 1$), we first remove
all the rows corresponding to the nodes with support entirely within $s_{i}$
(denoted by $\mathit{If}(s_{i})$). Then, we remove all rows that correspond to
nodes on the interface between $s_{i-1}$ and $s_{i}$ (denoted
$\mathit{If}(s_{i-1},s_{i})$). This results in the following complexity
formula for a single layer (excluding layers $s_{0}$, $s_{1}$ and $s_{r}$):

$\begin{array}[]{r @{\;}
l}T(s_{i})&=C_{r}(|\mathit{If}(s_{i})|,|\mathit{If}(s_{i})|+|\mathit{If}(s_{i-1},s_{i})|+|\mathit{If}(s_{i},s_{i}+1)|)\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
&+C_{r}(|\mathit{If}(s_{i-1},s_{i})|,|\mathit{If}(s_{i-1},s_{i})|+|\mathit{If}(s_{i},s_{i}+1)|)\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
&=|\mathit{If}(s_{i})|\left(|\mathit{If}(s_{i})|+2|\mathit{If}(s_{i-1},s_{i})|\right)^{2}+4|\mathit{If}(s_{i-1},s_{i})|^{3}.\end{array}$
(6)

To calculate the number of variables on the interface between layers, we
realize that this number stays the same from when the given layer is refined.
In other words, the number of variables on the interface between any layers is
equal to the number of variables on the interface between layers $s_{0}$ and
$s_{1}$:

$|\mathit{If}(s_{i-1},s_{i})|={(p+1)}^{d}-p^{d}.$ (7)

Similarly, we can calculate the number of ”internal” variables, which is equal
to the number of variables in a mesh of refinement level 1 ($G^{d}(1)$),
excluding the variables on the (hyper)faces not touching the corner with
singularity and minus the number of variables in a mesh of refinement level 0
($G^{d}(0)$):

$|\mathit{If}(s_{i})|=(2^{d}-1)p^{d}-|\mathit{If}(s_{i-1},s_{i})|=(2^{d}-1)p^{d}+p^{d}-(p+1)^{d}=(2p)^{d}-(p+1)^{d}$
(8)

Therefore, the complexity of the elimination of a single layer $s_{i}$, where
$1<i<r$, goes as follows:

$\begin{array}[]{r @{\;}
l}T(s_{i})&=|\mathit{If}(s_{i})|(|\mathit{If}(s_{i})|+2|\mathit{If}(s_{i-1},s_{i})|)^{2}+4|\mathit{If}(s_{i-1},s_{i})|^{3}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
&=((2p)^{d}-(p+1)^{d})((2p)^{d}-{(p+1)}^{d}+{(p+1)}^{d}-p^{d})^{2}+4({(p+1)}^{d}-p^{d})^{3}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
&=((2p)^{d}-(p+1)^{d})((2^{d}-1)p^{d})^{2}+4({(p+1)}^{d}-p^{d})^{3}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
&\leq(2p)^{d}((2p)^{d})^{2}+4(p+1)^{3d}=(2p)^{3d}+4(p+1)^{3d}=\mathcal{O}{\left((2p)^{3d}\right)}.\end{array}$
(9)

The computational complexity of the whole mesh is then given as follows:

$\displaystyle
T(G^{d}(r))=\mathcal{O}(r(2p)^{3d})=\mathcal{O}\left((2p)^{2d}\cdot
2^{d}rp^{d}\right)\mathcal{O}\left(2^{2d}p^{2d}N\right).$ (10)

With the assumption that the parameter $p$ is constant, and considering that
dimension $d$ is constant for a given problem, we deduce that the time
complexity of the whole algorithm is linear $\mathcal{O}(N)$.

###  A.3 Quasi-optimal h-adapted mesh towards a point singularity

To generalize the analysis above, we present a set of properties of the mesh
that guarantee that it is possible to create an ordering that leads to linear
solver execution time. We call a mesh that fulfills those properties quasi-
optimal point singularity mesh. The properties of quasi-optimal h-adapted mesh
towards a point singularity are the following:

  1. 1.

The mesh can be split into no more than $Kr$ consecutive layers, where $r$ is
the refinement level of the mesh, and $K$ is some constant.

  2. 2.

Each basis function can be assigned to one of the layers so that a pair of
basis functions have non-overlapping supports if they are more than $M$ shells
apart, where $M$ is some constant.

  3. 3.

Each layer has no more than $L$ basis functions assigned, where $L$ is some
constant;

If those properties are met, it is possible to create an ordering in which the
nodes are removed layer by layer. Removal of each successive layer requires no
more than $C_{r}(L,L(M+1))=O(L^{3}M^{2})$ subtractions, and the whole solver
requires no more than $\mathcal{O}(KrL^{3}M^{2})$ operations. The constructed
mesh has $N=\mathcal{O}(KrL)$ nodes, so the time complexity of the solver is
stated as $\mathcal{O}(N(LM)^{2})$. As $K$, $L$, and $M$ are constants, the
resulting solution complexity is $\mathcal{O}(N)$. That is, it grows linearly
with the number of nodes (which, in turn, increases linearly with the
refinement level).

###  A.4 Quasi-optimality of $h$-adaptive point singularity meshes

A final step of the proof presented in this section is to show that any
$h$-adaptive point singularity mesh created with Algorithm 1 is quasi-optimal.
It is easy to notice that:

  1. 1.

The elements can be grouped into $R$ layers by their refinement level, and
each basis function can be assigned to the layer of one of its elements.

  2. 2.

Because of the $1$-irregularity rule, two elements sharing a vertex cannot
differ by more than two refinement levels, so a single basis function cannot
span over two elements with a difference of more than two refinement levels.

  3. 3.

There are no more than $4^{d}$ elements of each refinement level. At the same
time, each element has no more than $(p+1)^{d}$ basis functions – both values
are constant.

Those observations lead to the conclusion that any point singularity
$h$-adaptive mesh is quasi-optimal:

  1. 1.

There are exactly $R$ layers.

  2. 2.

Nodes that are more than two layers apart never overlap.

  3. 3.

There are no more than $4^{d}(p+1)^{d}$ nodes at each layer.

Extending those observations and seeing that a mesh is refined towards more
than one is straightforward. Still, a finite number of point singularities is
quasi-optimal in the same way (each layer potentially has the number of basis
functions multiplied by the number of singularities).

##  Appendix B Computational complexity for h-adapted meshes towards a multi-
dimensional singularity

This Appendix analyzes the time complexity of direct solvers being run over
adaptive meshes refined around singularities of simple shapes of higher
dimensionality. In particular, we analyze meshes with singularities in the
form of lines, planes/faces, and hyperplanes/hyperfaces with three or more
dimensions, depending on the dimensionality of the space. In this case, for
the simplicity of derivation, we ignore the polynomial order $p$ factor.

###  B.1 Structure of the mesh with singularity

This section analyzes $d$-dimensional meshes refined towards $q$-dimensional
singularities. The mesh is considered refined until the refinement level $r$
towards that singularity if all elements overlapping any section containing
the singularity have been refined. Such refinement can again be achieved using
Algorithm 1. A process of such refinement is shown in Figure 10(a) and 10(b),
respectively for edge and face singularity in 3-dimensional space. It is
crucial to notice that regular meshes can also be analyzed as mesh refined
toward $d$-dimensional singularity in $d$-dimensional space.

(a) Example of edge singularity mesh construction

(b) Example of face singularity mesh construction

Figure 10: Examples of mesh singularities in 3D

###  B.2 Analysis: singularity placed on the boundary of the mesh

For the sake of simplicity, we start by analyzing the basic case of the
singularity placed on the boundary of the full $h$-adapted mesh. Let us denote
the dimensionality of the mesh as $d$, the dimensionality of the singularity
as $q$, and the singularity itself as $S_{q}$. The mesh is denoted as
$\mathcal{G}_{S_{q}}^{d}(r)$, where $r$ is its refinement level.

####  B.2.1 Properties of the mesh

The number of elements in such a singularity mesh can be calculated using the
following recursive formula:

$\begin{array}[]{r @{\;} l}N_{e}(\mathcal{G}_{S_{q}}^{d}(0))&=1\vskip 3.0pt
plus 1.0pt minus 1.0pt\\\
N_{e}(\mathcal{G}_{S_{q}}^{d}(r))&=N_{v}(\mathcal{G}_{S_{q}}^{d}(r-1))+2^{q(r-1)}(2^{d}-1)\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
N_{e}(\mathcal{G}_{S_{q}}^{d}(r))&=\begin{cases}(2^{d}-1)r+1&\text{if
}q=0\vskip 3.0pt plus 1.0pt minus 1.0pt\\\
\frac{(2^{d}-1)2^{rq}-(2^{d}-2^{q})}{2^{q}-1}&\text{if }q\geq
1\end{cases}\end{array}$ (11)

The formula expands to the values given in Table 8.

| Point ($q=0$) | Edge ($q=1$) | Face ($q=2$) | Hyperface ($q=3$)
---|---|---|---|---
1-D | $r+1$ | $2^{r}$ |  |
2-D | $3r+1$ | $3(2^{r})-2$ | $4^{r}$ |
3-D | $7r+1$ | $7(2^{r})-6$ | $\tfrac{1}{3}(7(4^{r})-4)$ | $8^{r}$
4-D | $15r+1$ | $15(2^{r})-14$ | $\tfrac{1}{3}(15(4^{r})-12)$ | $\tfrac{1}{7}(15(8^{r})-8)$
$d$-D | $(2^{d}-1)r+1$ | $(2^{d}-1)2^{r}-(2^{d}-2)$ | $\tfrac{1}{3}((2^{d}-1)4^{r}-(2^{d}-4))$ | $\tfrac{1}{7}((2^{d}-1)8^{r}-(2^{d}-8))$
Table 8: Number of elements for d–dimensional h–adapted meshes towards
g-dimensional singularities

For $q\geq 1$, the number of elements is approximated by the following lower
and upper bounds:

$2^{(d-1)}(2^{q})^{r}\leq N_{e}(\mathcal{G}_{S_{q}}^{d}(r))\leq
2^{d}(2^{q})^{r}.$ (12)

The number of variables is estimated to be between $p^{d}$ and $(p+1)^{d}$ per
element. Leading us to the following approximation for $q\geq 1$:

$2^{(d-1)}p^{d}(2^{q})^{r}\leq N_{v}(\mathcal{G}_{S_{q}}^{d}(r))\leq
2^{d}(p+1)^{d}(2^{q})^{r}.$ (13)

For set $d$, $q\geq 1$ and $p$, the approximations lead to the following
formulas:

$\displaystyle N_{e}(\mathcal{G}_{S_{q}}^{d}(r))$
$\displaystyle=\mathcal{O}((2^{q})^{r}),$ (14) $\displaystyle
N_{v}(\mathcal{G}_{S_{q}}^{d}(r))$ $\displaystyle=\mathcal{O}((2^{q})^{r}).$
(15)

In other words, both the number of elements and variables grow proportionally
to $\mathcal{O}(2^{qr})$ and this growth speed (understood in terms of
$\mathcal{O}$-notation) depends only on the dimensionality of the singularity
$q$, not on the dimensionality of the mesh $d$.

####  B.2.2 Time complexity of a solution with singularity built on mesh
boundary

To analyze the time complexity of the solver, we can again use the element
partition tree approach. The element partition tree is built using the
following recursive procedure:

  1. 1.

Create a root node of the element partition tree and attach all the elements
to that node.

  2. 2.

If there is just one element, finish the procedure, and the root node is the
sole node of the returned tree.

  3. 3.

Create a child node of the root node containing all the least refined
elements.

  4. 4.

Divide the remaining elements by $q$ parallel planes perpendicular to the
singularity and parallel to the boundaries of the mesh (let us denote those
planes as dividing planes), crossing the midpoint of the singularity. This
refinement creates $2^{q}$ sub meshes.

  5. 5.

For each sub mesh generated above, run this procedure recursively and attach
the resulting trees as subtrees of the second child node.

  6. 6.

Finish the procedure and return the tree stemming from the root node.

Figure 11 shows an example of such an element partition tree. We denote by
$s_{i}$ the element partition tree nodes from the refinement level $i$.
Even though the element partition tree nodes on each level have an analogous
set of elements, the order of elimination differs slightly for the tree nodes
that contain elements on the boundary of the mesh other than the boundary
containing the singularity. To simplify the analysis, we modify the order of
elimination slightly so that those tree nodes behave similarly to the others:
the variables corresponding to the basis functions on the boundary of the mesh
are eliminated at the root node $s_{r}$. This change increases the computation
time slightly. However, it has no impact on the time complexity.

Figure 11: Example element partition tree for one-dimensional boundary
singularity in 2-D.

To calculate the computational complexity of the solver using the ordering
generated from that element partition tree, we need to know two values for
each element partition tree node:

  * 1.

The number of variables removed in that element partition tree node $n_{r}$:
For the $s_{i}$ nodes, this number is proportional to the number of elements
of that node that are touching the dividing planes that are used to divide the
sub mesh further.

  * 2.

The total number of variables with support over the elements in this subtree
$n_{e}$: For the $s_{i}$ nodes, this number is proportional to the number of
elements on the dividing planes of the ancestral tree nodes.

It is not difficult to see that the cross-section of the mesh, through the
dividing planes, behaves as a mesh in the space of one less dimension built
over a singularity of one-dimensionality less than the original one. Thus, for
$q>1$, the following equalities hold true:

$\displaystyle n_{r}(s_{i})$
$\displaystyle=\mathcal{O}(2^{(q-1)r}q)=\mathcal{O}(2^{(q-1)r}),$ (16)
$\displaystyle n_{e}(s_{i})$
$\displaystyle=\mathcal{O}(2^{(q-1)r}q)=\mathcal{O}(2^{(q-1)r}).$ (17)

All the remaining nodes have a $\mathcal{O}(1)$ number of elements. Thanks to
those observation, we can calculate the time complexity of running the solver
for $q>1$ using the following equation:

$\begin{array}[]{r @{\;} l}T(s_{0})&=\mathcal{O}(1),\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\ T(s_{r})&=\mathcal{O}((2^{3(q-1)})^{r}),\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\
T(\mathcal{G}_{S_{q}}^{d}(r))&=\mathcal{O}((2^{3(q-1)})^{r}),\end{array}$ (18)

where the second equality in Equation (18) follows from the fact that:

$\begin{array}[]{r @{\;} l}T(s_{r})&=\displaystyle
2^{q}T(s_{r-1})+\mathcal{O}(C_{e}(n_{r},n_{r}+n_{e}))+\mathcal{O}(1)=2^{q}T(s_{r-1})+\mathcal{O}(2^{3(q-1)r})\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
&=\displaystyle\mathcal{O}\left(\sum_{h=0}^{r}2^{q(r-h)}\cdot
2^{3(q-1)h}\right)=\mathcal{O}\left(\sum_{h=0}^{r}(2^{q(r-h)+3(q-1)h})\right)=\mathcal{O}\left(\sum_{h=0}^{r}2^{q(r+2h)-3h}\right).\end{array}$

For $q=1$, analogous calculations give us the estimates:

$\displaystyle n_{r}(s_{i})$ $\displaystyle=\mathcal{O}(r),$ (19)
$\displaystyle n_{e}(s_{i})$ $\displaystyle=\mathcal{O}(r).$ (20)

Moreover, the time complexity of running the solver follows this equations:

$\begin{array}[]{r @{\;} l}T(s_{0})&=\mathcal{O}(1),\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\ T(s_{r})&=\mathcal{O}(2^{r}),\vskip 3.0pt plus 1.0pt minus
1.0pt\\\ T(\mathcal{G}_{S_{1}}^{d}(r))&=\mathcal{O}(2^{r}),\end{array}$ (21)

where the second equality in Equation (21) follows from:

$\begin{array}[]{r @{\;} l}T(s_{r})&=\displaystyle
2T(s_{r-1})+\mathcal{O}(C_{e}(n_{r},n_{r}+n_{e}))+\mathcal{O}(1)=2T(s_{r-1})+\mathcal{O}(r^{3})\vskip
3.0pt plus 1.0pt minus 1.0pt\\\
&=\displaystyle\mathcal{O}\left(\sum_{h=0}^{r}2^{r-h}\cdot
h^{3}\right)=\mathcal{O}(2^{r}\cdot 1^{3}+2^{r-1}\cdot
2^{3}+\dots+2^{1}\cdot(r-1)^{3}+1\cdot r^{3}).\end{array}$

Considering that the number of variables
$N_{v}(\mathcal{G}_{S_{q}}^{d}(r))=\mathcal{O}((2^{q})^{r})$, we estimate the
time complexity as a fuction of the number of variables $q\geq 1$ as:

$T(\mathcal{G}_{S_{q}}^{d}(r))=\mathcal{O}(N_{v}^{3\frac{q-1}{q}}).$ (22)

We sum up the analysis in Table 9.

Singularity type | Variables | Operations | Operations in $N_{v}$
---|---|---|---
Point | $\mathcal{O}(r)$ | $\mathcal{O}(r)$ | $\mathcal{O}(N_{v})$
Edge | $\mathcal{O}(2^{r})$ | $\mathcal{O}(2^{r})$ | $\mathcal{O}(N_{v})$
Face | $\mathcal{O}(4^{r})$ | $\mathcal{O}(8^{r})$ | $\mathcal{O}({N_{v}}^{\frac{3}{2}})$
Hyperface (3-D) | $\mathcal{O}(8^{r})$ | $\mathcal{O}(64^{r})$ | $\mathcal{O}({N_{v}}^{2})$
4-D | $\mathcal{O}(16^{r})$ | $\mathcal{O}(512^{r})$ | $\mathcal{O}({N_{v}}^{2.25})$
5-D | $\mathcal{O}(32^{r})$ | $\mathcal{O}(4096^{r})$ | $\mathcal{O}({N_{v}}^{2.4})$
$q$-D | $\mathcal{O}((2^{q})^{r})$ | $\mathcal{O}((2^{3(q-1)})^{r})$ | $\mathcal{O}({N_{v}}^{3\frac{(q-1)}{q}})$
Table 9: Number of variables and operations for different meshes and
singularity dimensions

###  B.3 Quasi-optimal h-adapted meshes towards a $q$ singularity

To generalize the analysis from the previous section, we can observe that a
broader class of meshes with singularities follow the same time complexity –
we define those meshes as quasi-optimal $q$-dimensional singularity meshes.
Examples of such meshes are presented in Figures 12 and 13.

Figure 12: Two-dimensional mesh refined towards non-regular “edge”
singularity. Figure 13: Three-dimensional mesh refined towards non-regular
“edge” singularity.

We say that a mesh is a quasi-optimal $q$-dimensional singularity mesh if it
has the following properties:

  1. 1.

The basis functions of the mesh can be assigned to tree nodes of a full
$K$-nary tree ($K\geq 2$) of height not larger than some
$R=\lceil{\log_{K}{N}}\rceil+R^{\prime}$, where $R^{\prime}$ is some constant.

  2. 2.

If a pair of basis functions have overlapping supports, they are assigned to
the same tree node, or one of them is assigned to an ancestor of the tree node
of the other one.

  3. 3.

Each tree node has not more than $QK^{h\frac{q-1}{q}}$ (if $q\geq 2$) or more
than $Qh$ (if $q\leq 1$) basis functions assigned, where $h$ is the height of
the subtree that given node is the root of and $Q$ is some arbitrary positive
constant. At the same time, the number of overlaps between basis functions
belonging to that tree with basis functions of ancestor tree nodes is limited
by the same number.

If such a tree is created, we can use it to define an ordering that would
follow the post-order traversal of that tree. If so, for a $q$-dimensional
singularity with $q\leq 1$, the complexity of removing all nodes belonging to
a tree node $S_{h}$, with height $h$, is bounded by the following equation:

$T(S_{h})\leq C_{r}(Qh,2Qh)=4Q^{3}h^{3}.$ (23)

The total complexity of the execution of the solver is in turn no more than:

$\begin{array}[]{l @{\;}
l}T(\mathcal{G}(r))&\displaystyle\leq\sum_{h=1}^{R}K^{R-h}\cdot
4(Q^{3}h^{3})=4Q^{3}K^{R}\sum_{i=1}^{R}K^{-h}h^{3}\vskip 3.0pt plus 1.0pt
minus 1.0pt\\\ &=\displaystyle
Q^{3}K^{R}\left(\frac{1^{3}}{K^{1}}+\frac{2^{3}}{K^{2}}+\dots+\frac{(R-1)^{3}}{K^{(}R-1)}+\frac{R^{3}}{K^{R}}\right)<Q^{3}K^{R}\cdot
26.\end{array}$ (24)

Thus,

$T(\mathcal{G}(r))=\mathcal{O}(K^{R})=\mathcal{O}(N).$ (25)

In the case of $q$-dimensional singularity with $q\geq 2$, the complexity of
removing all nodes belonging to the tree node $S_{h}$, with height $h$, is
bounded by the following equation:

$T(S_{h})\leq C_{r}(Q\cdot K^{h\frac{q-1}{q}},2Q\cdot
K^{h\frac{q-1}{q}})=4Q^{3}\cdot K^{3h\frac{q-1}{q}}$ (26)

Therefore, the total complexity of running the whole solver is no more than:

$\begin{array}[]{l @{\;}
l}T(\mathcal{G}(r))&\displaystyle=\sum_{h=1}^{R}K^{R-h}\cdot
4Q^{3}(K^{3h\frac{q-1}{q}})=4Q^{3}\sum_{h=1}^{R}K^{R-h+3h\frac{q-1}{q}}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ &=\displaystyle
4Q^{3}\sum_{h=1}^{R}K^{R}K^{h\frac{3(q-1)-q}{q}}=4Q^{3}\sum_{h=1}^{R}K^{R}(K^{\frac{2q-3}{q}})^{h}\vskip
3.0pt plus 1.0pt minus 1.0pt\\\ &=\displaystyle
4Q^{3}K^{R}\sum_{h=1}^{R}(K^{\frac{2q-3}{q}})^{h}<4Q^{3}K^{R}(K^{\frac{2q-3}{q}})^{R+1}.\end{array}$
(27)

Thus,

$\begin{array}[]{l @{\;}
l}T(\mathcal{G}(r))&=\displaystyle\mathcal{O}(K^{R}(K^{\frac{2q-3}{q}})^{R})=\mathcal{O}(K^{R+\frac{2q-3}{q}R})=\mathcal{O}(K^{(1+\frac{2q-3}{q})R})=\mathcal{O}(K^{3\frac{q-1}{q}R})\\\
&=\displaystyle\mathcal{O}(K^{3\frac{q-1}{q}\log_{K}{N}})=\mathcal{O}((K^{\log_{K}{N}})^{3\frac{q-1}{q}})=\mathcal{O}(N^{3\frac{q-1}{q}}).\end{array}$
(28)

Summarizing, both cases can be stated as:

$\displaystyle
T(\mathcal{G}(r))=\mathcal{O}\left(N^{\min\left\\{3\frac{q-1}{q},1\right\\}}\right).$
(29)

###  B.4 Quasi-optimality of $h$-adaptive meshes around singularities

To prove that every $h$-adaptive mesh adapted towards a singularity using the
Algorithm 1 is quasi-optimal, we propose the following tree generation
algorithm:

  1. 1.

Find a dividing plane that crosses the least amount of basis functions’
supports out of all planes perpendicular to the singularity that divide the
mesh so that no more than half of all basis functions lay solely on either one
of the sides of the plane. Create a tree node with all basis functions that
the chosen plane crosses the support of and remove them from the mesh.

  2. 2.

For each side of the plane, take the basis functions on that side and
recursively run the algorithm. The resulting trees become subtrees of the node
created in the previous procedure.

  3. 3.

Return the tree rooted in the node created in the first step.

Let us analyze if a tree generated that way proves the quasi-optimality of the
mesh:

  1. 1.

Every child of any node has at most half of the basis functions of its parent.
Because of that, the total height of the tree cannot be larger than
$\lceil{\log_{2}{N_{v}}}\rceil$ ($K=2$).

  2. 2.

There is no overlap between supports of basis functions of either side of the
dividing plane.

  3. 3.

Elements of refinement level $r$ cannot be produced farther than $2\sqrt{d}$
times the side of those elements. Thus the number of elements of a given
refinement level crossing the dividing plane is limited.

    1. (a)

In the case of $q=1$, all the dividing planes are parallel to each other.
Furthermore, the dividing plane of tree node $g$ nodes deep from the root has
a distance between the nearest planes of tree nodes higher in the gree of at
most $2^{-g}\cdot L$, where $L$ is the length of the singularity. Last means
that no elements larger than $2^{-g}\cdot L$ are eliminated at this node. In
other words, while each dividing plane crosses at most some $Q$ elements of
each refinement level ($Q$ is an arbitrary constant), the minimal refinement
level of elements removed at given tree node increases by $1$ every step down
the tree. Last means that the root node has at most $Q\cdot R$ elements, and
other nodes have at most $Q\cdot h$ elements.

    2. (b)

In case of $q\geq 2$, the limit of elements of refinement level $r$ crossed by
a dividing plane is $\mathcal{O}(2^{r(q-1)}\cdot L)$, where $L$ is the
length/area/volume of cross section between the singularity and the dividing
plane (assuming that the whole mesh has side of length $1$). Because the plane
of the smallest cross-section is chosen, $L$ decreases on average by
$2^{\frac{q-1}{q}}$. The total number of elements crossed are
$\mathcal{O}(2^{r(q-1)}\cdot((\frac{1}{2})^{\frac{q-1}{q}})^{R-h})$.
Considering that $R=rq$ (up to a constant), this is equivalent to
$\mathcal{O}(2^{r(q-1)}\cdot((\frac{1}{2})^{\frac{q-1}{q}})^{rq-h})=\mathcal{O}(2^{h\frac{q-1}{q}})$.

## References

  * [1] U. Langer and O. Steinbach (Eds.). Space-Time Methods: Applications to Partial Differential Equations (Vol. 25). Walter de Gruyter GmbH & Co KG. (2019)
  * [2] V. Thomée, Galerkin finite element methods for parabolic problems (Vol. 25). Springer Science & Business Media (2007).
  * [3] C. A. de Moura, C. S. Kubrusly, The Courant-Friedrichs-Lewy (CFL) Condition: 80 Years After Its Discovery, Birkhauser, (2013).
  * [4] J. Hartmanis and R. Stearns:  On the computational complexity of algorithms.  Transactions of the American Mathematical Society, (117) (1965) 285-306.
  * [5] I. S. Duff, J. K., Reid The multifrontal solution of indefinite sparse symmetric linear systems. ACM Transactions on Mathematical Software, 9 (1983) 302-325.
  * [6] I. S. Duff, J. K. Reid, The multifrontal solution of unsymmetric sets of linear systems., SIAM Journal on Scientific and Statistical Computing, 5 (1984) 633-641.
  * [7] J.W.H. Liu, The multifrontal method for sparse matrix solution: theory and practice, SIAM Review 34 (1992), 82-109.
  * [8] V. M. Calo, N. Collier, D. Pardo, M. Paszyński, Computational complexity and memory usage for multi-frontal direct solvers used in p finite element analysis, Procedia Computer Science 4 (2011) 1854-1861.
  * [9] M. Paszyński, V. M. Calo, D. Pardo, Direct solvers performance on h-adapted grids Computers & Mathematics with Applications 70 (3) (2015) 282-295.
  * [10] M. Yannakakis, Computing the minimum fill-in is NP-complete, SIAM Journal on Algebraic Discrete Methods, 2 (1981) 77-79.
  * [11] W. Hackbush, Hierarchical Matrices: Algorithms and Analysis, Springer (2015)
  * [12] G. Karypis, V. Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM Journal of Scientiffic Computing, 20, 1 (1998) 359-392.
  * [13] P. Heggernes, S.C. Eisenstat, G. Kumfert, A. Pothen, The Computational Complexity of the Minimum Degree Algorithm, ICASE Report No. 2001-42, (2001).
  * [14] J. Schulze, Towards a tighter coupling of bottom-up and top-down sparse matrix ordering methods, BIT, 41, 4 (2001) 800.
  * [15] P. R. Amestoy, T. A. Davis, I. S. Du, An Approximate Minimum Degree Ordering Algorithm, SIAM Journal of Matrix Analysis & Application, 17, 4 (1996) 886-905.
  * [16] G.W. Flake, R.E. Tarjan, K. Tsioutsiouliklis, Graph clustering and minimum cut trees, Internet Mathematics 1 (2003), 385-408.
  * [17] A. Paszyńska, Volume and neighbors algorithm for finding elimination trees for three dimensional h-adaptive grids, Computers & Mathematics with Applications, 68(10) (2014) 1467-1478.
  * [18] A. Paszyńska, M. Paszyński, K. Jopek, M. Woźniak, D. Goik, P. Gurgul, H. AbouEisha, M. Moshkov, V. M. Calo, A. Lenharth, D. Nguyen, K. Pingali, Quasi-optimal elimination trees for 2D grids with singularities, Scientific Programming, Article ID 303024 (2015), 1-18.
  * [19] H. AbouEisha, V. M. Calo, K. Jopek, M. Moshkov, A. Paszyńska, M. Paszyński, Bisections-Weighted-by-Element-Size-and-Order Algorithm to Optimize Direct Solver Performance on 3D hp-adaptive Grids, Lecture Notes in Computer Science, 10861 (2018) 760-772.
  * [20] Demkowicz, L., Kurtz, J., Pardo, D., Paszyński, M., Rachowicz, W., & Zdunek, A. (2007). Computing with hp-Adaptive Finite Elements, Vol. II. Frontiers. Three Dimensional Elliptic and Maxwell Problems with Applications. Chapman and Hall/Crc Applied Mathematics and Nonlinear Science.
  * [21] Y. Saad, Iterative Methods for Sparse Linear Systems, Society for Industrial and Applied Mathematics; 2nd edition (2003)
  * [22] M. Skotniczny, Computational Complexity of Hierarchically Adapted Meshes, Lecture Notes in Computer Science volume 12139 (2020) 226–239.
  * [23] D. Pardo, J. Álvarez-Aramberri, M. Paszyński, L. Dalcin, V. M. Calo, Impact of element-level static condensation on iterative solver performance, Computers & Mathematics with Applications 70 (10), (2015) 2331-2341.
  * [24] M. Woźniak, M. Paszyński, D. Pardo, L. Dalcin, V. M. Calo, Computational cost of isogeometric multi-frontal solvers on parallel distributed memory machines, Computer Methods in Applied Mechanics and Engineering 284 (2015) 971-987.
  * [25] M. Woźniak, K. Kuźnik, M Paszyński, V. M. Calo, D. Pardo, Computational cost estimates for parallel shared memory isogeometric multi-frontal solvers, Computers & Mathematics with Applications 67 (10) (2014) 1864-1883.