A massively parallel explicit solver for elasto-dynamic problems
exploiting octree meshes

unsw]Junqi Zhangcor1

[cor1]Corresponding author

unsw]Ankit Ankit

ude]Hauke Gravenkamp

unsw]Sascha Eisenträger

unsw]Chongmin Song

[unsw]School of Civil and Environmental Engineering, University of New
South Wales, Sydney NSW 2052, Australia

[ude]Department of Civil Engineering, University of Duisburg-Essen, 45141
Essen, Germany

Typical areas of application of explicit dynamics are impact, crash
test, and most importantly, wave propagation simulations. Due to the
numerically highly demanding nature of these problems, efficient automatic
mesh generators and transient solvers are required. To this end, a
parallel explicit solver exploiting the advantages of balanced octree
meshes is introduced. To avoid the hanging nodes problem encountered
in standard finite element analysis (FEA), the scaled boundary finite
element method (SBFEM) is deployed as a spatial discretization scheme.
Consequently, arbitrarily shaped star-convex polyhedral elements are
straightforwardly generated. Considering the scaling and transformation
of octree cells, the stiffness and mass matrices of a limited number
of unique cell patterns are pre-computed. A recently proposed mass
lumping technique is extended to 3D yielding
a well-conditioned diagonal mass matrix. This enables us to leverage
the advantages of explicit time integrator, i.e., it is possible to
efficiently compute the nodal displacements without the need for solving
a system of linear equations. We implement the proposed scheme together
with a central difference method (CDM) in a distributed computing
environment. The performance of our parallel explicit solver is evaluated
by means of several numerical benchmark examples, including complex
geometries and various practical applications. A significant speedup
is observed for these examples with up to one billion of degrees of
freedom and running on up to 16,384 computing cores.

Explicit dynamics; Parallel computing; Mass lumping; Octree mesh;
Scaled boundary finite element method

§ INTRODUCTION

Numerical methods for the solution of partial differential equations
(PDEs) have been widely used in engineering due to their feasibility
and reliability in handling problems with complex geometries and boundary
conditions. The finite element method (FEM) is one of the most popular
numerical methods, in which a computational domain is spatially discretized
into small subdomains of simple shapes, referred to as elements. As
a result, the governing PDEs are transformed into a system of linear
algebraic equations, which can be easily solved using modern computers.

Dynamic loading is ubiquitous in engineering practice and can be caused
by various physical processes, e.g., earthquakes [Dupros et al., 2010],
blasts [Brun et al., 2012], impacts [Bettinotti et al., 2017] and many
others [Zhang and Chauhan, 2019, Meister et al., 2020]. The FEM has been
used in structural dynamics to obtain the response history of structures,
which can be determined based on several methods, such as the modal
superposition technique, the Ritz-vector analysis or direct time integration
schemes [Cook et al., 1989]. Direct time integration schemes are often
based on finite difference approximations to achieve a temporal discretization.
Equilibrium equations are solved at each time step to compute the
unknown variables (displacements). Such integration schemes can be
classified further into two categories, namely, explicit and implicit
methods [Brun et al., 2012]. In an explicit method, the results of
the current time step are calculated based on one or more previous
time steps only, while in implicit methods the result depends not
only on previous time steps but also on the nodal velocity and acceleration
of the current step [Patzák et al., 2001]. In explicit methods, if
the mass and damping matrices are diagonal, the new nodal displacements
can be obtained efficiently without solving a system of linear equations [Talebi et al., 2012].
Furthermore, the internal force vector can be assembled independently,
i.e., element-by-element (EBE). Therefore, it is not necessary to
assemble the global stiffness matrix, which yields advantages with
respect to the required memory. Due to the possibility of employing
EBE techniques, we are in a position to exploit parallel processing
strategies, which can significantly improve the numerical efficiency
by making use of the capabilities of multi-processor computers. In
recent years, the parallel implementation of explicit dynamics has
also been investigated in other numerical approaches such as meshfree
methods [Li et al., 2014], the virtual element
method [Park et al., 2019] and isogeometric analysis [Anitescu et al., 2019].

Parallel processing is an attractive technique in high-performance
computing to speed up the algorithm by partitioning a serial job into
several smaller jobs and distributing them to individual processors
for execution [Rek and Němec, 2016]. A commonly used parallel processing
paradigm is single program, multiple data (SPMD)[Note that in Flynn's taxonomy (classification of computer architectures)
the same approach is referred to as single instruction, multiple data
(SIMD) which exploits data level parallelism.], i.e., different processors execute the same command on different
data sets simultaneously [Altman, 2014]. In explicit dynamics,
the element set can be partitioned and assigned to several processors
so that several local nodal force vectors can be assembled at the
same time, which will then be synchronized across all processors to
form a global nodal force vector [Ma et al., 2020]. To this end, the
workload should be distributed equally to the processors, and the
amount of data communication should be as small as possible, because
load-imbalance will hinder the data synchronization, which constitutes
the bottleneck in current parallel processing strategies and limits
the attainable speedup. It is worth mentioning that GPU-accelerated
parallelization also attracted a lot of research interest with the
development of modern GPU architectures such as CUDA [Tang et al., 2020].
In the current contribution, our focus is on parallelization based
on the use of multiple CPUs, and therefore, GPU-based schemes are
only mentioned for the sake of completeness but certainly constitute
a very promising direction for future applications.

Here, an explicit dynamic analysis methodology based on the use of
balanced octree meshes is presented. Due to the ease of creating such
meshes in the scaled boundary finite element method (SBFEM), our approach
utilizes this particular numerical method for the spatial discretization.
The SBFEM was first proposed by Song and Wolf, 1997 and initially developed
for dynamic analysis in unbounded domains and later extended to a
broad range of applications [Deeks and Wolf, 2002, Song, 2004, Bazyar and Song, 2008].
It has since been developed into a general-purpose numerical method
for the solution of PDEs [Ooi et al., 2020].
The basic idea is to divide the problem domain into subdomains satisfying
the so-called scaling requirement (Section <ref>).
In the SBFEM, a semi-analytical approach is adopted to construct an
approximate solution in a subdomain. To this end, only the boundary
of the subdomain is discretized in a finite element sense, and the
solution in the radial direction must be obtained analytically. This
feature of the SBFEM enables a straightforward derivation of polytopal
(polygonal [Chiong et al., 2014, Gravenkamp and Natarajan, 2018]/polyhedral [Talebi et al., 2016])
elements with an arbitrary number of nodes, edges and faces to be
used in the analysis. Due to the versatility in element shapes, the
difficulty in mesh generation can be greatly reduced. In three-dimensional
analyses, the surface discretization typically consists of quadrilateral
and triangular elements [Saputra et al., 2017]. However, it is well-known
that quadrilateral elements are more accurate, and therefore, in an
attempt to employ only quadrilateral elements in the surface discretization,
a novel transition element approach has been developed in Refs. [Saputra et al., 2020, Gravenkamp et al., 2020].
Based on this methodology, arbitrary element types can be coupled [Duczek et al., 2020].
In recent years, researchers have endeavored to apply this method
to address image-based analysis [Saputra et al., 2017, Saputra et al., 2017, Gravenkamp and Duczek, 2017],
acoustics [Liu et al., 2019, Liu et al., 2019], contact [Xing et al., 2018, Xing et al., 2019],
domain decomposition [Zhang and Song, 2019], wave propagation [Gravenkamp et al., 2012, Gravenkamp et al., 2012, Gravenkamp et al., 2017]
and many other problems [Zou et al., 2019, Zhang et al., 2020, Liu et al., 2020, Natarajan et al., 2020, Eisenträger et al., 2020].
The first application of SBFEM in the context of explicit dynamics
has been reported in Ref. [Qu et al., 2020].

In the proposed explicit method, a balanced octree mesh (a recursive
partition of the three-dimensional space into smaller octants where
adjacent elements can only differ by one octree level) is generated
to approximate the geometry. The octree mesh generation is highly
efficient. An octree based algorithm has been developed for automatic
mesh generation of digital images [Saputra et al., 2017] and stereolithography
(STL) models [Liu et al., 2017]. The application of octree meshes in
the standard finite element analysis is hindered by the existence
of hanging nodes. On the other hand, octree meshes are highly complementary
with the SBFEM [Song, 2018, Zhang et al., 2020]. An advantage of the
octree meshes benefiting parallel computing is that only a limited
number of octree patterns are present, and therefore, the element
stiffness and mass matrices can be pre-computed. The nodal force vector
is calculated element-wise without the need for assembling the global
stiffness matrix, while the global mass matrix is lumped into a well-conditioned
diagonal matrix based on a recently developed technique [Gravenkamp et al., 2020].
All elements of the same shape (i.e., same nodal pattern) are grouped,
and their nodal displacement vectors are assembled as a matrix. In
the next step, the nodal forces are calculated by a simple matrix
multiplication. A parallel version of the solver is developed, in
which the mesh is first partitioned into several parts of almost equal
size with minimum connection [Chan et al., 1994], and each
processor assembles the nodal force vector of one individual part.
Data exchange between individual processors is only needed at the
interfaces of the connected parts. The local nodal displacement vector
for the next time step is efficiently computed on each processor due
to the diagonal structure of the mass matrix before a global update
is conducted.

The remainder of this article is organized as follows: In Section <ref>,
a brief description of polyhedral elements constructed using the SBFEM
is provided. The central difference method, a widely used explicit
time integration scheme which forms the basis of our parallel solver,
and its implementation using octree mesh is presented in Section <ref>.
Our parallel processing strategy is discussed in Section <ref>.
In Section <ref>, the performance including the attainable
speedup is demonstrated by means of five numerical examples of increasing
complexity of the geometric models, number of degrees of freedom and
computing cores. General remarks and a summary are given in Section <ref>.

§ THE SCALED BOUNDARY FINITE ELEMENT METHOD

In this section, the construction of polyhedral elements using the
SBFEM is briefly discussed. For the sake of brevity, only the key
equations that are necessary for the implementation are presented.
Readers interested in a detailed derivation and additional explanations
are referred to the monograph by Song, 2018.

§.§ Polyhedral elements based on scaled boundary finite element method

The derivation of a polyhedron as a volume element in the SBFEM is
illustrated in Fig. <ref>. A so-called scaling center
(denoted as “$O$” in Fig. <ref>) is defined,
which should be directly visible from the entire boundary of the volume
element (condition of star-convexity). This requirement on the geometry
can always be satisfied by subdividing a concave polyhedron. One distinct
advantage of the SBFEM is that only the boundary of the volume element
needs to be discretized by means of surface elements. The volume of
the element is then obtained by scaling the discretized boundary towards
the scaling center. In this work, only four-node quadrilateral and
three-node triangular elements are used as surface elements.

A polyhedral element obtained by scaling the surface mesh toward the
scaling center $O$

The surface elements that are required to discretize the boundary
of the SBFEM volume element are conventional finite elements. Triangular
and quadrilateral elements in their local coordinate systems are shown
in Fig. <ref> and Fig. <ref>,
respectively. The shape functions of a surface element are assembled
in the matrix $\mathrm{\mathbf{N}}\!\left(\eta,\zeta\right)$, where
$\eta$ and $\zeta$ are the local coordinates of the element. As
illustrated in Fig. <ref>, a dimensionless radial
coordinate $\xi$ spanning $\left[0,1\right]$ from the scaling center
$O$ to the boundary is introduced for the purpose of scaling.

[A triangular element]

[A quadrilateral element]
[Radial coordinate $\xi$ emanating from the scaling center $O$ to
the boundary of the volume element]

Scaled boundary coordinates $\xi$, $\eta$ and $\zeta$

In the scaled boundary coordinates $\left(\xi,\eta,\zeta\right)$,
a point $\left(\hat{x},\hat{y},\hat{z}\right)$ inside the domain
can be defined as

\begin{align}
\hat{x}\!\left(\xi,\eta,\zeta\right) & =\xi x\!\left(\eta,\zeta\right)=\xi\mathrm{\mathbf{N}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{x}},\\
\hat{y}\!\left(\xi,\eta,\zeta\right) & =\xi y\!\left(\eta,\zeta\right)=\xi\mathrm{\mathbf{N}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{y}},\\
\hat{z}\!\left(\xi,\eta,\zeta\right) & =\xi z\!\left(\eta,\zeta\right)=\xi\mathrm{\mathbf{N}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{z}},
\end{align}

where $\mathbf{x},\mathbf{y},\mathbf{z}$ are
the nodal coordinate vectors of the surface element in Cartesian coordinates.
The corresponding Jacobian matrix on the boundary $\left(\xi=1\right)$

\begin{equation}
\mathrm{\mathbf{J}}\!\left(\eta,\zeta\right)=\left[\begin{array}{ccc}
x\!\left(\eta,\zeta\right) & y\!\left(\eta,\zeta\right) & z\!\left(\eta,\zeta\right)\\
x\!\left(\eta,\zeta\right)_{,\eta} & y\!\left(\eta,\zeta\right)_{,\eta} & z\!\left(\eta,\zeta\right)_{,\eta}\\
x\!\left(\eta,\zeta\right)_{,\zeta} & y\!\left(\eta,\zeta\right)_{,\zeta} & z\!\left(\eta,\zeta\right)_{,\zeta}
\end{array}\right].\label{eq:J_bd}
\end{equation}

The unknown displacement functions $\mathrm{\mathbf{u}}\!\left(\xi\right)$
are introduced on the radial lines connecting the scaling center and
the boundary nodes

\begin{equation}
\mathrm{\mathbf{u}}\!\left(\xi\right)=\underset{e}{\sum}\mathrm{\mathbf{u}}^{e}\!\left(\xi\right),\label{eq:u_ue}
\end{equation}

where $\underset{e}{\sum}$ indicates a standard finite element assembly
process. Within the volume element, the displacements at a point inside
a sector, formed by a surface element $e$, are interpolated as

\begin{equation}
\mathrm{\mathbf{u}}\!\left(\xi,\eta,\zeta\right)=\mathrm{\mathbf{N}}^{\mathrm{u}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{u}}^{e}\!\left(\xi\right),\label{eq:u_Nuue}
\end{equation}

where $\mathrm{\mathbf{N}}^{\mathrm{u}}$ is the shape functions written in matrix form.
In scaled boundary coordinates, the linear differential operator $\mathrm{\mathbf{L}}$
for the strain-displacement relationship
is expressed as [Song and Wolf, 1997]

\begin{equation}
\mathrm{\mathbf{L}}=\mathrm{\mathbf{b}}_{\mathrm{1}}\!\left(\eta,\zeta\right)\frac{\partial}{\partial\xi}+\frac{1}{\xi}\!\left(\mathrm{\mathbf{b}}_{\mathrm{2}}\!\left(\eta,\zeta\right)\frac{\partial}{\partial\eta}+\mathrm{\mathbf{b}}_{\mathrm{3}}\!\left(\eta,\zeta\right)\frac{\partial}{\partial\zeta}\right),\label{eq:L}
\end{equation}

where $\mathrm{\mathbf{b}}_{1}\!\left(\eta,\zeta\right)$, $\mathrm{\mathbf{b}}_{2}\!\left(\eta,\zeta\right)$
and $\mathrm{\mathbf{b}}_{3}\!\left(\eta,\zeta\right)$ are defined
in [Song, 2018].

Using Eqs. (<ref>)
and (<ref>), the strains are expressed

\begin{equation}
\mathrm{\boldsymbol{\varepsilon}}\!\left(\xi,\eta,\zeta\right)=\mathrm{\mathbf{B}}_{\mathrm{1}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{u}}^{e}\!\left(\xi\right),_{\xi}+\frac{1}{\xi}\mathrm{\mathbf{B}}_{2}\!\left(\eta,\zeta\right)\mathrm{\mathbf{u}}^{e}\!\left(\xi\right),\label{eq:eps_B1B2ue}
\end{equation}

\begin{align}
\mathrm{\mathbf{B}}_{\mathrm{1}}\!\left(\eta,\zeta\right)= & \mathrm{\mathbf{b}}_{\mathrm{1}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{N}}^{\mathrm{u}}\!\left(\eta,\zeta\right),\label{eq:B1}\\
\mathrm{\mathbf{B}}_{\mathrm{2}}\!\left(\eta,\zeta\right)= & \mathrm{\mathbf{b}}_{\mathrm{2}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{N}}^{\mathrm{u}}\!\left(\eta,\zeta\right),_{\eta}+\mathrm{\mathbf{b}}_{\mathrm{3}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{N}}^{\mathrm{u}}\!\left(\eta,\zeta\right),_{\zeta}.\label{eq:B2}
\end{align}

Assuming linear elasticity, the stresses are equal

\begin{equation}
\mathrm{\boldsymbol{\sigma}}\!\left(\xi,\eta,\zeta\right)=\mathrm{\mathbf{D}}\!\left(\mathrm{\mathbf{B}}_{\mathrm{1}}\!\left(\eta,\zeta\right)\mathrm{\mathbf{u}}^{e}\!\left(\xi\right),_{\xi}+\frac{1}{\xi}\mathrm{\mathbf{B}}_{2}\!\left(\eta,\zeta\right)\mathrm{\mathbf{u}}^{e}\!\left(\xi\right)\right),\label{eq:sig_B1B2ue}
\end{equation}

where $\mathrm{\mathbf{D}}$ is the elasticity matrix.

§.§ Solution of the scaled boundary finite element equation for elasto-statics

The scaled boundary finite element equation in displacements can be
derived based on Galerkin’s method or the virtual work
principle, which were reported in Refs. [Song and Wolf, 1997] and [Deeks and Wolf, 2002],
respectively. For each volume element, the equation can be written

\begin{equation}
\mathrm{\mathbf{E}_{0}}\xi^{2}\mathrm{\mathbf{u}}\!\left(\xi\right),_{\xi\xi}+\left(2\mathbf{E}_{0}-\mathbf{E}_{1}+\mathbf{E}_{1}^{\mathrm{T}}\right)\xi\mathrm{\mathbf{u}}\!\left(\xi\right),_{\xi}+\left(\mathbf{E}_{1}^{\mathrm{T}}-\mathrm{\mathbf{E}_{2}}\right)\mathrm{\mathbf{u}}\!\left(\xi\right)=\boldsymbol{0}.\label{eq:SBFE}
\end{equation}

The coefficient matrices $\mathbf{E}_{0}^{e}$, $\mathbf{E}_{1}^{e}$
and $\mathbf{E}_{2}^{e}$ of a surface element $e$ are obtained by

\begin{align}
\mathbf{E}_{0}^{e}= & \int_{e}\mathrm{\mathbf{B}_{1}^{\mathrm{T}}\mathrm{\mathbf{D}}\mathbf{B}_{1}}\mathrm{|\mathrm{\mathbf{J}}|d\eta d\zeta},\\
\mathbf{E}_{1}^{e}= & \int_{e}\mathrm{\mathbf{B}_{2}^{\mathrm{T}}\mathrm{\mathbf{D}}\mathbf{B}_{1}}\mathrm{|\mathrm{\mathbf{J}}|d\eta d\zeta},\\
\mathbf{E}_{2}^{e}= & \int_{e}\mathrm{\mathbf{B}_{2}^{\mathrm{T}}\mathrm{\mathbf{D}}\mathbf{B}_{2}}\mathrm{|\mathrm{\mathbf{J}}|d\eta d\zeta}.
\end{align}

They are assembled to form the coefficient matrices
$\mathbf{E}_{0}$, $\mathbf{E}_{1}$ and $\mathbf{E}_{2}$ of a volume
element in Eq. (<ref>) according to the element connectivity
data. This procedure is identical to that known from the FEM. According
to Song, 2004, the internal nodal forces $\mathrm{\mathbf{q}}\!\left(\xi\right)$
on a surface with a constant $\xi$ can be formulated as

\begin{equation}
\mathrm{\mathbf{q}}\!\left(\xi\right)=\xi\!\left(\mathbf{E}_{0}\xi\mathrm{\mathbf{u}}\!\left(\xi\right),_{\xi}+\mathbf{E}_{1}^{\mathrm{T}}\mathrm{\mathbf{u}}\!\left(\xi\right)\right).\label{eq:q_E0E1}
\end{equation}

The order of the differential equations in Eq. (<ref>) can
be reduced by one.
The solution obtained
from eigenvalue decomposition is given in
the following form

𝐮(ξ)     =Φ_u1ξ^-λ-0.5𝐈𝐜_1,

𝐪(ξ)     =Φ_q1ξ^-λ+0.5𝐈𝐜_1,

where $-\boldsymbol{\lambda}$ is the
eigenvalues, $\mathbf{\mathbf{\mathrm{\Phi}}_{\mathrm{u1}}}$ and
$\mathbf{\mathbf{\mathrm{\Phi}}_{\mathrm{q1}}}$ are the sub-matrices
of the eigenvector matrix, while $\mathbf{c}_{1}$ and $\mathbf{c}_{2}$
denote the integration constants. On the boundary $\left(\xi=1\right)$, the stiffness
matrix of the volume element ($\mathrm{\mathbf{q}}\!\left(\xi=1\right)=\mathbf{K}\mathrm{\mathbf{u}}\!\left(\xi=1\right)$)
is obtained by eliminating the integration constants $\mathbf{c}_{1}$

\begin{equation}
\mathbf{K}=\mathbf{\mathbf{\mathrm{\Phi}}_{\mathrm{q1}}}\mathbf{\mathbf{\mathrm{\Phi}}}_{\mathbf{\mathrm{u1}}}^{-1}.\label{eq:K}
\end{equation}

In order to obtain the element mass matrix, a coefficient matrix $\mathbf{M}_{0}$
is introduced

\begin{align}
\mathbf{M}_{0} & =\underset{e}{\sum}\mathbf{M}_{0}^{e},
\end{align}

where the element coefficient matrix $\mathbf{M}_{0}^{e}$ is defined

\begin{align}
\mathbf{M}_{0}^{e} & =\int_{e}\mathrm{\mathbf{N}}^{\mathrm{u}}\!\left(\eta,\zeta\right)^{\mathrm{T}}\rho\mathrm{\mathbf{N}}^{\mathrm{u}}\!\left(\eta,\zeta\right)|\mathrm{\mathbf{J}}|\mathrm{d}\eta\mathrm{d}\zeta.
\end{align}

An abbreviation $\mathbf{m}_{0}$ is introduced as

\begin{align}
\mathbf{m}_{0} & =\mathbf{\mathrm{\Phi}}_{\mathrm{u1}}^{\mathrm{T}}\mathbf{M}_{0}\mathbf{\mathrm{\Phi}}_{\mathrm{u1}}.
\end{align}

To integrate analytically along the radial direction, a matrix $\mathbf{m}$
is constructed, the components of which can be computed analytically
based on the matrix $\mathbf{m}_{0}$

\begin{align}
\mathrm{m}_{ij} & =\int_{0}^{1}\xi^{\lambda_{i}}\mathrm{m}_{0ij}\xi^{\lambda_{j}}\xi^{2}\mathrm{d}\xi\nonumber \\
& =\frac{\mathrm{m}_{0ij}}{\lambda_{i}+\lambda_{j}+3},
\end{align}

and the consistent mass matrix $\mathbf{M}$ is written as

\begin{align}
\mathbf{M} & =\mathbf{\mathrm{\Phi}}_{\mathrm{u1}}^{-\mathrm{T}}\mathbf{m}\mathbf{\mathrm{\Phi}}_{\mathrm{u1}}^{-1}.
\end{align}

§ EXPLICIT DYNAMICS

In this section, the central difference equations for explicit dynamics
are presented, including the discussion of an SBFEM-specific mass
lumping technique and the application of octree cell patterns. The
explicit solver is first presented in a serial version, which lays
the foundation for the parallel explicit solver introduced in Section <ref>.

§.§ Explicit time integration

In our implementation of the central difference method (CDM), a fixed
time step size is assumed for the sake of simplicity of the algorithm.
However, it should be kept in mind that adaptive time-stepping techniques
would also be possible. The displacement vector of all degrees of
freedom (DOFs) at the $n$-th time step is referred to as $\mathbf{U}_{n}$.
The velocity and acceleration vectors are the derivatives of the displacement
vector with respect to time and are denoted by $\dot{\mathbf{U}}_{n}$
and $\ddot{\mathbf{U}}_{n}$, respectively. The equilibrium equation
at the $n$-th time step using Newton's second law can be written

\begin{equation}
\mathbf{M}\ddot{\mathbf{U}}_{n}+\mathbf{C}\dot{\mathbf{U}}_{n}+\mathbf{R}_{n}^{\mathrm{int}}=\mathbf{R}_{n}^{\mathrm{ext}},\label{eq:Equ}
\end{equation}

where $\mathbf{M}$ is the mass matrix, $\mathbf{C}$ denotes the
damping matrix, $\mathbf{R}_{n}^{\mathrm{ext}}$ represents the external
load vector, and the internal force vector for linear analyses is
given as

\begin{equation}
\mathbf{R}_{n}^{\mathrm{int}}=\mathbf{K}\mathbf{U}_{n}.\label{eq:Rint}
\end{equation}

Using the CDM, the derivatives of $\mathbf{U}_{n}$ can be approximated

\begin{align}
\dot{\mathbf{U}}_{n} & =\frac{1}{2\Delta t}\!\left(\mathbf{U}_{n+1}-\mathbf{U}_{n-1}\right),\\
\ddot{\mathbf{U}}_{n} & =\frac{1}{\Delta t^{2}}\!\left(\mathbf{U}_{n+1}-2\mathbf{U}_{n}+\mathbf{U}_{n-1}\right),
\end{align}

where $\Delta t$ denotes the time step. In the
next step, we substitute Eqs. (<ref>) into Eq. (<ref>),

\begin{equation}
\frac{1}{\Delta t^{2}}\mathbf{M}\!\left(\mathbf{U}_{n+1}-2\mathbf{U}_{n}+\mathbf{U}_{n-1}\right)+\frac{1}{2\Delta t}\mathbf{C}\!\left(\mathbf{U}_{n+1}-\mathbf{U}_{n-1}\right)+\mathbf{R}_{n}^{\mathrm{int}}=\mathbf{R}_{n}^{\mathrm{ext}},\label{eq:EquDn}
\end{equation}

which can be solved for the displacement at the next time step, i.e.,
$\mathbf{U}_{n+1}$, expressed as

\begin{align}
\left(\frac{1}{\Delta t^{2}}\mathbf{M}+\frac{1}{2\Delta t}\mathbf{C}\right)\mathbf{U}_{n+1} & =\mathbf{R}_{n}^{\mathrm{ext}}-\mathbf{R}_{n}^{\mathrm{int}}+\frac{2}{\Delta t^{2}}\mathbf{M}\mathbf{U}_{n}-\!\left(\frac{1}{\Delta t^{2}}\mathbf{M}-\frac{1}{2\Delta t}\mathbf{C}\right)\mathbf{U}_{n-1}.\label{eq:Dn1}
\end{align}

If mass-proportional damping is considered,

\begin{equation}
\mathbf{C}=\alpha\mathbf{M},\label{eq:C}
\end{equation}

where $\alpha$ is a constant, Eq. (<ref>) can be written

\begin{align}
\mathbf{U}_{n+1} & =\frac{\Delta t^{2}\mathbf{M}^{-1}\!\left(\mathbf{R}_{n}^{\mathrm{ext}}-\mathbf{R}_{n}^{\mathrm{int}}\right)+2\mathbf{U}_{n}-\!\left(1-\frac{\alpha}{2}\Delta t\right)\mathbf{U}_{n-1}}{1+\frac{\alpha}{2}\Delta t}.\label{eq:DnM}
\end{align}

We observe that the displacement at the next step is only related
to the displacement of the current and previous time steps. To start
the procedure, a fictitious point $\mathbf{U}_{-1}$ is added, which
is calculated as

\begin{equation}
\mathbf{U}_{-1}=\mathbf{U}_{0}-\Delta t\dot{\mathbf{U}}_{0}+\frac{\Delta t^{2}}{2}\ddot{\mathbf{U}}_{0}.\label{eq:D0}
\end{equation}

One important drawback of explicit methods is that they are generally
only conditionally stable, which means there exists a critical time
step $\Delta t_{\mathrm{cr}}$ which should not be exceeded, otherwise
the results will diverge. Thus, they are often only used in high-frequency
problems where a small time step size is already demanded by the physics
of the problem. The exact value of $\Delta t_{\mathrm{cr}}$ involves
the calculation of the maximum natural frequency of the whole structure,
which is time and memory consuming. In our implementation, the stability
limit is estimated by the maximum frequency on element level, i.e.,
the critical time step (in the CDM) is

\begin{equation}
\Delta t_{\mathrm{cr}}=\frac{2}{\omega_{\max}},\label{eq:tcr}
\end{equation}

where $\omega_{\max}$ is the maximum frequency of all elements [Cook et al., 1989, de Béjar and Danielson, 2016].

§.§ Implementation using octree mesh

In this section, we introduce the mesh generation technique based
on octree algorithm, which is highly complementary to the SBFEM. Unique
octree cell patterns are discussed, which enables an efficient calculation
of nodal force vector without the assembly of global stiffness matrix.

§.§.§ Balanced octree mesh for SBFEM

The octree algorithm is a highly-efficient hierarchical-tree based
technique for mesh generation. It discretizes a problem domain into
cells by recursively bisecting the cell edges until some specified
stopping criteria are met [Saputra et al., 2017]. An example of octree
mesh in shown in Fig. <ref>, in which each octree
cell is modeled as a scaled boundary finite element. An octree mesh
is usually balanced by limiting the length ratio between two adjacent
cells to $2$ ($2:1$ rule). All the possible node locations of a
balanced octree cell is indicated in Fig. <ref>.
There are $2^{12}=4096$ patterns in total <cit.>.

[A balanced octree mesh]

[A balanced octree cell indicating all possible corner ($\bullet$),
mid-edge ($\circ$) and center ($+$) nodes]

Example of balanced octree mesh and octree cell patterns [Song, 2018]

When considering rotations, the faces of the cubic octree cells have
only six unique node arrangements as shown in Fig. <ref>,
depending on the number of mid-edge nodes $n_{\mathrm{h}}$ and their
locations. The faces are treated as polygons in the SBFEM and discretized
using rectangular and isosceles triangular elements, which are considered
as high-quality elements.

[$n_{\mathrm{h}}=0$][$n_{\mathrm{h}}=1$][$n_{\mathrm{h}}=2$ adjacent edges]

[$n_{\mathrm{h}}=2$ opposite edges][$n_{\mathrm{h}}=3$][$n_{\mathrm{h}}=4$]

Discretization of the polygonal faces on the boundary of the cubic
octree volume elements. The number of 'hanging nodes' on the edges
for each discretization pattern is indicated by $n_{\mathrm{h}}$.

§.§.§ Unique octree cell patterns

The 4096 patterns can be further transformed into 144 unique cases
using rotation and mirroring operators. The 90 degree rotation matrices
about the $x$, $y$ and $z$ axes are

\begin{equation}
\mathbf{T}_{\mathrm{x}}=\left[\begin{array}{ccc}
1 & 0 & 0\\
0 & 0 & 1\\
0 & -1 & 0
\end{array}\right],\enskip\mathbf{T}_{\mathrm{y}}=\left[\begin{array}{ccc}
0 & 0 & 1\\
0 & 1 & 0\\
-1 & 0 & 0
\end{array}\right],\enskip\mathbf{T}_{\mathrm{z}}=\left[\begin{array}{ccc}
0 & 1 & 0\\
-1 & 0 & 0\\
0 & 0 & 1
\end{array}\right],\label{eq:Rotate}
\end{equation}

and the mirror matrices are

\begin{equation}
\mathbf{R}_{\mathrm{x}}=\left[\begin{array}{ccc}
-1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{array}\right],\enskip\mathbf{R}_{\mathrm{y}}=\left[\begin{array}{ccc}
1 & 0 & 0\\
0 & -1 & 0\\
0 & 0 & 1
\end{array}\right],\enskip\mathbf{R}_{\mathrm{z}}=\left[\begin{array}{ccc}
1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & -1
\end{array}\right].\label{eq:Reflect}
\end{equation}

These transformation matrices can be arranged and combined to form
48 unique transformations of an octree cell in 3D. Alternatively,
the 48 transformations can be constructed directly by filling $-1$
or $1$ to each row and column of a $3\times3$ empty matrix, i.e.,
there is exactly one nonzero value ($1$ or $-1$) in each row and
each column. The number of transformations is equal to $3\times2\times1\times2^{3}=48$.
The 48  transformations involve only permutations of nodal numbers in an element, which can be easily implemented using an indexing vector.
Each of the 4096 patterns can be transformed into one of the 144 unique
patterns using one of the 48 transformations. All the unique octree
patterns are presented in <ref>.

The 144 patterns are determined solely by the presence or absence of mid-edge nodes, which is independent of the polyhedral element formulation selected. In our implementation, the octree meshes with triangulated boundaries are handled by the SBFEM, however it is important to note that the same approach can also be applied to other polyhedral elements with well-conditioned lumped mass matrices. The use of element patterns significantly improves computational efficiency, in terms of time and memory requirements, especially on a modern high-performance computing system.

Two sample transformations are illustrated in Fig. <ref>.
The octree cell in Fig. <ref> can be transformed
into the corresponding master cell $\mathrm{c}_{1}$ in Fig. <ref>
by rotating 90 degrees about the $x$ axis and then rotating 90 degrees
about the $z$ axis (following the right hand rule as positive rotation
direction). The rotation matrix can be constructed as $\mathbf{Q}_{1}=\mathbf{T}_{\mathrm{x}}\mathbf{T}_{\mathrm{z}}$.
Similarly, the octree cell in Fig. <ref> can be transformed
into $\mathrm{c}_{2}$ in Fig. <ref> by mirroring
about the $y$-$z$ plane using $\mathbf{Q}_{2}=\mathbf{R}_{\mathrm{x}}$.

[Octree master cell: $\mathrm{c}_{1}$] [Octree cell with a pattern identical to $\mathrm{c}_{1}$][Transformation matrix $\mathbf{Q}_{1}=\mathbf{T}_{\mathrm{x}}\mathbf{T}_{\mathrm{z}}$]
0 & 1 & 0\\
0 & 0 & 1\\
1 & 0 & 0
\end{array}\right]$

[Octree master cell: $\mathrm{c}_{2}$] [Octree cell with a pattern identical to $\mathrm{c}_{2}$][Transformation matrix $\mathbf{Q}_{2}=\mathbf{R}_{\mathrm{x}}$]
-1 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 1
\end{array}\right]$

Transformation of octree cells

When calculating the nodal force vector of an octree cell, the nodal
displacement vector is transformed using the same transformation matrix
$\mathbf{Q}$, which is then multiplied by the stiffness matrix of
the corresponding master cell. The obtained nodal force vector must
be transformed back using the inverse of transformation matrix, $\mathbf{Q}^{-1}$.

§.§.§ Implementation exploiting octree cell patterns

As an alternative to Eq. (<ref>), $\mathbf{R}_{n}^{\mathrm{int}}$
can be calculated EBE so that the assembly of the global stiffness
matrix $\mathbf{K}$ is not required.

\begin{equation}
\mathbf{R}_{n}^{\mathrm{int}}=\stackrel[e=1]{n_{\mathrm{e}}}{\sum}\mathbf{K}_{e}\mathbf{U}_{e,n},\label{eq:Rinte}
\end{equation}

where $\sum$ represents the assembly process according to the element
connectivity, $n_{\mathrm{e}}$ is the total number of elements, $\mathbf{K}_{e}$
is the element stiffness matrix of volume element $e$, and $\mathbf{U}_{e,n}$
is the nodal displacement vector of the element at the $n$-th time

In an octree mesh, elements of the same shape (see <ref>)
can be grouped together, and therefore, their nodal force vectors
can be calculated in one step using a single matrix-multiplication.
Eq. (<ref>) can be rewritten as

\begin{equation}
\mathbf{R}_{n}^{\mathrm{int}}=\stackrel[p=1]{n_{\mathrm{p}}}{\sum}\mathbf{K}_{p}\!\left[\begin{array}{cccc}
\mathbf{U}_{p1} & \mathbf{U}_{p2} & \cdots & \mathbf{U}_{pk}\end{array}\right]_{n}\mathrm{diag}\!\left(\mathbf{S}\right),\label{eq:Rintep}
\end{equation}

where $n_{\mathrm{p}}$ is the number of unique patterns in the octree
mesh, $\mathbf{K}_{p}$ is the stiffness matrix of a master element
of pattern $p$ with unit size, $\mathbf{U}_{pi}$ is the nodal displacement
vector of the $i$-th element of pattern $p$, and $\mathrm{diag}\!\left(\mathbf{S}\right)$
is a diagonal matrix containing the edge lengths of each element of
pattern $p$. Each column obtained from the right-hand side will be
assembled into $\mathbf{R}_{n}^{\mathrm{int}}$ using the element
connectivity data. As mentioned before, the spatial discretization
is based on a balanced octree mesh. This holds two important advantages:
(i) all surface elements are well-shaped and (ii) for a cubic volume
element only a limited number of surface patterns exist. These are
the prerequisites to exploit a pre-computation technique. To this
end, the stiffness and mass matrices of all master elements corresponding
to the different patterns $p$ can be computed beforehand. For the
pre-computation strategy, a unit size of the master element is assumed.
Due to the limited surface element shapes, the stiffness and mass
matrices for the master element only need to be scaled by suitable
factors including information on the element size and material properties.

The stiffness matrix of an element is proportional to its Young's
modulus and edge length. On the other hand, the mass matrix of an
element is proportional to the mass density and the cube of edge length.
Here, we define the ratios between Young's modulus, mass density and
edge length as $\eta_{E}=\nicefrac{E_{2}}{E_{1}}$, $\eta_{\rho}=\nicefrac{\rho_{2}}{\rho_{1}}$
and $\eta_{L}=\nicefrac{L_{2}}{L_{1}}$ (see Fig. <ref>).
Therefore, the stiffness and mass matrices, $\mathbf{K}_{2}$ and
$\mathbf{M}_{2}$, of an element can be derived easily from a master
element of same pattern and Poisson's ratio as

\begin{equation}
\mathbf{K}_{2}=\eta_{E}\eta_{L}\mathbf{K}_{1},
\end{equation}

\begin{equation}
\mathbf{M}_{2}=\eta_{\rho}\eta_{L}^{3}\mathbf{M}_{1},
\end{equation}

where $\mathbf{K}_{1}$ and $\mathbf{M}_{1}$ are
stiffness and mass matrices of the master element, respectively. In
cases where materials with different Poisson's ratios are present,
the master element matrices have be computed for each material separately.
The critical time step is estimated based on the maximum eigenfrequency
of the element in Eq. (<ref>), therefore, considering Eq. (<ref>),
$\Delta t_{\mathrm{cr}}{}_{2}$ of an element can be calculated from
the known time step of the master element $\Delta t_{\mathrm{cr}}{}_{1}$,

\begin{equation}
\Delta t_{\mathrm{cr}}{}_{2}=\sqrt{\frac{\eta_{\rho}}{\eta_{E}}}\eta_{L}\Delta t_{\mathrm{cr}}{}_{1}.\label{eq:tcr_scale}
\end{equation}

[Master element][Scaled element]

Scaling of an octree element

§.§ Mass lumping in the SBFEM

It can be inferred from Eq. (<ref>) that if the mass matrix
$\mathbf{M}$ is given in diagonal form, it is possible to advance
in time by simple matrix-vector products without the need for solving
a system of linear equations. The reason for this increased efficiency
lies in the trivial inversion of a diagonal matrix. In the context
of the finite and spectral element methods (SEM), different mass lumping
schemes have been reported and thoroughly assessed in Refs. [Duczek and Gravenkamp, 2019, Duczek and Gravenkamp, 2019].

Different from the FEM and SEM, the displacement interpolations inside
the scaled boundary finite elements are obtained semi-analytically,
and the $x$, $y$ and $z$ displacement components are coupled. In
our proposed mass lumping scheme, only those components of the mass
matrix that correspond to DOFs in the same coordinate direction are
summed up. The basic approach is illustrated in Fig. <ref>
for the DOFs related to $x$ direction. The procedure is repeated
in the same fashion for the components corresponding to $y$ and $z$
directions. This method results in a block-diagonal matrix as shown
in Fig. <ref>, from which a diagonal mass matrix
$\mathbf{M}_{\mathrm{diag}}$ is then extracted (Fig. <ref>).
The off-diagonal terms, which sum up to zero, are ignored. The total
mass of the system is conserved in the formulation of the SBFEM. A
detailed explanation and verification of this mass lumping approach
has been reported by Gravenkamp et al., 2020.

[Fully populated mass matrix $\mathbf{M}$ ]    [Block-diagonal mass matrix]

[Diagonal mass matrix $\mathbf{M}_{\mathrm{diag}}$]

Lumping of the mass matrix

Due to the mass matrix $\mathbf{M}$ being available in diagonal format
$\mathbf{M}_{\mathrm{diag}}$, we can re-write the expression for
the displacement vector in the next time step as

\begin{equation}
\mathbf{U}_{n+1}=\frac{\Delta t^{2}\mathbf{M}_{\mathrm{diag}}^{-1}\!\left(\mathbf{R}_{n}^{\mathrm{ext}}-\mathbf{R}_{n}^{\mathrm{int}}\right)+2\mathbf{U}_{n}-\!\left(1-\frac{\alpha}{2}\Delta t\right)\mathbf{U}_{n-1}}{1+\frac{\alpha}{2}\Delta t}.\label{eq:DnMdiag}
\end{equation}

Lumping the mass matrix avoids solving a system of linear equations
and as a result, each component in $\mathbf{U}_{n+1}$ can be calculated
independently of the other entries. Lumped mass matrices of several
octree patterns mentioned in Section <ref> are presented
below as ratios to the total mass of the octree cells. Poisson's ratios
are 0 and 0.3 in Fig. <ref> and Fig. <ref>,

[Element with 8 nodes][Element with 11 nodes]

[Element with 13 nodes][Element with 26 nodes]

Lumped mass matrices with Poisson's ratio $\nu=0$

[Element with 8 nodes][Element with 11 nodes]

[Element with 13 nodes][Element with 26 nodes]

Lumped mass matrices with Poisson's ratio $\nu=0.3$

§.§ Explicit solver in serial

Based on the previous sections, an explicit solver can be developed
making use of lumped mass matrices and the unique octree patterns.
The data structure of the solver is presented in Table <ref>.
Element stiffness and mass matrices of the 144 unique octree patterns
can be pre-computed and stored as $\mathbf{K}_{e}$ and $\mathbf{M}_{e}$.
It is noted that if multiple materials with different Poisson's ratios
are present in the model, each material will be assigned its own $\mathbf{K}_{e}$
and $\mathbf{M}_{e}$ matrices. For each specific problem, a matrix
$\mathbf{V}$ is constructed by concatenating the DOF vector of each
element, which is already arranged according to the transformation
matrix to the master cells. Furthermore, $n_{\mathrm{e}}$ is the
number of elements, and $n_{\mathrm{dof}}$ denotes the total number
of DOFs. The vector $\mathrm{\boldsymbol{\varGamma}_{D}}$ stores
the DOFs which are fixed (Dirichlet boundary conditions), and $\mathrm{\boldsymbol{\varGamma}_{N}}$
is a sparse vector containing the location and magnitude of the external
nodal forces (Neumann boundary conditions). $\mathrm{\mathbf{T}}$
is an integer vector indicating the element type, while $\mathrm{\mathbf{S}}$
is a vector containing the edge length of each element. The variables
mentioned above are constant for each problem during the calculation.
In each time step, the displacement vectors corresponding to the previous
step, current step and the next step, $\mathbf{U}_{\mathrm{p}}$,
$\mathbf{U}_{\mathrm{c}}$ and $\mathbf{U}_{\mathrm{n}}$, are updated
accordingly. $\mathbf{R}^{\mathrm{ele}}$ is a matrix containing the
nodal force vector of each element, which are grouped by element type.
$\mathbf{R}^{\mathrm{int}}$ is the internal force vector which is
assembled from $\mathbf{R}^{\mathrm{ele}}$, while $\mathbf{R}^{\mathrm{ext}}$
is the external force vector calculated at the time $t$. The pseudo-code
of the explicit solver is presented in Algorithm <ref>.

Data structure of explicit solver

Variable     Data type     Variable type     Size     Description20.1Constant     $\mathbf{K}_{e}$     float     hyper-matrixHyper-matrix: This term denotes that each component of the matrix
is a matrix itself. In MATLAB this can be for example realized as
a cell-array.     144     Element stiffness matrices of the 144 unique octree cell patterns2-6 3-6 4-6 5-6 6-6
$\mathbf{M}_{e}$     float     hyper-matrix     144     Diagonal mass matrices of the 144 unique octree cell patterns60.1Pre-compute     $\mathbf{V}$     integer     hyper-vector     $n_{\mathrm{e}}$     DOF vector of each element2-6 3-6 4-6 5-6 6-6
$\boldsymbol{\varGamma}_{\mathrm{D}}$     integer     vector     < $n_{\mathrm{dof}}$     Dirichlet boundary condition (the DOFs which are fixed)2-6 3-6 4-6 5-6 6-6
$\boldsymbol{\varGamma}_{\mathrm{N}}$     float     vector     $n_{\mathrm{dof}}$     Neumann boundary condition (the location and magnitude of external
forces)2-6 3-6 4-6 5-6 6-6
$\mathbf{T}$     integer     vector     $n_{\mathrm{e}}$     Pattern of each element, from 1 to 1442-6 3-6 4-6 5-6 6-6
$\mathbf{S}$     float     vector     $n_{\mathrm{e}}$     Edge length of each element2-6 3-6 4-6 5-6 6-6
$\mathbf{M}_{\mathrm{diag}}$     float     matrix     $n_{\mathrm{dof}}$     Global diagonal mass matrix30.1Time stepping     $\mathbf{U}_{\mathrm{p}}$, $\mathbf{U}_{\mathrm{c}}$, $\mathbf{U}_{\mathrm{n}}$     float     vector     $n_{\mathrm{dof}}$     Displacement vector of previous step, current step and next step2-6 3-6 4-6 5-6 6-6
$\mathbf{R}^{\mathrm{ele}}$     float     hyper-vector     $n_{\mathrm{e}}$     Internal nodal force vector of each element2-6 3-6 4-6 5-6 6-6
$\mathbf{R}^{\mathrm{ext}}$, $\mathbf{R}^{\mathrm{int}}$     float     vector     $n_{\mathrm{dof}}$     External and internal force vectors

[1]

Input mesh

Load pre-computed master element stiffness and mass matrices
$\mathbf{K}_{e}$ and $\mathbf{M}_{e}$

Identify element patterns $\mathrm{\mathbf{T}}$ and edge
lengths $\mathbf{S}$, and link to corresponding master element

Assemble global $\mathbf{V}$ matrix and group based on element

Assemble global diagonal mass matrix $\mathbf{M}_{\mathrm{diag}}$

Apply boundary conditions

Initialize nodal displacement vectors $\mathbf{U}_{\mathrm{p}}$,
$\mathbf{U}_{\mathrm{c}}$ and $\mathbf{U}_{\mathrm{n}}$

$t$ in $\left[0,T\right]$

Initialize global nodal force vector: $\mathbf{R}^{\mathrm{int}}=\left\{ \mathbf{0}\right\} $

Assemble $\mathbf{U}_{\mathrm{c}p}$ of $p$-th pattern using
Eq. (<ref>)

Calculate $\mathbf{R}_{p}^{\mathrm{ele}}$ of $p$-th pattern:

Assemble the $e$-th column of $\mathbf{R}_{p}^{\mathrm{ele}}$
into $\mathbf{R}^{\mathrm{int}}$

Calculate external force $\mathbf{R}_{t}^{\mathrm{ext}}$
at current time

Calculate next step displacement $\mathbf{U}_{\mathrm{n}}$
using Eq. (<ref>)

$\mathbf{U}_{\mathrm{p}}=\mathbf{U}_{\mathrm{c}}$; $\mathbf{U}_{\mathrm{c}}=\mathbf{U}_{\mathrm{n}}$

Explicit solver in serial

§ PARALLEL PROCESSING

Due to the fact that in explicit dynamics, the nodal force vector
can be assembled in an element-wise fashion, it is natural to exploit
parallel processing techniques to improve the efficiency of the solver.
The basic idea is to partition the mesh into different parts and assign
these parts to multiple processors. On each processor the nodal force
vector is calculated, and only the interface DOFs of each part need
to be synchronized. In this section, details regarding the partitioning
strategy and the parallel solver are presented.

§.§ Mesh partition

Mesh partition serves a vital role in parallel processing as the quality
of partition directly affects the efficiency of the solver. There
are mainly two types of mesh partitioning approaches, namely node-
and element-cut strategies [Krysl and Bittnar, 2001]. In this work, as most
of the computational effort is spent on the element processing stage,
the node-cut strategy is adopted. There are two main considerations
that must be taken into account in order to achieve a good partitioning
of the mesh:

* Size of the parts: An identical size of the individual parts ensures
that a balanced workload is achieved, i.e., all processors complete
their tasks at the same time. If that is not the case, some processors
are idle, and important resources are wasted as the speedup entirely
depends on the slowest processor. Although different patterns of elements
require a different amount of calculation, it is reasonable to estimate
the amount of work on each processor using the number of elements.

* Interface between parts: The interface between individual parts should
be as small as possible. While data communication between processors
cannot be avoided in explicit dynamics because the nodal force vectors
need to be synchronized in each time step, the target is to minimize
the amount of data communication which is the main bottleneck in parallel
processing in terms of the attainable speedup. The amount of data
that needs to be transferred is measured by the number of nodes on
the interface.

In order to partition a mesh, it is first converted into its dual
graph. Each element in the mesh is represented by a vertex $\mathrm{v}$
in the graph, while each interface shared by two elements is represented
by an edge $\mathrm{e}$ connecting two vertices. The set of vertices
is denoted as $V$ and the set of edges is referred to as $E$. A
typical quadtree mesh in 2D and its dual graph $G=\left\{ V,E\right\} $
are depicted in Fig. <ref>. An optimal partition can
be obtained by cutting through nodes $\left\{ \mathrm{n}_{4},\mathrm{n}_{12},\mathrm{n}_{11},\mathrm{n}_{19}\right\} $
and dividing the mesh into two parts with six elements each. Cutting
through $\left\{ \mathrm{n}_{3},\mathrm{n}_{8},\mathrm{n}_{11},\mathrm{n}_{12},\mathrm{n}_{15},\mathrm{n}_{20}\right\} $
is sub-optimal since more nodes are involved in the cutting, and thus
more results need to be communicated.

[A sample quadtree mesh][Dual graph and its optimal partition (shaded line)]

A quadtree mesh and its dual graph

However, in general, it is difficult to find the optimal partition,
which is known as an NP-hard problem (non-deterministic polynomial-time
hardness [Knuth, 1974]). There are several techniques in practice
that can provide a reasonably good approximation of the optimal partition,
including spectral partitioners and geometric separators. The spectral
partition method is based on the topology of the graph [Cvetkovic et al., 1980],
while the geometric separator is based on the coordinates of the vertices [Gilbert et al., 1998],
which can be represented by the geometric center of the elements in
an octree mesh. It is noted that spectral methods are often significantly
slower than geometric separators if the mesh size is large, as it
involves an eigenvalue decomposition. The algorithms presented here
only partition the mesh into two parts, but they can be applied recursively
to obtain partitions with an integer power of two parts.

Many mesh partition packages have been developed in recent years based
on these techniques [Karypis and Kumar, 1998, Ogawa, 2003, LaSalle and Karypis, 2016]. The
mesh partition package meshpart developed by Li, 2018
is used in this work, in which a spectral partitioner is used for
examples with less than 5 million DOFs, while a geometric separator
is used for all other cases.

§.§ Explicit solver in parallel

The explicit solver provided in Algorithm <ref>
can be easily parallelized based on the structure shown in Fig. <ref>.
The mesh generation and mesh partition are executed on the head node
(the main computer node which distributes jobs to individual processors),
and the element matrices of the 144 unique master cells are stored
as a static library, which is loaded every time the code runs. The
mesh partition package will generate an array of size $n_{\mathrm{e}}$,
which stores a flag for each element, indicating the part to which
it belongs. The element DOF vectors are grouped based on this array,
and they are assigned to the different processors. The displacement
vectors are initialized and distributed to the processors, where each
processor only stores the chunk associated with its own nodes. At
each time step, the nodal force vector is calculated using Eq. (<ref>)
on each processor locally. This local version of the nodal force vector
can be divided into two parts, the one corresponding to the nodes
inside each part (internal DOFs) and those shared with other parts
(external DOFs). The values of the latter part must be updated (corrected)
by the contributions of the adjacent part(s). The values of the external
nodes in the local nodal force vectors are summed and synchronized
for all the processors so that the contribution from all elements
connected to the interface is included. As the size of the interface
is relatively small compared to the total size, this operation can
be done in a relatively short time. After the updated nodal force
vector is distributed to each processor, the local portion of displacement
of the next step $\mathbf{U}_{n+1}$ is calculated.

Parallel code structure

At the beginning of the program, the octree mesh is generated and
partitioned. Pre-computed stiffness and mass matrices of the master
cells are imported and linked to the individual elements using corresponding
transformation matrices. The vector $\mathbf{V}$ is constructed and
distributed to individual processors based on the partitioning result.
The initial displacement vector $\mathbf{U}_{0}$ and fictitious $\mathbf{U}_{-1}$
are calculated and distributed to processors. These steps are considered
as an initialization stage, and are not included in the timing as
they take constant time independent of the number of time steps. The
speedup of the parallel solver is measured based on the time spent
after the initialization stage, including the calculation of the nodal
force vector, the synchronization between processors and the calculation
of the nodal displacement vector.

A sample timing of a parallel computation using 32 processors is shown
in Fig. <ref>. It can be observed that the
speed of these processors is usually different, as well as the amount
of work distributed to each processor. The light gray bar represents
the time spent by each processor to calculate the nodal force vector
locally, and the dark gray bar indicates the time for each processor
to wait for the slower ones to finish their job and communicate the
results. The black bar represents the time required to calculate the
displacement vector for the next time step once the overall nodal
force is obtained.

A sample of time distribution of parallel processing

In parallel computing, the accurate measurement of computational time
is difficult since different processors usually start working at different
time. Here we use the reasonable assumption that all processors commence
the calculation when the first (earliest) processor starts, and all
of them finish when the slowest processor completes the assigned task.
In the numerical examples, the average calculation time for all processors
(including light gray and black bars in Fig. <ref>)
and communication and waiting time (dark gray bars) are used to indicate
the performance of the proposed method.

§ NUMERICAL EXAMPLES

In this section, five numerical examples are presented to verify the
proposed procedure and to demonstrate the attainable speedup of the
novel parallel explicit solver with increasing number of computing
cores. The efficiency $\eta$ of the parallel computation is measured
by the ratio of the speedup and the number of cores (strong scaling).
The different transient amplitude functions that are used for our
numerical examples to describe the time-dependent excitations are
depicted in Fig. <ref>. Depending on the problem, we choose
either a Ricker wavelet (see Fig. <ref>), a triangular
pulse (see Fig. <ref>), or a sine-burst signal (see Fig. <ref>).
The analytical form of the Ricker wavelet in the time domain is

\begin{equation}
\end{equation}

where $P_{0}$ is the maximum load. The triangular impact load is
expressed as

\begin{equation}
\!\left(1-\!\left|\frac{t}{t_{1}}-1\right|\right)P_{0} & \textrm{if }0\leq t\leq2t_{1}\\
0 & \textrm{if }t>2t_{1}
\end{cases},\label{eq:Triangle}
\end{equation}

while the sine-burst (sine excitation modulated by a Hann window)
is expressed as

\begin{equation}
\sin\!\left(\frac{2\pi nt}{t_{1}}\right)\sin^{2}\!\left(\frac{\pi t}{t_{1}}\right)P_{0} & \textrm{if }0\leq t\leq t_{1}\\
0 & \textrm{if }t>t_{1}
\end{cases},\label{eq:Lamb}
\end{equation}

where $n$ is the number of cycles. The central frequency $f_{\mathrm{m}}$
is calculated as

\begin{equation}
\end{equation}

To obtain the frequency spectrum of the given excitation functions
we apply the Fourier transform. The frequency spectrum of the Ricker
wavelet is depicted in Fig. <ref> and can be
expressed as

\begin{equation}
\end{equation}

where the central frequency $f_{\mathrm{m}}$ is

\begin{equation}
f_{\mathrm{m}}=\frac{5}{\sqrt{2}\pi t_{1}}.\label{eq:Ricker_fm}
\end{equation}

Similarly, the amplitude spectrum of the Fourier transform of the
triangular pulse is shown in Fig <ref> and can
be written as

\begin{equation}
A\!\left(f\right)=\!\left(\frac{\sin\!\left(\pi ft_{1}\right)}{\pi ft_{1}}\right)^{2},\label{eq:Triangle_A}
\end{equation}

A critical frequency $f_{1}$ is identified such that the area under
the curve on the left of $f_{1}$ is equal to a certain percentage
of the total area, e.g., $95\%$. This approach ensures that most
of the energy content of the signal is accounted for. For a Ricker
wavelet, $f_{1}$ is approximately $2.2/t_{1}$, while in the case
of triangular impact $f_{1}$ is $1.8/t_{1}$. For the sine-burst
excitation, $f_{1}$ is identified corresponding to $85\%$ of the
area, which is $5.8/t_{1}$ when $n=5$. The determined critical frequency
is an important parameter that is closely related to the required
spatial discretization. Depending on the value of $f_{1}$ the wave
velocities of the propagating modes are calculated from which the
number of nodes per wavelength is derived.

[Load history of a Ricker wavelet][Spectrum of a Ricker wavelet]

[Load history of a triangular pulsed][Spectrum of a triangular pulse]

[Load history of a sine-burst signal][Spectrum of a sine-burst signal]

Dynamic loading for the examples

In infinite elastic media the only propagating waves are bulk waves
which can be distinguished in two types, i.e., dilatational/pressure/primary
and shear/secondary waves. The dilatational wave velocity $v_{\mathrm{P}}$
and shear wave velocity $v_{\mathrm{S}}$ for an isotropic, homogeneous,
linear-elastic medium are calculated as

\begin{align}
v_{\mathrm{P}} & =\sqrt{\frac{E\!\left(1-\nu\right)}{\rho\!\left(1+\nu\right)\!\left(1-2\nu\right)}},\\
v_{\mathrm{S}} & =\sqrt{\frac{E}{2\rho\!\left(1+\nu\right)}},
\end{align}

where $E$ is Young's modulus, $\nu$ denotes Poisson's
ratio, and $\rho$ stands for the mass density of the material. Therefore,
the minimum dilatational wavelength $L_{\mathrm{P}}$ and shear wavelength
$L_{\mathrm{S}}$ considered in the structure are

\begin{align}
L_{\mathrm{P}} & =\frac{v_{\mathrm{P}}}{f_{1}},\\
L_{\mathrm{S}} & =\frac{v_{\mathrm{S}}}{f_{1}}.
\end{align}

A common assumption on the spatial discretization
is that the maximum element size in the mesh should be limited to
have at least 10 (linear) finite elements per wavelength. This is
a bare minimum which will result in deviations of more than 1% compared
to a reference solution. Therefore, Willberg et al., 2012 and Gravenkamp, 2014
recommend roughly 30 nodes per wavelength when using linear shape
function in numerical methods for wave propagation analysis.

The in-house solver is written in Python 3.7, while the inter-processor
communication is implemented using MPI4Py [Dalcín et al., 2005, Dalcín et al., 2008, Dalcin et al., 2011],
which provides bindings for the Message Passing Interface (MPI) standard
with OpenMPI (4.0.2) as the back end. The distributed computing environment
used in this research is the supercomputer Gadi maintained by the
national computing infrastructure (NCI) in Canberra (Australia). A
typical compute node on Gadi has two sockets, each equipped with a
24-core Intel Xeon Scalable ‘Cascade Lake’
processor and connected to $190\,\unit{GB}$ RAM. The compute nodes
and the storage systems are connected through a high-speed network
with a data transfer rate of up to $200\,\unit{Gbps}$. In order to
maximize the parallel performance, the MPI processes are mapped by
sockets within each compute node, which provides an optimal achievable
memory bandwidth to each MPI process on both sockets [Chang et al., 2018].
For multi-node computations, the MPI processes are evenly distributed
among the compute nodes to fully exploit the computational capability
of Gadi.

§.§ Eigenfrequencies of a cube

In the first example, we verify the lumped mass matrix by investigating
the eigenfrequencies of a cube with a Young's modulus of $E=1000\,\unit{Pa}$
and a Poisson's ratio of $\nu=0.3$, while the mass density of the
cube is chosen as $\rho=1\,\unitfrac{kg}{m^{3}}$. The edge length
of the cube is $L=8\,\unit{m}$. The displacements perpendicular to
the surfaces are constrained as shown in 2D in Fig. <ref>.

Geometry and boundary conditions of the cube

Three types of meshes are generated based on three typical octree
patterns depicted in Fig. <ref>. Note that in the
coarsest mesh the edge length of an element is $1\,\unit{m}$. The
eigenfrequencies based on these meshes using either a consistent mass
matrix (CMM) or a lumped mass matrix (LMM) are calculated, as well
as the FEM result calculated using commercial software ABAQUS with
type 1 mesh. More detailed investigations regarding the use of LMMs
in the framework of the SBFEM are discussed in Ref. [Gravenkamp et al., 2020].
A convergence study of all element types is performed, and the $L_{2}$
norm of the relative error $e$ in the first 100 frequencies is presented
in Fig. <ref>. The error $e$ is calculated as

\begin{equation}
\end{equation}

where $\omega_{i}^{\mathrm{num}}$ is the numerical solution of the
$i$-th eigenfrequency, and $\omega_{i}^{\mathrm{ref}}$ is the analytical
(reference) solution. The analytical solution is given in the form

\begin{equation}
\omega^{\mathrm{ref}}=\sqrt{\frac{E\pi^{2}}{\rho\left(1-\nu^{2}\right)L^{2}}\!\left(l^{2}+m^{2}+n^{2}\right)\gamma},\label{eq:Cube_f}
\end{equation}

where $l$, $m$ and $n$ are non-negative integers, and the parameter
$\gamma$ is defined as

\begin{equation}
\gamma=\!\begin{cases}
1 & \textrm{if two of \ensuremath{\left\{  l,m,n\right\} } are 0 }\\
\!\left\{ 1,\frac{1-\nu}{2}\right\}  & \textrm{if one of \ensuremath{\left\{  l,m,n\right\} } is 0}\\
\!\left\{ 1,\frac{1-\nu}{2},\frac{1-\nu}{2}\right\}  & \textrm{if none of \ensuremath{\left\{  l,m,n\right\} } is 0}
\end{cases},\label{eq:Cube_gama}
\end{equation}

Note that if all of $\left\{ l,m,n\right\} $ are zeros, there is
no solution to the equation.

Three types of meshes of the cube

Relative error of first 100 eigenfrequencies calculated from different
types of meshes using different methods

The numerical results depicted in Fig. <ref> highlight
that with the proposed approach, an optimal convergence rate (2/3)
is achieved independent of using either an LMM- or CMM-based formulation
in the SBFEM or the FEM formulation. Furthermore, we observe that
the different octree pattern types have no influence on the slope
of the curves. Although the same rates of convergence are attainable,
we have to note that in general the CMM-based approach yields slightly
more accurate results for the same number of DOFs, and the FEM solution
is less accurate than both LMM- or CMM-based approach, which results
in an offset of the error curves. This difference is, however, of
no concern for the following discussions as optimal convergence is
still observed.

§.§ Wave propagation in a beam

In the second example, a cantilever beam is investigated, whose dimensions
are $16\,\unit{m}\times1\,\unit{m}\times1\,\unit{m}$. Its material
properties are: Young's modulus $E=10\,\unit{kPa}$, Poisson's ratio
$\nu=0$, and mass density $\rho=1\,\unitfrac{kg}{m^{3}}$. The beam
is fixed on one end, and a uniform pressure load is applied on the
other end, which follows the Ricker wavelet depicted in Fig. <ref>
with the parameters $t_{1}=15\,\unit{ms}$ and $P_{0}=1\,\unit{Pa}$.
Therefore, the maximum frequency of interest $f_{1}$ is $150\,\unit{Hz}$.
The dilatational wave speed, calculated using Eq. (<ref>),
is $100\,\unitfrac{m}{s}$, and the wavelength in the beam is $0.67\,\unit{m}$
according to Eq. (<ref>). This problem is equivalent
to a 1D wave propagation problem; therefore, Duhamel's integral is
applied to obtain a reference solution for the displacement response,
the details of which can be found in Ref. [Gravenkamp et al., 2020].

Geometry and boundary conditions of the beam

Again, three types of meshes are generated based on three typical
octree patterns as shown in Fig. <ref>, and the convergence
behavior is examined. In the coarsest group of meshes, the element
size is set to $0.0625\,\unit{m}$; as a result, there are approximately
10 elements in a wavelength. The time steps in the three coarsest
meshes are $0.5897\,\unit{ms}$, $0.3141\,\unit{ms}$ and $0.3380\,\unit{ms}$,
which guarantees there are at least 10 time steps[The Shannon-Nyquist theorem establishes a sufficient condition for
the sampling rate of a signal. According to this theorem, at least
two sampling points are required for one period of the signal associated
with the highest frequency. This is, however, not sufficient in structural
dynamics, and therefore, it is often suggested to use 10 to 20 sampling
points. For highly accurate simulations, this should be increased
to more than 100 sampling points.] during a period of the highest frequency. The displacement and acceleration
history of four selected points are plotted in Fig. <ref>.
The $L_{2}$ norm of the error of the selected points at all time
steps is calculated, and the convergence curves of these meshes are
presented in Fig. <ref>.

Three types of meshes of the beam

[Displacement history of four selected points][Acceleration history of four selected points]

Displacement and acceleration histories of four selected points on
the beam

[Convergence curve of displacement][Convergence curve of acceleration]

Convergence curves of different types of meshes

As in the modal analysis example of the previous section, we observe
that the theoretical rates of convergence for both displacement and
acceleration are recovered by our approach, independent of the chosen
mesh type. These results confirm the correctness of the implementation
and provide confidence in the performance of our method. Thus, more
complex examples that are of practical interest are tackled in the
following sections.

§.§ Multi-story building induced ground vibration

To further verify the proposed technique and evaluate its computational
performance, the ground vibration induced by a dynamic load applied
to a multi-story building is analyzed. The multi-story building and
the adjacent ground are depicted in Fig. <ref>. The
building consists of columns, beams, slabs and a foundation with a
basement, and the dimensions are illustrated in Fig. <ref>.

Geometry of the frame and ground

Dimension of the frame (Unit: mm)

The material properties of the frame are: Young's modulus $E=17\,\unit{GPa}$,
Poisson's ratio $\nu=0.2$ and mass density $\rho=2400\,\unitfrac{kg}{m^{3}}$.
The dilatational wave speed is equal to $2805\,\unitfrac{m}{s}$,
and the shear wave speed is $1718\,\unitfrac{m}{s}$. A region of
$400\,\unit{m}\times400\,\unit{m}\times200\,\unit{m}$ of the ground
adjacent to the building is included in the model. The coordinates
of point $P$ depicted in Fig. <ref> are $\left(-183\,\unit{m},-187\,\unit{m},-200\,\unit{m}\right)$,
and the origin of the coordinate system is denoted as $O$. The ground
is made of three layers, and the shapes of these interfaces between
the individual layers are defined as (unit: m)

\begin{align}
z_{12} & =-60+20\times\sin\!\left(\frac{x+183}{200}\pi\right),\\
z_{23} & =-130+20\times\cos\!\left(\frac{x+183}{200}\pi\right).
\end{align}

From the top to the bottom layers, the Young's
moduli are equal to $0.5\,\unit{GPa}$, $0.8\,\unit{GPa}$ and $1\,\unit{GPa}$,
respectively. All three layers have the same Poisson's ratio of $\nu=0.2$
and mass density of $\rho=2400\,\unitfrac{kg}{m^{3}}$. The corresponding
dilatational wave velocities are equal to $481\,\unitfrac{m}{s}$,
$609\,\unitfrac{m}{s}$ and $680\,\unitfrac{m}{s}$, and the shear
wave velocities are $295\,\unitfrac{m}{s}$, $373\,\unitfrac{m}{s}$
and $417\,\unitfrac{m}{s}$.

An area of $\unit[1]{m^{2}}$ at the side of the top of the building
($2.5\,\unit{m}\leq x\leq3.5\,\unit{m}$, $y=2.5\,\unit{m}$ and $39\,\unit{m}\leq z\leq40\,\unit{m}$),
indicated by region $A$ in Fig. <ref>, is subjected
to an impact load described by Fig. <ref> with the peak
pressure of $P_{0}=10^{7}\,\unit{Pa}$. The half duration of the impact
is $t_{1}=250\,\unit{ms}$, and the maximum frequency of interest
is $8\,\unit{Hz}$.

A reference solution is obtained from a convergence study using the
commercial finite element package ABAQUS. To limit the size of the
problem, only the building including its foundation with basement
is modeled. The displacements perpendicular to the side faces below
the ground level and the bottom of the basement are constrained. The
displacement and acceleration responses at Points $B$, $C$ and $D$
($B\!\left(2.5\,\unit{m},22.5\,\unit{m},40\,\unit{m}\right)$, corner
of the building; $C\!\left(22.5\,\unit{m},7.5\,\unit{m},40\,\unit{m}\right)$,
middle of a beam; $D\!\left(17\,\unit{m},17\,\unit{m},15\,\unit{m}\right)$,
center of a slab) indicated in Fig. <ref> are considered.
A series of analyses using uniform meshes of hexahedral elements (C3D8)
are performed. From the results of the convergence study we can infer
that an acceptable accuracy is obtained for meshes with an element
size of $0.0625\,\unit{m}$. Such a mesh contains 9,242,624 elements
and 32,186,598 DOFs.

Generating a non-uniform hexahedral mesh for an efficient simulation
of wave propagation in the complete system including both building
and ground requires considerable human interventions. To compare with
the present technique using octree meshes that are generated fully
automatically (shown as a part of Fig. <ref>), non-structured
tetrahedral meshes are used in ABAQUS. Both types of meshes are generated
according to specified element sizes. Using the building model and
the reference solution, we identify a tetrahedral mesh and an octree
mesh that yield a similar level of accuracy in the displacement and
acceleration response histories as plotted in Fig. <ref>.
The minimum and maximum element sizes are $0.0625\,\unit{m}$ and
$0.25\,\unit{m}$, respectively. Using our SBFEM-based code in conjunction
with an octree mesh, the number of elements in the building model
is 6,018,480 and therefore, the number of DOFs is 25,789,146. On the
other hand, the tetrahedral mesh generated in ABAQUS has a minimum
element size of $0.05\,\unit{m}$ while the maximum element size is
$0.25\,\unit{m}$. Hence, there are 27,274,592 elements and 21,632,154
DOFs in the mesh.

To model the ground, both the tetrahedral and octree meshes are generated
with a maximum element size of $4\,\unit{m}$, which is about 1/10
of the shear wavelength of $36.9\,\unit{m}$ in the top layer at the
maximum frequency of interest. There are 7,339,939 elements and 30,482,412
DOFs in the octree mesh of the whole system (Fig. <ref>),
while there are 38,847,472 elements and 26,192,571 DOFs in the tetrahedral
mesh. The time step used in this example is $\Delta t=0.018\,\unit{ms}$
which was estimated using Eq. (<ref>). Damping coefficient
in Eq. (<ref>) is chosen as $\alpha=0.1$. A total duration
of $1\,\unit{s}$ is analyzed with 55,556 steps.

Octree mesh of the frame and ground

Three points $E\!\left(0\,\unit{m},-30\,\unit{m},0\,\unit{m}\right)$,
$F\!\left(0\,\unit{m},-60\,\unit{m},0\,\unit{m}\right)$ and $G\!\left(0\,\unit{m},-90\,\unit{m},0\,\unit{m}\right)$
indicated in Fig. <ref> and located at a certain distance
from the building are selected for comparison purposes. The corresponding
displacement and acceleration response histories are plotted in Fig. <ref>.
It is observed that similar results are obtained with ABAQUS using
the tetrahedral mesh and with our in-house code utilizing the octree
mesh. The process of wave propagation in the ground is illustrated
in Fig. <ref> at four selected time instances.

[Displacement of points $B$, $C$ and $D$][Acceleration of points $B$]

[Acceleration of points $C$][Acceleration of points $D$]

The displacement and acceleration history of the building using our
in-house code and ABAQUS

[Displacement at points $E$, $F$ and $G$][Acceleration at point $E$]

[Acceleration at point $F$][Acceleration at point $G$]

Displacement and acceleration history at selected points on the ground
surface obtained with in-house code and ABAQUS

Contour of displacement magnitude in the ground at various time instances

Since the computational time increases linearly with the number of
steps of time integration, the CPU time required to complete 1000
steps is selected in measuring the computational performance. This
choice is made to facilitate the evaluation of the speedup of parallel
computing, starting from a small number of cores. The computational
times taken by the in-house code implementing the present technique
and by ABAQUS are listed in Table <ref> and are plotted
in a log-log scale in Fig. <ref>. To limit the
wall-clock time required by the ABAQUS analyses, we opt for using
four or more computer cores. It is observed that our in-house code
is significantly faster than the commercial software ABAQUS when using
a tetrahedral mesh to obtain an accuracy similar to our octree-based
solver. The speedup compared to a single core simulation and the efficiency
of the proposed method are shown in Fig. <ref>.

[Computational time taken by ABAQUS and the in-house code with different
number of cores][Speedup and efficiency of in-house code at different number of cores]

Comparison of computational time of the in-house code using octree
mesh and ABAQUS using tetrahedral mesh

Computational time, speedup and efficiency of in-house code and ABAQUS

2*$N_{\mathrm{cores}}$     ABAQUS     5c|in-house code     2*Ratio $t_{\mathrm{T\textrm{-ABAQUS}}}/t_{\mathrm{T\textrm{-in-house}}}$2-7 3-7 4-7 5-7 6-7 7-7
$t_{\mathrm{T}}$ [s]     $t_{\mathrm{C}}$ [s]     $t_{\mathrm{W}}$ [s]     $t_{\mathrm{T}}$ [s]     $s$     $\eta$     1     -     6565.47     -     6566.12     1.00     1.00     -2     -     3078.62     155.46     3234.37     2.03     1.02     -4     4550     1651.08     92.94     1744.20     3.76     0.94     2.618     2264     894.91     58.47     953.49     6.89     0.86     2.3716     1407     634.49     45.90     680.48     9.65     0.60     2.0732     1154     355.03     94.84     449.92     14.59     0.46     2.5664     534     141.43     63.28     204.73     32.07     0.50     2.61128     334     52.18     41.66     93.85     69.97     0.55     3.568l$t_{\mathrm{C}}$: Calculation time8l$t_{\mathrm{W}}$: Waiting and communication time8l$t_{\mathrm{T}}$: Total time8l$s$: Speedup factor8l$\eta$: Efficiency

We observe that with an increasing number of cores the computational
time is significantly reduced in both our in-house research code and
in ABAQUS. Additionally, it is noted that the proposed explicit solver
outperforms ABAQUS by a factor of at least two, meaning that for this
example only half of the computational time or less is required (Table <ref>).
The parallel efficiency for the frame example stays in a reasonable
range around 50% for 128 cores. For 128 cores, only around 230,000
DOFs are distributed to each core, which might be one source for the
decreasing efficiency with an increasing number of cores. The results
depicted in Fig. <ref> indicate that it is preferable
to partition the mesh such that more than one million DOFs are distributed
to each core. This is related to a more favorable ratio between time
spent on calculating nodal forces and displacements in contrast to
the communication and waiting time; therefore, the computational resources
can be used more efficiently.

§.§ Mountainous region

In this example, wave propagation in a mountainous region of size
$21.8\,\unit{km}\times18.0\,\unit{km}$ is investigated as an indication
of the practical applicability of the proposed parallel explicit solver.
The location of the region is shown in a satellite image from Google
Maps in Fig. <ref>, while the terrain information
(Fig. <ref>) and the STL model (Fig. <ref>)
are generated using the online service Terrain2STL at  [Chamberlin, 2020].
For illustration purposes, the model is assumed to be homogeneous
with a Young's modulus $E=55\,\unit{GPa}$, Poisson's ratio $\nu=0.2$
and mass density $\rho=2400\,\unitfrac{kg}{m^{3}}$. The dilatational
wave speed is $5046\,\unitfrac{m}{s}$, while the shear wave speed
is $3090\,\unitfrac{m}{s}$. The boundary conditions of the model
are shown in Fig. <ref>, where the displacements
perpendicular to the side faces and the bottom are constrained. An
excitation in the form of a Ricker wavelet (Fig. <ref>)
with $t_{1}=0.5\,{\textstyle \mathrm{s}}$ and $P_{0}=1.157\times10^{15}\,\unit{N}$
is applied at the center of the bottom of the model (point $A$) in
the $y$ direction. The maximum frequency of interest is $5\,\unit{Hz}$,
therefore, the wavelength of the shear wave is $618\,\unit{m}$.

[Satellite image ][Terrain ][STL model ]Geometry of the mountainous region

Boundary condition and dimensions of the mountainous region

An octree mesh is generated with a maximum element size of $68.125\,\unit{m}$
and a minimum element size of $8.516\,\unit{m}.$ Hence, the mesh
consists of 7,608,933 elements and 32,052,075 DOFs. The mesh is partitioned
using the algorithm proposed in Ref. [Gilbert et al., 1998]. The obtained
mesh partitioning is depicted in <ref>. The time
step is chosen as $\Delta t=1.38\,\unit{ms}$. The response is analyzed
for a duration of $10\,\unit{s}$ using 7246 time steps. The results
at the points $B\!\left(9000\,\unit{m},10000\,\unit{m},4000\,\unit{m}\right)$,
$E\!\left(11772\,\unit{m},9042\,\unit{m},7336\,\unit{m}\right)$ and
$F\!\left(9895\,\unit{m},4095\,\unit{m},4265\,\unit{m}\right)$ are
selected to evaluate the results and to illustrate the time-dependent
displacements. The origin of the coordinate system is point $O$ in
Fig. <ref>. The displacement and acceleration in
$y$ direction of these points are plotted against time in Fig. <ref>.
The contour of displacement magnitude at four time instances are presented
in Fig. <ref>, illustrating the wave propagation

[Displacements ][Accelerations ]

Displacements and accelerations in $y$ direction at five selected
points in the mountainous region model

Contour of displacement magnitude of the mountain region model at
four selected time instances $t$

The computational time and speedup for 1000 steps using different
numbers of cores are compared in Fig. <ref> and
Table <ref>. The presented results in terms of
the speedup factor and the parallel efficiency are in good agreement
with those reported for the frame example in Section <ref>.
We again observe that the efficiency of the parallel explicit solver
decreases when fewer DOFs are distributed to each core. However, in
this example, a higher parallel efficiency of roughly 60% is achieved
due to the shape of the model, which is more regular and allows for
a better partitioning.

[Computational time ][Speedup and efficiency ]

Computational time, speedup and efficiency of the mountainous region
simulation using different numbers of cores

Computational time, speedup and efficiency of the mountainous region

$N_{\mathrm{cores}}$     $t_{\mathrm{C}}$ [s]     $t_{\mathrm{W}}$ [s]     $t_{\mathrm{T}}$ [s]     $s$     $\eta$1     7273.36     -     7273.94     1.00     1.002     3650.19     46.03     3696.50     1.97     0.994     1788.73     36.26     1825.15     3.99     1.008     868.96     74.28     943.34     7.71     0.9616     504.91     65.58     570.55     12.7     0.7932     280.71     101.89     382.62     19.0     0.5964     136.02     47.31     183.34     39.7     0.62128     51.72     34.58     86.30     84.3     0.666l$t_{\mathrm{C}}$: Calculation time6l$t_{\mathrm{W}}$: Waiting and communication time6l$t_{\mathrm{T}}$: Total time6l$s$: Speedup factor6l$\eta$: Efficiency

§.§ Sandwich panel

In the last example, a sandwich panel with two steel cover sheets
and a foam-like aluminum core is modeled. The digital image, shown
in Fig. <ref>, is obtained by X-ray CT scans[The CT data of the aluminum foam core has been provided to the authors
by Prof. A. Düster and Dr. Meysam Joulaian [Joulaian, 2017, Mossaiby et al., 2019]
from Hamburg University of Technology.]. The size of the panel is $288\,\unit{mm}\times288\,\unit{mm}\times57.6\,\unit{mm}$,
while the thickness of the cover sheet is $4.8\,\unit{mm}$. Young's
modulus, Poisson's ratio and the mass density of aluminum are $E_{\mathrm{a}}=70\,\unit{GPa}$,
$\nu_{\mathrm{a}}=0.3$ and $\rho_{\mathrm{a}}=2700\,\unitfrac{kg}{m^{3}}$,
while the corresponding material parameters of steel are taken as
$E_{\mathrm{s}}=210\,\unit{GPa}$, $\nu_{\mathrm{s}}=0.3$ and $\rho_{\mathrm{s}}=8050\,\unitfrac{kg}{m^{3}}$.
An octree mesh is generated with a maximum element size of $0.96\,\unit{mm}$
and a minimum element size of $0.24\,\unit{mm}$, as shown in Fig. <ref>.
The mesh consists of 13,520,194 elements and 66,645,588 DOFs (for
comparison, a uniform voxel based mesh would contain 89,862,518 elements
and 288,049,833 DOFs). In order to test our approach for massively
parallel computing of large-scale wave propagation problems and examine
the increase of CPU time with the problem-size, we mirror the panel
in the $x$ and $y$ directions twice. Hence, the edge length of the
panel becomes $576\,\unit{mm}$ ($2\times2$ panel, Fig. <ref>)
and $1152\,\unit{mm}$ ($4\times4$ panel, Fig. <ref>),
respectively. In the $2\times2$ panel, there are 54,080,776 elements
and 266,295,885 DOFs, and in the $4\times4$ panel the number of elements
is 216,323,104, and the number of DOFs is 1,064,602,902.

In the simulation, we fix the displacement DOFs of the side faces
at $x=0$ and $y=0$ in orthogonal direction, as indicated in Fig. <ref>.
Additionally, Neumann boundary conditions are applied to approximate
the wave excitation that is usually achieved by means of piezoelectric
transducers being attached to the plate. To this end, two distributed
line loads, $P_{\mathrm{x}}$ and $P_{\mathrm{y}}$, are applied along
two line segments to simulate the excitation of a shear transducer.
This excitation is used to avoid the adverse effects of a concentrated
force and to achieve an excitation that is more realistic without
the need of including multi-physics elements with piezoelectric capabilities.
The fundamental idea for this excitation is based on the pin-force
model which is often adopted in the simulation of guided waves in
the context of structural health monitoring [Nieuwenhius et al., 2005, Greve et al., 2005, Giurgiutiu, 2008].
$P_{\mathrm{x}}$ is applied in $x$ direction on the line segment
($x=19.2\,\unit{mm}$, $0\,\unit{mm}\leq y\leq19.2\,\unit{mm}$ and
$z=57.6\,\unit{mm}$), while $P_{\mathrm{y}}$ is applied in $y$
direction, and the coordinate range is ($0\,\unit{mm}\leq x\leq19.2\,\unit{mm}$,
$y=19.2\,\unit{mm}$ and $z=57.6\,\unit{mm}$). These lines describe
the outer contour of a piezoelectric shear transducer which is often
applied in wave propagation analyses. Both of the distributed loads
vary over time following the sine-burst signal depicted in Fig. <ref>
with $t_{1}=0.050\,{\textstyle \mathrm{ms}}$, $n=5$ and $P_{0}=2.1\times10^{4}\,\unit{\unitfrac{N}{mm}}$.
In aluminum, the dilatational wave speed is $6198\,\unitfrac{m}{s}$,
and the shear wave speed is $3122\,\unitfrac{m}{s}$. In steel, the
dilatational wave speed is $5926\,\unitfrac{m}{s}$, and the shear
wave speed is $3168\,\unitfrac{m}{s}$. The center frequency is $100\,\unit{kHz}$
while the maximum frequency of interest is $116\,\unit{kHz}$, and
the shear wavelengths at the maximum frequency in the aluminum and
steel parts are $26.9\,\unit{mm}$ and $27.3\,\unit{mm}$, respectively.
The time step used in this example is $\Delta t=2.98\times10^{-5}\,{\textstyle \mathrm{ms}}$
according to Eq. (<ref>), and a total duration of $0.4\,{\textstyle \mathrm{ms}}$
is analyzed. Note that in the cover plate Lamb waves are excited which
will interact with the core and hence, energy will leak gradually
from the top plate into the bottom plate. For more detailed discussions
on guided waves and Lamb waves in particular the interested reader
is referred to Refs. [Viktorov, 1967, Rose, 1999, Giurgiutiu, 2008].

[Geometry of the sandwich panel]

[Mesh of the sandwich panel]

Geometry and mesh of the sandwich panel

[$2\times2$ panel]

[$4\times4$ panel]

Mirroring operation of the sandwich panel

Boundary condition of the sandwich panel with an enlarged view of
the loading area (left)

The mesh corresponding to the $4\times4$ panel can be naturally divided
into 16 parts by the mirror planes, and each part is further divided
into 1024 parts, as shown in <ref>. The top and bottom
views of the magnitude of the displacement response at five time instances
are presented in Fig. <ref>, illustrating the wave
propagation process. It can be observed that the wave in the top panel
disperses into the core and eventually to the bottom panel. The displacement
and acceleration of four selected points shown in Fig. <ref>
(coordinates are $D\!\left(288,0,57.6\right)$, $F\!\left(576,0,57.6\right)$,
$G\!\left(864,0,57.6\right)$ and $H\!\left(1152,0,57.6\right)$,
unit: mm) are plotted in Fig. <ref>. The speedup
of the multiple cores is illustrated in Fig. <ref>
and Table <ref>.

[Top view, $t=0.04\,\unit{ms}$][Bottom view, $t=0.04\,\unit{ms}$]

[Top view, $t=0.12\,\unit{ms}$][Bottom view, $t=0.12\,\unit{ms}$]

[Top view, $t=0.2\,\unit{ms}$][Bottom view, $t=0.2\,\unit{ms}$]

[Top view, $t=0.28\,\unit{ms}$][Bottom view, $t=0.28\,\unit{ms}$]

Contour of displacement magnitude at various time instances $t$

[Top view, $t=0.36\,\unit{ms}$][Bottom view, $t=0.36\,\unit{ms}$]

Contour of displacement magnitude at various time instances $t$

[Displacement of four points in the sandwich panel][Acceleration of four points in the sandwich panel]

Response history of selected points of the sandwich panel

[Computational time of $1\times1$ panel][Speedup and efficiency of $1\times1$ panel]

[Computational time of $2\times2$ panel][Speedup and efficiency of $2\times2$ panel]

[Computational time of $4\times4$ panel][Speedup and efficiency of $4\times4$ panel]

Time and speedup using multiple cores for the sandwich panel

Computational time of the sandwich panel

20.07$N_{\mathrm{cores}}$     3c|$1\times1$ panel     3c|$2\times2$ panel     3c|$4\times4$ panel2-10 3-10 4-10 5-10 6-10 7-10 8-10 9-10 10-10
$t_{\mathrm{C}}$ [s]     $t_{\mathrm{W}}$ [s]     $t_{\mathrm{T}}$ [s]     $t_{\mathrm{C}}$ [s]     $t_{\mathrm{W}}$ [s]     $t_{\mathrm{T}}$ [s]     $t_{\mathrm{C}}$ [s]     $t_{\mathrm{W}}$ [s]     $t_{\mathrm{T}}$ [s]4     3790.78     12.54     3803.67     -     -     -     -     -     -8     1937.86     36.19     1974.24     -     -     -     -     -     -16     1175.07     55.04     1230.24     5317.17     232.99     5550.67     -     -     -32     709.82     108.81     818.73     3253.78     514.92     3769.13     -     -     -64     283.48     81.35     364.87     1320.66     347.40     1668.21     6316.70     1780.27     8097.55128     119.05     41.04     160.11     620.23     204.36     824.66     2699.81     1010.44     3710.49256     44.49     38.07     82.57     258.58     115.88     374.49     1184.53     516.15     1700.79512     21.61     13.25     34.87     120.45     43.02     163.48     553.94     297.26     851.261024     9.92     6.46     16.37     47.45     37.05     84.63     264.50     115.35     379.872048     -     -     -     19.60     17.82     37.42     103.23     111.36     214.604096     -     -     -     9.81     6.92     16.73     42.76     45.70     88.468192     -     -     -     -     -     -     18.30     28.78     47.0816384     -     -     -     -     -     -     9.81     7.60     17.4110l$t_{\mathrm{C}}$: Calculation time10l$t_{\mathrm{W}}$: Waiting and communication time10l$t_{\mathrm{T}}$: Total time

Speedup and efficiency of the sandwich panel

20.07$N_{\mathrm{cores}}$     2c|$1\times1$ panel     2c|$2\times2$ panel     2c|$4\times4$ panel2-7 3-7 4-7 5-7 6-7 7-7
$s$     $\eta$     $s$     $\eta$     $s$     $\eta$4     4.00     1.00     -     -     -     -8     7.71     0.96     -     -     -     -16     12.37     0.77     16.00     1.00     -     -32     18.58     0.58     23.56     0.74     -     -64     41.70     0.65     53.24     0.83     64.00     1.00128     95.03     0.74     107.69     0.84     139.67     1.09256     184.27     0.72     237.15     0.93     304.71     1.19512     436.36     0.85     543.25     1.06     608.80     1.191024     929.21     0.91     1049.42     1.02     1364.27     1.332048     -     -     2373.29     1.16     2414.93     1.184096     -     -     5309.59     1.30     5858.25     1.438192     -     -     -     -     13028.57     1.5916384     -     -     -     -     29775.04     1.827l$s$: Speedup factor7l$\eta$: Efficiency

It can be observed that the CPU time increases with the number of
DOFs slightly more than linearly. For a fixed number of cores, better
efficiency can be achieved when the scale of the problem increases.
In the analysis of the $1\times1$ panel, a speedup factor of 929
is obtained when 1024 cores are used. In the $2\times2$ panel, the
maximum speedup is 5308 when using 4096 cores. In the $4\times4$
panel, a speedup factor of 29760 is obtained with 16384 cores. It
is noted that these values are calculated based on the assumption
that $100\%$ efficiency is achieved when using the lowest number
of cores, e.g., 4, 16 and 64, respectively. It is observed that for
a fixed problem scale, the efficiency can increase with the number
of cores in some cases, and some of the values are even higher than
1. This super-linear speedup is caused by the cache effect, i.e.,
as the mesh is partitioned and distributed to more cores, more data
can be fit into the cache memory of the compute nodes, which has significantly
higher data processing speed than RAM.

§ SUMMARY

In this article, a parallel explicit solver for structural dynamics
based on the CDM is developed, exploiting the inherent advantages
of octree meshes. In our implementation, the octree meshes are efficiently
handled by the SBFEM, which is our method of choice. However, it is
important to note that exactly the same approach can also be applied
to FEA when special transition elements with well-conditioned lumped
mass matrices are deployed. These could be based on transfinite mappings [Provatidis, 2019, Duczek et al., 2020],
virtual elements [Park et al., 2019], etc. A special mass lumping technique
is deployed by extending the work reported in Ref. [Gravenkamp et al., 2020]
to 3D, yielding a well-conditioned diagonal mass matrix. The advantages
of an explicit time integrator are leveraged, which enables the nodal
displacements to be computed without the need for solving a system
of linear equations. The pre-computation approach for the stiffness
and mass matrices, relying on a limited number of unique octree cells
in balanced meshes, further improves the efficiency of the algorithm.

The performance of the massively parallel explicit solver using the
octree mesh generation and the SBFEM is demonstrated by means of five
numerical examples, including academic benchmark problems with simple
geometry as well as image and STL based models. Large-scale problems
with up to one billion DOFs are simulated while up to 16,384 computing
cores are utilized. The octree algorithm provides a fully automatic
approach for mesh generation of complex geometric model. A Ricker
wavelet, a triangular pulse load and a sine-burst signal are considered
in these examples to excite the structures and to illustrate realistic
and high-resolution engineering applications. Compared to a conventional
(serial) implementation, a significant speedup is observed using the
parallelization approach on multi-processor machines, and the in-house
code can achieve an improved efficiency compared to the commercial
software ABAQUS. In the current approach, the imbalance of workload
of individual processors hinders a better scalability, as some processors
have to wait for the slower processors to finish the job so that the
results can be synchronized. In the future applications, a dynamic
partitioning approach and a load distribution of the job based on
the availability of individual processors may further improve the
performance of the proposed algorithm. In addition, non-linear analyses,
such as damage and contact problems can be studied making use of the
novel methodology.

§ ACKNOWLEDGMENTS

The work presented in this paper is partially supported by the Australian
Research Council through Grant Number DP180101538 and DP200103577.
This research includes results that have been obtained using the computational
cluster Katana supported by Research Technology Services at UNSW Sydney.
Furthermore, the provided resources from the National Computational
Infrastructure (NCI Australia), an NCRIS enabled capability supported
by the Australian Government, are gratefully acknowledged. The authors
would like to thank the NCI user support team for their technical
support regarding the use of NCI's supercomputing facility. The authors
would also like to thank Dr. Meysam Joulaian and Prof. Alexander
D�ster from Hamburg University of Technology for providing the X-ray
CT scan data, which was used in Section <ref>.

[Dupros et al., 2010]
authorF. Dupros, authorF. De Martin,
authorE. Foerster, authorD. Komatitsch,
authorJ. Roman,
titleHigh-performance finite-element simulations of
seismic wave propagation in three-dimensional nonlinear inelastic geological
journalParallel Computing volume36
(year2010) pages308–325.

[Brun et al., 2012]
authorM. Brun, authorA. Batti,
authorA. Limam, authorA. Gravouil,
titleExplicit/implicit multi-time step co-computations for
blast analyses on a reinforced concrete frame structure,
journalFinite Elements in Analysis and Design
volume52 (year2012) pages41–59.

[Bettinotti et al., 2017]
authorO. Bettinotti, authorO. Allix,
authorU. Perego, authorV. Oancea,
authorB. Malherbe,
titleSimulation of delamination under impact using a
global-local method in explicit dynamics,
journalFinite Elements in Analysis and Design
volume125 (year2017) pages1–13.

[Zhang and Chauhan, 2019]
authorJ. Zhang, authorS. Chauhan,
titleFast explicit dynamics finite element algorithm for
transient heat transfer,
journalInternational Journal of Thermal Sciences
volume139 (year2019) pages160–175.

[Meister et al., 2020]
authorF. Meister, authorT. Passerini,
authorV. Mihalef, authorA. Tuysuzoglu,
authorA. Maier, authorT. Mansi,
titleDeep learning acceleration of total Lagrangian
explicit dynamics for soft tissue mechanics,
journalComputer Methods in Applied Mechanics and
Engineering volume358 (year2020)

[Cook et al., 1989]
authorR. Cook, authorD. Malkus,
authorM. Plesha, titleConcepts and applications of
finite element analysis, publisherWiley, year1989.

[Patzák et al., 2001]
authorB. Patzák, authorD. Rypl,
authorZ. Bittnar,
titleParallel explicit finite element dynamics with
nonlocal constitutive models,
journalComputers & Structures volume79
(year2001) pages2287–2297.

[Talebi et al., 2012]
authorH. Talebi, authorC. Samaniego,
authorE. Samaniego, authorT. Rabczuk,
titleOn the numerical stability and mass-lumping schemes
for explicit enriched meshfree methods,
journalInternational Journal for Numerical Methods in
Engineering volume89 (year2012)

[Li et al., 2014]
authorB. Li, authorM. Stalzer,
authorM. Ortiz,
titleA massively parallel implementation of the optimal
transportation meshfree method for explicit solid dynamics,
journalInternational Journal for Numerical Methods in
Engineering volume100 (year2014)

[Park et al., 2019]
authorK. Park, authorH. Chi, authorG. H.
titleOn nonconvex meshes for elastodynamics using virtual
element methods with explicit time integration,
journalComputer Methods in Applied Mechanics and
Engineering volume356 (year2019)

[Anitescu et al., 2019]
authorC. Anitescu, authorC. Nguyen,
authorT. Rabczuk, authorX. Zhuang,
titleIsogeometric analysis for explicit elastodynamics
using a dual-basis diagonal mass formulation,
journalComputer Methods in Applied Mechanics and
Engineering volume346 (year2019)

[Rek and Němec, 2016]
authorV. Rek, authorI. Němec,
titleParallel computation on multicore processors using
explicit form of the finite element method and C++ standard libraries,
journalJournal of Mechanical Engineering
volume66 (year2016) pages67–78.

[Altman, 2014]
authorY. M. Altman, titleAccelerating MATLAB
Performance: 1001 tips to speed up MATLAB programs, publisherCRC
Press, year2014.

[Ma et al., 2020]
authorZ. Ma, authorY. Lou, authorJ. Li,
authorX. Jin,
titleAn explicit asynchronous step parallel computing
method for finite element analysis on multi-core clusters,
journalEngineering with Computers volume36
(year2020) pages443–453.

[Tang et al., 2020]
authorJ. Tang, authorY. Zheng,
authorC. Yang, authorW. Wang,
authorY. Luo,
titleParallelized implementation of the finite particle
method for explicit dynamics in GPU,
journalComputer Modeling in Engineering & Sciences
volume122 (year2020) pages5–31.

[Song and Wolf, 1997]
authorC. Song, authorJ. P. Wolf,
titleThe scaled boundary finite-element method–alias
consistent infinitesimal finite-element cell method–for elastodynamics,
journalComputer Methods in Applied Mechanics and
Engineering volume147 (year1997)

[Deeks and Wolf, 2002]
authorA. J. Deeks, authorJ. P. Wolf,
titleAn h-hierarchical adaptive procedure for the scaled
boundary finite-element method,
journalInternational Journal for Numerical Methods in
Engineering volume54 (year2002)

[Song, 2004]
authorC. Song,
titleA super-element for crack analysis in the time
journalInternational Journal for Numerical Methods in
Engineering volume61 (year2004)

[Bazyar and Song, 2008]
authorM. H. Bazyar, authorC. Song,
titleA continued-fraction-based high-order transmitting
boundary for wave propagation in unbounded domains of arbitrary geometry,
journalInternational Journal for Numerical Methods in
Engineering volume74 (year2008)

[Ooi et al., 2020]
authorE. T. Ooi, authorM. D. Iqbal,
authorC. Birk, authorS. Natarajan,
authorE. H. Ooi, authorC. Song,
titleA polygon scaled boundary finite element formulation
for transient coupled thermoelastic fracture problems,
journalEngineering Fracture Mechanics
(year2020) pages107300.

[Chiong et al., 2014]
authorI. Chiong, authorE. T. Ooi,
authorC. Song, authorF. Tin-Loi,
titleScaled boundary polygons with application to fracture
analysis of functionally graded materials,
journalInternational Journal for Numerical Methods in
Engineering volume98 (year2014)

[Gravenkamp and Natarajan, 2018]
authorH. Gravenkamp, authorS. Natarajan,
titleScaled boundary polygons for linear elastodynamics,
journalComputer Methods in Applied Mechanics and
Engineering volume333 (year2018)

[Talebi et al., 2016]
authorH. Talebi, authorA. A. Saputra,
authorC. Song,
titleStress analysis of 3D complex geometries using the
scaled boundary polyhedral finite elements,
journalComputational Mechanics volume58
(year2016) pages697–715.

[Saputra et al., 2017]
authorA. Saputra, authorH. Talebi,
authorD. Tran, authorC. Birk,
authorC. Song,
titleAutomatic image-based stress analysis by the scaled
boundary finite element method,
journalInternational Journal for Numerical Methods in
Engineering volume109 (year2017)

[Saputra et al., 2020]
authorA. A. Saputra, authorS. Duczek,
authorH. Gravenkamp, authorC. Song,
titleImage-based 3D homogenisation using the scaled
boundary finite element method,
journalComputers & Structures volume237
(year2020) pages106263.

[Gravenkamp et al., 2020]
authorH. Gravenkamp, authorA. A. Saputra,
authorS. Eisenträger,
titleThree-dimensional image-based modeling by combining
SBFEM and transfinite element shape functions,
journalComputational Mechanics  (year2020)

[Duczek et al., 2020]
authorS. Duczek, authorA. Saputra,
authorH. Gravenkamp,
titleHigh order transition elements: The xNy-element
concept-Part I: Statics,
journalComputer Methods in Applied Mechanics and
Engineering volume362 (year2020)

[Saputra et al., 2017]
authorA. A. Saputra, authorV. Sladek,
authorJ. Sladek, authorC. Song,
titleMicromechanics determination of effective material
coefficients of cement-based piezoelectric ceramic composites,
journalJournal of Intelligent Material Systems and
Structures volume29 (year2017)

[Gravenkamp and Duczek, 2017]
authorH. Gravenkamp, authorS. Duczek,
titleAutomatic image-based analyses using a coupled
quadtree-SBFEM/SCM approach,
journalComputational Mechanics volume60
(year2017) pages559–584.

[Liu et al., 2019]
authorL. Liu, authorJ. Zhang,
authorC. Song, authorC. Birk,
authorW. Gao,
titleAn automatic approach for the acoustic analysis of
three-dimensional bounded and unbounded domains by scaled boundary finite
element method,
journalInternational Journal of Mechanical Sciences
volume151 (year2019a)

[Liu et al., 2019]
authorL. Liu, authorJ. Zhang,
authorC. Song, authorC. Birk, authorA. A.
Saputra, authorW. Gao,
titleAutomatic three-dimensional acoustic-structure
interaction analysis using the scaled boundary finite element method,
journalJournal of Computational Physics
volume395 (year2019b)

[Xing et al., 2018]
authorW. Xing, authorC. Song,
authorF. Tin-Loi,
titleA scaled boundary finite element based node-to-node
scheme for 2D frictional contact problems,
journalComputer Methods in Applied Mechanics and
Engineering volume333 (year2018)

[Xing et al., 2019]
authorW. Xing, authorJ. Zhang,
authorC. Song, authorF. Tin-Loi,
titleA node-to-node scheme for three-dimensional contact
problems using the scaled boundary finite element method,
journalComputer Methods in Applied Mechanics and
Engineering volume347 (year2019)

[Zhang and Song, 2019]
authorJ. Zhang, authorC. Song,
titleA polytree based coupling method for non-matching
meshes in 3D,
journalComputer Methods in Applied Mechanics and
Engineering volume349 (year2019)

[Gravenkamp et al., 2012]
authorH. Gravenkamp, authorJ. Prager,
authorA. A. Saputra, authorC. Song,
titleThe simulation of Lamb waves in a cracked plate
using the scaled boundary finite element method,
journalThe Journal of the Acoustical Society of America
volume132 (year2012a)

[Gravenkamp et al., 2012]
authorH. Gravenkamp, authorC. Song,
authorJ. Prager,
titleA numerical approach for the computation of
dispersion relations for plate structures using the scaled boundary finite
element method,
journalJournal of Sound and Vibration
volume331 (year2012b)

[Gravenkamp et al., 2017]
authorH. Gravenkamp, authorA. A. Saputra,
authorC. Song, authorC. Birk,
titleEfficient wave propagation simulation on quadtree
meshes using SBFEM with reduced modal basis,
journalInternational Journal for Numerical Methods in
Engineering volume110 (year2017)

[Zou et al., 2019]
authorD. Zou, authorX. Teng, authorK. Chen,
authorJ. Liu,
titleA polyhedral scaled boundary finite element method
for three-dimensional dynamic analysis of saturated porous media,
journalEngineering Analysis with Boundary Elements
volume101 (year2019) pages343–359.

[Zhang et al., 2020]
authorJ. Zhang, authorJ. Eisenträger,
authorS. Duczek, authorC. Song,
titleDiscrete modeling of fiber reinforced composites
using the scaled boundary finite element method,
journalComposite Structures volume235
(year2020) pages111744.

[Liu et al., 2020]
authorL. Liu, authorJ. Zhang,
authorC. Song, authorK. He, authorA. A.
Saputra, authorW. Gao,
titleAutomatic scaled boundary finite element method for
three-dimensional elastoplastic analysis,
journalInternational Journal of Mechanical Sciences
volume171 (year2020) pages105374.

[Natarajan et al., 2020]
authorS. Natarajan, authorP. Dharmadhikari,
authorR. K. Annabattula, authorJ. Zhang,
authorE. T. Ooi, authorC. Song,
titleExtension of the scaled boundary finite element
method to treat implicitly defined interfaces without enrichment,
journalComputers & Structures volume229
(year2020) pages106159.

[Eisenträger et al., 2020]
authorJ. Eisenträger, authorJ. Zhang,
authorC. Song, authorS. Eisenträger,
titleAn SBFEM approach for rate-dependent inelasticity
with application to image-based analysis,
journalInternational Journal of Mechanical Sciences
volume182 (year2020) pages105778.

[Qu et al., 2020]
authorY. Qu, authorD. Zou, authorX. Kong,
authorX. Yu, authorK. Chen,
titleSeismic cracking evolution for anti-seepage face
slabs in concrete faced rockfill dams based on cohesive zone model in
explicit SBFEM-FEM frame,
journalSoil Dynamics and Earthquake Engineering
volume133 (year2020) pages106106.

[Liu et al., 2017]
authorY. Liu, authorA. A. Saputra,
authorJ. Wang, authorF. Tin-Loi,
authorC. Song,
titleAutomatic polyhedral mesh generation and scaled
boundary finite element analysis of STL models,
journalComputer Methods in Applied Mechanics and
Engineering volume313 (year2017)

[Song, 2018]
authorC. Song, titleThe Scaled Boundary Finite Element
Method: Introduction to Theory and Implementation, publisherJohn
Wiley & Sons, year2018.

[Zhang et al., 2020]
authorJ. Zhang, authorS. Natarajan,
authorE. T. Ooi, authorC. Song,
titleAdaptive analysis using scaled boundary finite
element method in 3D,
journalComputer Methods in Applied Mechanics and
Engineering volume372 (year2020)

[Gravenkamp et al., 2020]
authorH. Gravenkamp, authorC. Song,
authorJ. Zhang,
titleOn mass lumping and explicit dynamics in the scaled
boundary finite element method,
journalComputer Methods in Applied Mechanics and
Engineering volume370 (year2020)

[Chan et al., 1994]
authorT. F. Chan, authorJ. R. Gilbert,
authorS.-H. Teng,
titleGeometric spectral partitioning,

[Song, 2004]
authorC. Song,
titleA matrix function solution for the scaled boundary
finite-element equation in statics,
journalComputer Methods in Applied Mechanics and
Engineering volume193 (year2004)

[de Béjar and Danielson, 2016]
authorL. A. de Béjar, authorK. T. Danielson,
titleCritical time-step estimation for explicit
integration of dynamic higher-order finite-element formulations,
journalJournal of Engineering Mechanics
volume142 (year2016) pages04016043.

[Duczek and Gravenkamp, 2019]
authorS. Duczek, authorH. Gravenkamp,
titleMass lumping techniques in the spectral element
method: On the equivalence of the row-sum, nodal quadrature, and diagonal
scaling methods,
journalComputer Methods in Applied Mechanics and
Engineering volume353 (year2019a)

[Duczek and Gravenkamp, 2019]
authorS. Duczek, authorH. Gravenkamp,
titleCritical assessment of different mass lumping schemes
for higher order serendipity finite elements,
journalComputer Methods in Applied Mechanics and
Engineering volume350 (year2019b)

[Krysl and Bittnar, 2001]
authorP. Krysl, authorZ. Bittnar,
titleParallel explicit finite element solid dynamics with
domain decomposition and message passing: dual partitioning scalability,
journalComputers & Structures volume79
(year2001) pages345–360.

[Knuth, 1974]
authorD. E. Knuth,
titlePostscript about NP-hard problems,
journalSIGACT News volume6
(year1974) pages15–16.

[Cvetkovic et al., 1980]
authorD. M. Cvetkovic, authorM. Doob,
authorH. Sachs, titleSpectra of Graphs,
volume volume10, publisherAcademic Press, New York,

[Gilbert et al., 1998]
authorJ. R. Gilbert, authorG. L. Miller,
authorS.-H. Teng,
titleGeometric mesh partitioning: Implementation and
journalSIAM Journal on Scientific Computing
volume19 (year1998) pages2091–2110.

[Karypis and Kumar, 1998]
authorG. Karypis, authorV. Kumar,
titleA fast and high quality multilevel scheme for
partitioning irregular graphs,
journalSIAM Journal on scientific Computing
volume20 (year1998) pages359–392.

[Ogawa, 2003]
authorT. Ogawa,
titleParallelization of an adaptive Cartesian mesh flow
solver based on the 2N-tree data structure,
in: editorK. Matsuno, editorA. Ecer,
editorN. Satofuka, editorJ. Periaux,
editorP. Fox (Eds.), booktitleParallel Computational
Fluid Dynamics 2002, publisherNorth-Holland,
addressAmsterdam, year2003, pp.

[LaSalle and Karypis, 2016]
authorD. LaSalle, authorG. Karypis,
titleA parallel hill-climbing refinement algorithm for
graph partitioning,
in: booktitle2016 45th International Conference on
Parallel Processing (ICPP), organizationIEEE,
year2016, pp. pages236–241.

[Li, 2018]
authorY. Li, titlemeshpart,

[Willberg et al., 2012]
authorC. Willberg, authorS. Duczek,
authorJ. Vivar Perez, authorD. Schmicker,
authorU. Gabbert,
titleComparison of different higher order finite element
schemes for the simulation of Lamb waves,
journalComputer Methods in Applied Mechanics and
Engineering volume241 (year2012)

[Gravenkamp, 2014]
authorH. Gravenkamp, titleNumerical Methods for the
Simulation of Ultrasonic Guided Waves, volume volume116 of
publisherBAM-Berlin, year2014.

[Dalcín et al., 2005]
authorL. Dalcín, authorR. Paz,
authorM. Storti,
titleMPI for Python,
journalJournal of Parallel and Distributed Computing
volume65 (year2005) pages1108–1115.

[Dalcín et al., 2008]
authorL. Dalcín, authorR. Paz,
authorM. Storti, authorJ. DElía,
titleMPI for Python: Performance improvements and
MPI-2 extensions,
journalJournal of Parallel and Distributed Computing
volume68 (year2008) pages655–662.

[Dalcin et al., 2011]
authorL. D. Dalcin, authorR. R. Paz,
authorP. A. Kler, authorA. Cosimo,
titleParallel distributed computing using Python,
journalAdvances in Water Resources volume34
(year2011) pages1124–1139.

[Chang et al., 2018]
authorJ. Chang, authorK. Nakshatrala,
authorM. G. Knepley, authorL. Johnsson,
titleA performance spectrum for parallel computational
frameworks that solve PDEs,
journalConcurrency and Computation: Practice and
Experience volume30 (year2018)

[Chamberlin, 2020]
authorT. Chamberlin, titleTerrain2stl,

[Joulaian, 2017]
authorM. Joulaian, titleThe Hierarchical Finite Cell
Method for Problems in Structural Mechanics, publisherVDI
Fortschritt-Berichte Reihe 18 Nr. 348, year2017.

[Mossaiby et al., 2019]
authorF. Mossaiby, authorM. Joulaian,
authorA. Düster,
titleThe spectral cell method for wave propagation in
heterogeneous materials simulated on multiple GPUs and CPUs,
journalComputational Mechanics volume63
(year2019) pages805–819.

[Nieuwenhius et al., 2005]
authorJ. H. Nieuwenhius, authorJ. Neumann,
authorD. W. Greve, authorI. J. Oppenheim,
titleGeneration and detection of guided waves using PZT
wafer transducers,
journalIEEE Transactions on Ultrasonics, Ferroelectrics,
and Frequency Control volume52 (year2005)

[Greve et al., 2005]
authorD. W. Greve, authorJ. J. Neumann,
authorJ. H. Nieuwenhius, authorI. J. Oppenheim,
authorN. L. Tyson,
titleUse of Lamb waves to monitor plates: Experiments
and simulations,
in: booktitleProceedinggs of SPIE 5765,

[Giurgiutiu, 2008]
authorV. Giurgiutiu, titleStructural Health Monitoring
with Piezoelectric Active Wafer Sensors: Fundamentals and Applications,
publisherElsevier, year2008.

[Viktorov, 1967]
authorI. A. Viktorov, titleRayleigh and Lamb Waves,
Ultrasonic Technology, publisherPlenum Press.,

[Rose, 1999]
authorJ. L. Rose, titleUltrasonic Waves in Solid Media,
publisherCambridge University Press, year1999.

[Provatidis, 2019]
authorC. G. Provatidis, titlePrecursors of Isogeometric
Analysis, publisherSpringer International Publishing,

§ UNIQUE OCTREE PATTERNS IN 3D

The unique octree patterns with hanging nodes on zero to six edges
are listed. The remaining cases (seven to twelve hanging nodes) can
be obtained from these cases, e.g., the first case of an octree cell
with seven hanging nodes can be obtained from the first case with
five hanging nodes by swapping the edges with and without hanging
nodes, as shown in Fig. <ref>. The hanging nodes
at the center of faces are not shown, which will be added following
the discretization patterns.

Octree cell without hanging node

Octree cell with hanging node on one edge

Octree cells with hanging nodes on two edges

Octree cells with hanging nodes on three edges

Octree cells with hanging nodes on four edges

Octree cells with hanging nodes on five edges

Octree cells with hanging nodes on six edges

Octree cells with hanging nodes on seven edges (right) can be obtained
from those with five edges (left)

§ MESH PARTITION OF THE MOUNTAINOUS REGION

[two parts][four parts]

[eight parts][16 parts]

[32 parts][64 parts]

Mesh partition of the mountainous region

§ MESH PARTITION OF THE SANDWICH PANEL

[128 parts][256 parts]

[512 parts][1024 parts]

Mesh partition of the $1\times1$ sandwich panel