cross/.style=cross out, draw=black, minimum size=2*(#1-), inner sep=0pt, outer sep=0pt,cross/.default=1pt

§ ABSTRACT

Efficient solution strategies for wave propagation problems with complex geometries are essential in many engineering fields. Immersed boundary methods simplify mesh generation by embedding the domain of interest into an extended domain that is easy to mesh, introducing the challenge of dealing with cells that intersect the domain boundary. We consider the finite cell method that extends the weak form from the physical into the extended domain, scaled by a small value, to stabilize badly cut cells. Combined with explicit time integration schemes, the finite cell method introduces a lower bound for the critical time step size, which is restrictive for applications that require a small amount of stabilization not to compromise the accuracy. Explicit transient analyses commonly use the spectral element method due to its natural way of obtaining diagonal mass matrices through nodal lumping. The resulting method is very efficient since nodal lumping renders the solution of the equation systems trivial while not reducing the accuracy or the critical time step size. The combination of the spectral element and finite cell methods is called the spectral cell method. Unfortunately, a direct application of nodal lumping in the spectral cell method is impossible due to the special quadrature necessary to treat the discontinuous integrand inside the cut cells. The existing approaches to lump the mass matrices of cut cells significantly reduce the accuracy of the approximation.

We analyze an implicit-explicit (IMEX) time integration method to exploit the advantages of the nodal lumping scheme for uncut cells on one side and the unconditional stability of implicit time integration schemes for cut cells on the other. In this hybrid, immersed Newmark IMEX approach, we use explicit second-order central differences to integrate the uncut degrees of freedom that lead to a diagonal block in the mass matrix and an implicit trapezoidal Newmark method to integrate the remaining degrees of freedom (those supported by at least one cut cell). The immersed Newmark IMEX approach preserves the high-order convergence rates and the geometric flexibility of the finite cell method while retaining the efficiency of the nodal lumping scheme for uncut cells without compromising its critical time step size. We analyze a simple system of spring-coupled masses to highlight some of the essential characteristics of Newmark IMEX time integration. We then solve the scalar wave equation on two- and three-dimensional examples with significant geometric complexity to show that our approach is more efficient than state-of-the-art time integration schemes when comparing accuracy and runtime. While we focus on the finite and spectral cell methods, we expect other immersed methods to benefit equally from Newmark IMEX time integration.

implicit-explicit time integration, spectral elements, finite/spectral cell method, wave equation, immersed boundary methods

§ INTRODUCTION

Various technical and scientific applications depend on efficient methods to gain information about the propagation of waves. The most important fields include mechanical and civil engineering, medical engineering, and geophysics. For example, structural health monitoring (SHM) and non-destructive testing (NDT) require efficient and accurate simulations of waves propagating in plate-like structures or solids. In guided wave tomography [1, 2, 3, 4], Lamb waves are excited in thin plates for damage or corrosion detection. With its origins in seismic exploration [5, 6, 7], the full waveform inversion (FWI) employs simulated wave fields for the reconstruction of material models and anomalies in mechanical, civil, and medical engineering [8, 9, 10, 11, 12, 13]. The underlying mathematical models of these applications, among others, are most commonly one-, two-, or three-dimensional scalar or elastic wave equations.

Many combinations of spatial and temporal discretization methods exist to approximate wave equations. Popular spatial discretizations are the finite difference method (FDM) [14, 15] and variants of the finite element method (FEM) [16, 17], in particular, its $p$-extension [18] and isogeometric analysis (IGA) [19, 20] that uses B-Splines to generate approximations with higher continuity. For solutions with sufficient regularity, such high-order approaches benefit from higher convergence rates and a reduced dispersion error, enabling the solution of large-scale problems. A discussion of $p$- and $hp$-methods and their applications in electromagnetism can be found in [21]. Typical time integration approaches are basic finite difference schemes [22], the Newmark family [23], and Runge-Kutta methods [24, 25]. Particularly efficient approaches combine explicit time integration schemes with spatial discretizations that generate a diagonal mass matrix, as the solution of the equation systems of a time step becomes trivial. The spectral element method (SEM)

defines Lagrange polynomials of order $p$ on each element of a conforming mesh and combines them to $C^0$ continuous finite element basis functions. The SEM is predominantly used for explicit time integration, where lumping techniques diagonalize the mass matrix and thus render the solution of the equation system trivial. The widely adopted “nodal lumping” technique uses Gauss-Lobatto-Legendre (GLL) points as Lagrange nodes in the shape function definition and then again in the numerical quadrature of the finite element integrals to obtain a diagonal mass matrix due to the Kronecker delta property of the Lagrange shape functions. Nodal lumping does not suffer from a reduced accuracy compared to the consistent variant [26, 27]. The standard IGA does not naturally lump the mass matrix; however, it gains efficiency by exploiting the elevated continuity of the solution space across element interfaces to obtain higher accuracy with a reduced number of degrees of freedom (dofs) [28, 29]. Finding methods that diagonalize the mass matrix in IGA without losing much accuracy remains challenging. For instance, classical row-sum lumping leads to a breakdown in $p$-convergence [19], rendering it unsuitable for wave simulations. The most promising lumping approaches currently use approximate dual basis functions [30, 31].

Wave propagation problems often benefit from an automatic way of incorporating highly complex geometries and different types of geometric representations into the finite element simulation when a manual mesh generation is tedious or unfeasible. One possible remedy is to embed the computational geometry, called the physical domain, into an extended domain with a simple shape that can be meshed easily. Later, the influence of this newly added, fictitious domain onto the solution in the original, physical domain is eliminated during the assembly. The major challenges of such immersed boundary methods are the stabilization and numerical quadrature of cut elements (i.e., elements that intersect the domain boundary), the imposition of boundary conditions, and the solution of the (often still) ill-conditioned equation system. The finite cell method (FCM) is an immersed boundary method that retains the definiteness of the problem by defining the original physical model also in the fictitious domain but multiplied with a small number $\alpha_f$ to minimize the introduced model error [32, 18]. The extended domain is often chosen as the bounding box of the geometry and is meshed with a Cartesian grid whose finite elements are then called finite cells. The quadrature of the finite cell integrals on the cut cells is challenging as their discontinuous nature renders standard Gauss-Legendre rules unsuitable. A robust (yet expensive) strategy is to create a spacetree-based subdivision of the cut cells that is refined towards the domain boundary and distribute Gauss-Legendre points on each of these newly created integration cells. Alternatively, the computation of custom weights for a fixed set of quadrature points via moment fitting can drastically reduce the assembly effort [33].

The combination of the spectral element and finite cell methods is called the spectral cell method (SCM). Like the SEM, the SCM uses Lagrange polynomials with GLL nodes; however, the standard quadrature methods to address the discontinuity of the finite cell integrals introduce non-diagonal blocks in the mass matrix. Any additional quadrature point that is not a Lagrange node (GLL point) breaks the natural lumping as it adds off-diagonal terms due to the interaction of multiple nonzero basis functions. We call this approach the consistent SCM [34]. The HRZ-lumped SCM [35, 36] uses HRZ lumping to diagonalize the mass matrices of cut cells. Alternatively, a moment fitting with fixed GLL points recovers the nodal lumping by computing custom weights for the GLL points [37, 38]. Unfortunately, lumped SCM methods suffer from reduced accuracy [39], and we observe a significantly restricted critical time step size in the presence of badly cut cells.

As a result, the cell with the lowest volume fraction dominates the stability behavior of the entire system. As shown in [40], eigenvalue stabilization can mitigate the severe limitation of the critical time step size. The development of explicit versions of the isogeometric finite cell analysis (IGA-FCM) [41, 42] that combine IGA with the FCM faces similar challenges. IGA-FCM has been applied to various transient problems [43, 44, 45, 13], and while the critical time step of its lumped version appears to be independent of the cut cells, the reduction in accuracy is severe [44, 45]. An in-depth discussion of the critical time step size of immersed boundary methods, particularly in the context of IGA-FCM, can be found in [46].

A similar restriction on the overall time step size exists in dynamic fluid-structure and soil-structure problems [47]. This limitation arises because the structure mesh introduces a stability limit significantly lower than the medium mesh when using explicit time integration. To address this challenge, Belytschko and Mullen [48] and Hughes and Liu [49] suggest integrating the part of the mesh belonging to the structure implicitly in time and the part of the mesh belonging to the fluid or soil explicitly. This combination results in an implicit-explicit (IMEX) time integration scheme. The analysis in [49, 50] shows that the IMEX integration approach inherits the stability characteristics of the explicit subsystem while preserving the accuracy of both individual methods. We transfer this scheme to systems discretized using the SCM. The part of the system corresponding to cut cells — hence strongly restricting the stability and contributing to a non-diagonal block in the global mass matrix — is integrated implicitly over time using the trapezoidal Newmark method [51]. In contrast, the part of the system corresponding to uncut cells — thus contributing to a diagonal block in the global mass matrix — is integrated explicitly using the second-order central difference method (CDM) [17]. The proposed approach, denoted as immersed Newmark IMEX method, preserves the accuracy of the consistent SCM, eliminates the dependency of the time step size on the cut cells, and exploits the efficiency of explicit time stepping for the diagonal block of the mass matrix belonging to uncut cells.

The central aspects of our paper are as follows.

* We discuss the Newmark IMEX time integration and examine its stability and accuracy on a simplified system of ten spring-coupled masses.

* We apply this approach to the spectral cell method, where we integrate the dofs of cut cells implicitly and the remaining dofs explicitly. We evaluate the efficiency and accuracy on a two-dimensional perforated plate and compare it to conventional time integration schemes.

* We demonstrate the applicability of the immersed Newmark IMEX method to large-scale problems using a complex three-dimensional example. Our approach performs better, even when dealing with an unfavorable ratio between the implicit and explicit subsystems.

The remainder of the paper is structured as follows. In sec:theory, we start with a discussion of the critical time step size to motivate the development of our immersed Newmark IMEX approach. We then introduce the scalar wave equation and outline its spatial and temporal discretization. sec:NumericExamples applies the proposed approach to examples of increasing complexity. First, we demonstrate the convergence and stability characteristics of the implicit-explicit time integration scheme using a system of spring-coupled masses. Second, we compare the efficiency of the immersed Newmark IMEX method to conventional time integration schemes on a two-dimensional example. Third, we validate the improved performance on a three-dimensional problem. Finally, sec:conclusion concludes our paper and highlights some of the limitations.

§ IMEX FOR IMMERSED BOUNDARY METHODS

§.§ Critical time step analysis of cut cells

Explicit time integration is computationally efficient but comes with the drawback of conditional stability. A critical time step size $\Delta t_\text{crit}$ defines the stability limit, beyond which the time integration becomes unstable, and the solution grows infinitely [17]. In the context of the spectral cell method (SCM), the critical time step size of the global system is predominantly determined by cells intersected by the boundaries of the geometry, even if the number of cut cells is small. For the explicit second-order central difference method (CDM), the systems's critical time step size $t_\text{crit}$ is equal to [17]

\begin{equation}
\Delta t_\text{crit} = \frac{2}{\omega_\text{max}}\text{,}
\end{equation}

where $\omega_\text{max}$ is the highest eigenfrequency of the system. For a spatially discretized system, we obtain

\begin{equation}
\omega_\text{max} = \sqrt{\lambda_\text{max}(\v{K}, \v{M})}
\end{equation}

by computing the largest eigenvalue $\lambda_\text{max}$ of the generalized eigenproblem corresponding to the stiffness matrix $\v{K}$ and mass matrix $\v{M}$.

When dealing with cut cells in the SCM, the maximum eigenfrequency of one cut cell depends on the ratio between its physical and fictitious parts. fig:fillRatio shows a motivating example that we use to analyze the behavior of the eigenfrequency for badly cut cells.

Cell cut by a horizontal boundary

A horizontal boundary with varying heights $d_\text{p}$ cuts a single quadrilateral cell. Below the cut, the material's density and wave speed are constant, i.e., $\rho_\text{p} = c_\text{p} = 1$. Above the boundary, the material is fictitious, and its density is scaled to $\rho_\text{f} = 10^{-6}$, while the wave speed remains untouched, i.e., $c_\text{f} = 1$. The fill ratio of the cell, depending on the cut position, is $\eta = \frac{d_\text{p}}{l_\text{y}}$.
The maximum eigenfrequency of the cell is computed using spectral basis functions of polynomial orders $p = 1, 2, $ and $3$, and a quadtree of depth $13$ for the integration of the mass and stiffness matrices. fig:fillRatioHZ plots the maximum eigenfrequencies of the cells compared to their fill ratios. Cells with smaller fill ratios feature higher maximum eigenfrequencies, which leads to a decisive restriction on the critical time step size. This effect is more pronounced the higher the polynomial degree of a cell.

Highest eigenfrequency for different fill ratios $\eta$ and polynomial degrees $p$

In addition to the cut ratio, factors such as the cut's shape also significantly affect the maximum eigenfrequency of a cut cell. A detailed discussion regarding critical time step sizes of immersed methods can be found in [46]. Implicit-explicit time integration overcomes this limitation of the SCM by integrating cut cells implicitly in time.

§.§ Spatial Discretization

For a given spatial dimension $d = 1, 2,$ or $3$, the scalar wave equation is defined on a domain $\Omega \subset \mathbb{R}^d$ for the time interval $\mathcal{T} = \left[0, T\right]$, where $T$ is the simulation duration. By introducing the wave field $u(\v{x}, t))$, its acceleration $\ddot{u}(\v{x}, t))$, the external force term $f(\v{x}, t))$, the heterogeneous density $\rho(\v{x})$ and wave speed $c(\v{x})$, we obtain the following partial differential equation (PDE):

\begin{equation}
\rho(\v{x}) \ddot{u}(x, t) - \nabla \cdot (\rho(\v{x}) c^2(\v{x}) \nabla u(\v{x}, t)) = f(\v{x}, t), \;
\v{x} \in \Omega, \quad t \in \mathcal{T}\text{.}
\label{eq:scalarWaveEq}
\end{equation}

The initial conditions are $u(\v{x}, 0) = u_\text{0}(\v{x})$ and $\dot{u}(\v{x}, 0) = \dot{u}_\text{0}(\v{x})$ on $\v{x} \in \Omega$. The Dirichlet boundary conditions are $u(\v{x}, t) = g(\v{x}, t)$ on $\v{x} \in \Gamma_\text{D}$, and the Neumann boundary conditions are ${\v{n} \cdot \nabla u(\v{x}, t) = h(\v{x}, t)}$ on $\v{x} \in \Gamma_\text{N}$, where $\partial \Omega = \Gamma = \Gamma_\text{N} \cup \Gamma_\text{D}$ and $\Gamma_\text{N} \cap \Gamma_\text{D} = \emptyset$ is satisfied. The unit normal vector $\v{n}$ on the boundary $\partial \Omega$ points outwards.

The finite element discretization introduces $n^\text{dof}$ basis functions $N_i(\v{x})$ and their corresponding time-dependent coefficients $\hat{u}_i(t)$, for $i = 1, 2, \ldots, n^\text{dof}$ to approximate the wave field as

\begin{equation}
u(\v{x}, t) \approx \sum_{i = 1}^{n^\text{dof}} N_i(\v{x}) \hat{u}_i(t) \text{.}
\end{equation}

The hat notation indicates the coefficients of a discretized quantity. An in-depth introduction of the finite element method (FEM) for transient problems can be found in [17]. Like other finite element methods, the SEM subdivides the computational domain into a set of elements that support tensor-products of Lagrange polynomials of degree $p$. These shape functions are connected with those of adjacent elements to form the basis functions $N_i(\v{x})$. The $p + 1$ Lagrange polynomials in the reference interval $[-1, 1]$ are defined as follows

\begin{equation}
N_i^\text{Lag,p}(\xi) = \prod_{j = 1, j \neq i}^{p + 1} \frac{\xi - \xi_j}{\xi_i - \xi_j} \quad \text{for } i = 1, 2, \ldots, p + 1 \text{.}
\end{equation}

We choose Gauß-Lobatto-Legendre (GLL) points for the definition of $N_i^\text{Lag,p}$ to obtain a diagonal mass matrix. We refer to [26, 27, 7] for a discussion of GLL points combined with the SEM. The $p + 1$ GLL points of one element are given by

\begin{equation}
\xi_j = \begin{cases}
-1 \qquad &\text{if } j = 1 \\
\xi_{0, j-1}^{\text{Lo}, p-1} \qquad &\text{if } 2 \leq j < p + 1 \\
1 \qquad &\text{if } j = p+1 \\
\end{cases}.
\end{equation}

where $\xi_{0, k}^{\text{Lo}, p-1}$ with $k = 1, 2, \ldots, p-1$ are the roots of Lobatto polynomials of order $p-1$. Combining the SEM basis functions with an immersed boundary approach, in particular, the finite cell method (FCM) [32, 18], leads to the spectral cell method (SCM) [34, 35], a non-boundary-conforming adaption of the SEM. In both FCM and SCM, the physical domain $\Omega_\text{p}$ is embedded within a larger computational domain $\Omega$, which is of simple geometric shape and can be readily meshed using a Cartesian grid. As illustrated in fig:FCM, the computational domain $\Omega$ is subdivided into the physical domain $\Omega_\text{p}$ and the fictitious domain $\Omega_\text{f}$, such that

\begin{equation}
\Omega_\text{p} \cup \Omega_\text{f} = \Omega \text{.}
\end{equation}

In the context of FCM, the elements generated by the Cartesian grid are referred to as cells, and those cells that are intersected by the boundaries of the physical domain are termed cut cells.

[Physical domain ]
[Meshed embedding domain ]
[Integration mesh defined by a quadtree ]
Principles of the FCM/SCM [39]

We introduce an indicator function $\alpha(\v{x})$ to recover the original geometry and differentiate between the physical domain $\Omega_\text{p}$ and the fictitious domain $\Omega_\text{f}$. Within the fictitious domain, $\alpha(\v{x})$ is set to a value close to zero, while it remains one in the physical domain

\begin{equation}
\alpha(\v{x}) = \begin{cases}
1 \qquad &\v{x} \in \Omega_\text{p} \\
\alpha_\text{f} = 10^{-\beta} \qquad &\v{x} \in \Omega_\text{f}
\end{cases}.
\end{equation}

In practice, $\beta$ is commonly assigned some value between five and eight. According to [34, 35], the entire weak form of the scalar wave equation is multiplied by $\alpha(\v{x})$. This scaling is investigated in [12] is equivalent to scaling the system's density, while keeping the wave speed untouched, i.e.,

\begin{equation}
\rho^*(\v{x}) = \alpha(\v{x}) \rho(\v{x}) \text{.}
\end{equation}

This scaling turns out to be most favourable for treating voids within the computational domain. A generalization to the elastic wave equation and a detailed mathematical analysis is given in [52]. The entries of the corresponding mass matrix $\v{M}$, stiffness matrix $\v{K}$ and time-dependent force vector $\v{f}(t)$ are obtained by

M_ij    = ∫_Ω α(x̌) ρ(x̌) N_i(x̌) N_j(x̌) dΩ,

K_ij    = ∫_Ω α(x̌) ρ(x̌) c(x̌)^2 (∇N_i(x̌))^T ∇N_j(x̌) dΩ,

f̂_i(t)    = ∫_Ω α(x̌) f(x̌, t) N_i(x̌) dΩ.

This spatial discretization yields the second-order system

\begin{equation}
\v{M} \ddot{\dv{u}}(t) + \v{K} \dv{u}(t) = \dv{f}(t) \text{,}
\label{eq:discretizedWaveEquation}
\end{equation}

where the time-dependent coefficient vector $\dv{u}(t)$ contains the unknowns of the system.
In wave propagation applications, it is often possible to split the force term into a spatial contribution and a temporal excitation

\begin{equation}
\dv{f}(t) = f_\text{t}(t) \cdot \dv{f}_{\text{x}}.
\end{equation}

The spatial force vector $\dv{f}_{\text{x}}$ entails the geometric information of the excitation, e.g., the source positions. The term $f_\text{t}(t)$ contains the excitation signal in time.

In the finite and spectral cell methods, the boundaries of the physical domain are recovered by accurately computing the discontinuous integrands in Equations (<ref>), (<ref>) and (<ref>). There are several approaches for integrating cut cells, among them are: space-trees [32, 53], smart octrees [54], and moment-fitting [33, 37]. We use space-trees (i.e. binary-, quad-, and octrees in one, two, and three dimension) that are refined towards the domain boundary for the integration of cut cells. fig:fcmQuadtree shows an exemplary integration mesh for a quadtree of depth $3$. The leaves of the spacetree are used to accumulate Gauss-Legendre points close to the domain boundary.

In the SCM, the choice of a cell's quadrature rule depends on whether a cell is cut or not. Cells not intersected by the boundaries of the physical domain (from now on referred to as uncut cells) are integrated using the Gauß-Lobatto quadrature with GLL points as integration points. The latter causes the element mass matrices to be diagonal for uncut elements, since in each uncut element at all GLL points exactly one Lagrange basis function is unequal to zero. However, if a cell is cut, we use a space-tree integration scheme. Since the diagonality of the element mass matrix is lost on cut elements, we prefer to distribute the more accurate Gauß-Legendre quadrature on the spacetree integration mesh. The element stiffness matrices and the load vectors are integrated using Gauß-Legendre quadrature regardless of whether the cell is cut or uncut.

§.§ Partitioning of the Degrees of Freedom

The degrees of freedom of basis functions supported exclusively in uncut cells contribute to exactly one entry in the diagonal mass matrix block $\v{M}^\text{dd}$. Conversely, dofs with support in at least one cut cell contribute to a non-diagonal block of the mass matrix $\v{M}^\text{cc}$. We partition the mass matrix into diagonal and a non-diagonal blocks

\begin{equation}
\label{eq:MStructure}
\v{M} = \begin{bmatrix} \v{M}^\text{dd} & \v{0} \\ \v{0} & \v{M}^\text{cc} \end{bmatrix} \text{.}
\end{equation}

In the following, the dofs, the stiffness matrix, and the force vector are partitioned based on this distinction. The off-diagonal blocks in the mass matrix equal $\v{0}$ and do not introduce any coupling. Sub-vectors or sub-matrices corresponding to the dofs with the diagonal mass matrix block are denoted by the superscript $\star^\text{d}$ for `diagonal', while all other dofs corresponding to the non-diagonal block of the mass matrix are denoted by $\star^\text{c}$ for `cut'. fig:MStructure illustrates an exemplary structure of an obtained global mass and stiffness matrix.

[Mass matrix $\v{M}$ ]
[Stiffness matrix $\v{K}$ ]
Typical matrix structures for a consistent SCM discretization

The stiffness matrix $\v{K}$ couples diagonal and cut dofs, as demonstrated in fig:KStructure. Nevertheless, it is possible to decompose the stiffness matrix into

\begin{equation}
\label{eq:KStructure}
\v{K} = \begin{bmatrix} \v{K}^\text{dd} & \v{K}^\text{dc} \\ \v{K}^\text{cd} & \v{K}^\text{cc} 	\end{bmatrix}.
\end{equation}

Given the symmetric nature of the stiffness matrix, $\v{K}^\text{dc} = (\v{K}^\text{cd})^T$ holds. For future use, the sub-matrices $\v{K}^\text{d}$ and $\v{K}^\text{c}$ are introduced

Ǩ^d    =  Ǩ^dd     Ǩ^dc

Ǩ^c    =  Ǩ^cd     Ǩ^cc 	.

The spatial force vector $\dv{f}$ is partitioned accordingly

\begin{equation}
\dv{f}_\text{x} = \begin{bmatrix} \bigl(\dv{f}_\text{x}^\text{d}\bigr)^T, & \bigl(\dv{f}_\text{x}^\text{c}\bigr)^T \end{bmatrix}^T \text{.}
\label{eq:f_xStructure}
\end{equation}

§.§ Time Integration

Flowchart for the trapezoidal Newmark method adapted from [55]

Flowchart for the Newmark IMEX method

Since numerical time integration exclusively deals with discretized quantities, the $\hat{\star}$ notation is omitted in this section.
The primary advantage of explicit methods is their potentially higher computational efficiency per time step compared to implicit methods. On the other hand, they are only conditionally stable. However, as long as the time step size limitation is not too restrictive, explicit methods, particularly the CDM, offer a good trade-off between accuracy and computational cost [17]. The following derivations closely follow [39]. Newmark's time integration formulas [23] are well established and equal

ǔ_n+1    = ǔ_n + Δt ǔ̇_n + Δt^2 ((1/2 - β) ǔ̈_n + βǔ̈_n+1)

ǔ̇_n+1    = ǔ̇_n + Δt ((1-γ) ǔ̈_n + γǔ̈_n+1) .

The two parameters $\beta$ and $\gamma$ serve to define the weights for the acceleration calculation, $\beta$ for the displacement estimation, and $\gamma$ for the velocity estimation, respectively. Depending on the choice of these parameters, the method is either explicit or implicit. The explicit CDM, is obtained by setting $\beta = 0$ and $\gamma = \frac{1}{2}$. The formulas then simplify to

ǔ_n+1    = ǔ_n + Δt ǔ̇_n + Δt^2/2 ǔ̈_n

ǔ̇_n+1    = ǔ̇_n + Δt/2(ǔ̈_n + ǔ̈_n+1) .

The CDM is second-order accurate and exhibits no numerical dissipation [49]. In the context of wave propagation, one can often omit the computation of the velocities $\dot{\v{u}}$. As a consequence, the velocity can be eliminated by substitution, resulting in a two-step scheme, where

\begin{equation}
\v{u}_{n+1} = 2 \v{u}_n - \v{u}_{n-1} + \Delta t^2 \ddot{\v{u}}_n \text{.}
\label{eq:CDMstep}
\end{equation}

is referred to as the CDM sum step. The accelerations $\ddot{\v{u}}$ are computed from the spatially discretized equations of motion

\begin{equation}
\v{M} \ddot{\v{u}}_n = f_\text{t}(t_n) \v{f}_\text{x} - \v{K} \v{u}_n \text{.}
\label{eq:CDMacc}
\end{equation}

The computation of the acceleration together with the CDM sum step will hereafter be called CDM step. Since the mass matrix remains constant, it can be factorized once outside the time integration loop, e.g., with a Cholesky-decomposition. Inside the loop, the system of equations is solved using backward and forward substitution. For a diagonal mass matrix, this factorization is omitted since the system of equations can directly be solved by element-wise inversion of the diagonal mass matrix. The CDM is extensively used, e.g., in [35, 44, 12, 13].

The CDM is highly efficient for uncut cells due to its moderate time step size constraint and the low computational cost per time step. However, the situation changes in the presence of cut cells, especially for small cuts. The very small critical time step size for explicit integration of cut cells significantly slows down the overall computation process. One way to overcome this limitation is to introduce substepping on cut cells, commonly referred to as leap-frog method [56, 37]. By introducing a ratio $m$, the cut cells are integrated in time with a finer time step size $\Delta t_{\text{fine}}$:

\begin{equation}
m = \frac{\Delta t_{\text{coarse}}}{\Delta t_{\text{fine}}} \text{.}
\end{equation}

The parameter $\Delta t_{\text{coarse}}$ denotes the coarse time step size of the overall system. In this way, only $\Delta t_{\text{fine}}$ must satisfy the highly restrictive stability constraints for CDM of cut cells. Implementations details of the leap-frog method applied to the SCM can be found in [37, 57].

Another way to overcome the restrictive time step size limitation of explicit methods in the presence of cut cells is to use unconditionally stable implicit time integration schemes. When using implicit time integration, larger time step sizes can compensate the higher computational costs per time step. The increased computational costs arises from the necessity of solving a system of equations in each time step. One of the most frequently employed implicit methods is the Newmark trapezoidal method [51], which we now briefly discuss as a starting point for the suggested implicit explicit time integration scheme.

The Newmark trapezoidal method is unconditionally stable, non-dissipative, and second-order accurate [55, 17]. It can be implemented as a predictor-corrector scheme and obtained from the general Newmark formulas (equations (<ref>) to (<ref>)) by setting $\beta = \frac{1}{4}$ and $\gamma = \frac{1}{2}$. fig:implicitNewmark illustrates the procedure of the method.

The implicit-explicit Newmark method composed of the implicit trapezoidal Newmark method and the explicit CDM combines the advantages of both methods in the context of the spectral cell method. Firstly, we employ the CDM for the uncut dofs associated with the diagonal mass matrix block. The acceleration of these diagonal dofs can directly be computed using Equation (<ref>), without the need to solve a system of equations. Subsequently, by incorporating this acceleration to the CDM sum step (Equation (<ref>)), we obtain the solution for the uncut dofs at the next time step. Secondly, as the cut dofs correspond to a non-diagonal block in the mass matrix and may cause a very small critical time step size, we integrate them in time using the implicit trapezoidal Newmark method. To predict the displacement and velocity in the trapezoidal Newmark method, we utilize the solution of the previous CDM step for the uncut dofs. In the remainder of this paper, the combination of both algorithms is referred to as the Newmark IMEX method. Its flowchart is depicted in fig:TrapCDMIMEX. The symbols $I_\text{d}$ and $I_\text{c}$ denote the index sets of the diagonal and cut dofs, respectively. Furthermore, note that an update in a dof-subgroup, e.g., $\v{u}^\text{d}$, is automatically assembled into the full vector $\v{u}$ at the corresponding indices $I_\text{d}$.

§.§ Lumping of Cut Cells

We now introduce a lumped scheme to compare it with our implicit-explicit time integration scheme.
A fully diagonal mass matrix that is favorable for explicit time integration is obtained by additionally lumping the cut elements. In the works of [34, 35, 58], Hinton-Rock-Zienkiewicz (HRZ) lumping [59] is employed to approximate the consistent element mass matrices of cut cells using diagonally lumped element mass matrices. In HRZ lumping, the diagonal terms of the cut element mass matrix are scaled to preserve its total mass, while off-diagonal entries are set to zero, i.e.,

\begin{equation}
\tilde{M}_{ii}^e = \frac{m_e}{\sum_{k=1}^{n_e} M_{kk}^e} M_{ii}^e
\end{equation}

\begin{equation}
\tilde{M}_{ij}^e = 0 \qquad \text{if } i \neq j \text{.}
\end{equation}

According to [60], HRZ lumping leads to a reduction of the convergence rates. On the contrary, unlike row-summing (see [17, 45]), it guarantees positive diagonal components and captures the material distribution in the cell more accurately. In the paper at hand, the CDM is used when HRZ lumping is applied for the cut cells. Nicoli et al. [37, 38] lump the cut cells by using moment fitting to adjust the integration weights while retaining the GLL points as integration points. On this lumped SCM system, they integrate in time using a leap-frog algorithm with smaller time step sizes in the cut cells. Another way to mitigate the time step size constraint of SCM, is the use of eigenvalue stabilization, which is investigated in [40]. A comprehensive comparative study of various lumping techniques applied to the SCM is given in [39]. From now on, the SCM with no further lumping on the cut cells is referred to as consistent SCM, and the SCM with additional HRZ lumping on the cut cells is denoted as HRZ-lumped SCM.

§ NUMERICAL EXAMPLES

We now demonstrate fundamental principles of implicit-explicit time integration using a basic system consisting of ten spring-coupled masses before applying the immersed Newmark IMEX method to more complex examples in two and three dimensions.

§.§ Spring-coupled masses

fig:springCoupledMasses shows a system of ten spring-coupled masses that we analyze in this section to discuss key attributes concerning the accuracy and stability of the Newmark IMEX method. The results transfer to systems obtained from spatial discretizations of the wave equation, such as the SCM.

System composed of ten spring-coupled masses

Convergence curves for time integration of the spring-coupled masses with critical time step sizes indicated

Elastic energy for the Newmark IMEX method with a time step size smaller and higher than the critical time step size of explicit subsystem

We assign the first eight masses the value $m_\text{1} = \SI{1}{\kilo\gram}$ and the remaining two the value $m_\text{2} = \SI{e-3}{\kilo\gram}$. The stiffness of all springs is $k = \SI{1}{\newton\per\meter}$. The system is subjected to a harmonic excitation

\begin{equation}
f_\text{t}(t) = \sin \left( 2 \pi f_\text{s} t \right)
\end{equation}

applied to the third last mass with an excitation frequency of $f_\text{s} = \SI{e-1}{\hertz}$. The equations of motion for all masses can be written as a single second-order system

\begin{equation}\label{eq:springCoupledSystem}
\v{M} \ddot{\dv{u}} + \v{K} \dv{u} = \v{f} \text{,}
\end{equation}

where $\dv{u}$ are the dofs of the system, representing the displacements of the point masses. The solution of the corresponding Helmholtz equation yields an amplitude vector whose product with the harmonic excitation $f_\text{t}(t)$ forms our time-dependent reference solution. The equation system (<ref>) is similar to (<ref>), and its structure is identical to a mass-lumped SCM discretization of the wave equation. We emphasize this similarity by referring to the first eight masses as uncut dofs and the remaining two as badly cut dofs. The behavior of the badly cut dofs dominates the stability of the global system.

In the following, we compare the convergence of the explicit CDM, implicit trapezoidal Newmark, and Newmark IMEX time-integration schemes on the time interval $T = \SI{50}{\second}$. The Newmark IMEX method integrates the eight uncut dofs explicitly and the remaining two badly cut dofs implicitly. Our convergence analysis starts with an initial $N_\text{t} = 10$ time steps, resulting in a $\Delta t$ of $\SI{5}{\second}$, and subsequently logarithmically refines the temporal discretization until we reach $N_\text{t} = 10^4$ time steps with a $\Delta t$ of $\SI{5e-3}{\second}$. The initial displacements and velocities are determined from the reference solution outlined above.

fig:springCoupledMassesConvergence compares the convergence of the three methods, using the $l_2$-norm of the difference between the computed and reference displacements as an error measure. The green line indicates second-order convergence, the vertical dashed line on the left highlights the critical time step size $\Delta t_\text{crit}^\text{global} = \SI{3.9086e-2}{\second}$ of the global system, and the dashed line on the right highlights the critical time step size $\Delta t_\text{crit}^\text{explicit} = \SI{1.0154}{\second}$ of the explicit subsystem. The CDM and the trapezoidal Newmark method show the expected second-order convergence, with a slightly lower convergence constant for the CDM. The CDM becomes unstable beyond the global critical time step size $\Delta t_\text{crit}^\text{global}$, resulting in no data points to the right of the corresponding dashed line. In contrast, the trapezoidal Newmark method is unconditionally stable. The Newmark IMEX method preserves the convergence rates of the CDM and trapezoidal Newmark methods, which are its basic building blocks. Here, the convergence constant even reaches that of the CDM. As expected, the Newmark IMEX method remains stable up to the critical time step size $\Delta t_\text{crit}^\text{explicit}$ of the explicit subsystem (corresponding to the uncut dofs), which is much less restrictive than $\Delta t_\text{crit}^\text{global}$.

fig:springCoupledMassesElastic shows the evolution of the elastic energy $E^\text{elastic} = \frac{1}{2} \; \dv{u}^\text{T} \v{K} \dv{u}$ for two Newmark IMEX simulations with a time step size slightly below and slightly above the critical time step size of the explicit subsystem. The elastic energy for $\Delta t = 0.99 \; \Delta t_\text{crit}^\text{explicit}$, plotted in black, is stable and remains bounded. However, as the red curve shows, a time step size slightly larger than $\Delta t_\text{crit}^\text{explicit}$ leads to an unstable simulation with an unboundedly growing elastic energy. We refer to [50, 49, 61] for further discussion of the theory, derivations, and examples of implicit-explicit time integration.

§.§ Perforated plate

$t = \SI{1.5}{\second}$

$t = \SI{3}{\second}$

$t = T = \SI{10}{\second}$

CDM solution of the perforated plate with $\Delta t = \SI{1}{\milli\second}$

Histogram of fill ratios $\eta$ for the mesh shown in fig:RandomHolesForward

We now consider the scalar wave equation to simulate the propagation of acoustic waves in a perforated plate with dimensions $l_x \times l_y = \SI{10}{\meter} \times \SI{4}{\meter}$. We subtract ten circular holes with randomly sampled center points $(x_{c, i}, y_{c, i})$ and radii $r_i$ from uniform distributions:

\begin{equation}
x_{\text{c}, i} \in \mathcal{U}(\SI{2}{\meter}, \SI{10}{\meter}), \quad
y_{\text{c}, i} \in \mathcal{U}(\SI{0}{\meter}, \SI{4}{\meter}), \quad
r_i \in \mathcal{U}(\SI{0.2}{\meter}, \SI{0.6}{\meter}).
\end{equation}

We choose constant values $\rho = \SI{1}{\kg\per\meter\squared}$ and $c = \SI{1}{\meter\per\second}$ for the density and wave speed and scale the density by a factor $\alpha_f = 10^{-6}$ inside the fictitious domain. Starting from the initial condition $u(x, 0) = 0$, we run our simulation until $T = \SI{10}{\second}$, so that waves can travel through the domain exactly once. We assume homogeneous Neumann boundary conditions on all four sides, representing a perfectly reflective behavior. The excitation $f(x, y, t)$ is the product of the spatial function $f_\text{x}(x, y)$ and the temporal function $f_{\text{t}}(t)$. We define $f_\text{t}(t)$ as the derivative of a Gaussian function in time

\begin{equation}
f_\text{t}(t) = \frac{-(t - t_\text{0})}{\sqrt{2 \pi} \sigma_\text{t}^3} e^{\left(\frac{-(t-t_\text{0})^2}{2 \sigma_\text{t}^2}\right)} ,\label{eq:ft}
\end{equation}

centered at $t_\text{0} = \frac{1}{f_\text{s}}$, with a temporal width of $\sigma_\text{t} = \frac{1}{2 \pi f_\text{s}}$. The dominant frequency of the excitation signal is set to $f_\text{s} = \SI{2}{\hertz}$, corresponding to a dominant wavelength of $\lambda = \frac{c}{f} = \SI[parse-numbers=false]{\frac{1}{2}}{\meter}$. The spatial distribution $f_\text{x}$ of the excitation is also a Gaussian bell, centered at the position $(x_s, y_s) = (\SI{1}{\meter}, \SI{2}{\meter})$ and a width of $\sigma_s = \SI{6e-2}{\meter}$:

\begin{equation}
f_\text{x}(x, y) = 10 \cdot e^{\left(-\frac{(x - x_s)^2 + (y - y_s)^2}{2 \sigma_s^2}\right)}\text{.}
\label{eq:fx}
\end{equation}

Cell-wise critical time step sizes

We discretize the domain with a Cartesian SCM grid of $40 \times 16$ cells in the $x$ and $y$ directions and assign a polynomial order of $p = 5$. We integrate cut cells by distributing grids with $p + 1$ Gauss-Legendre points per direction on a quadtree with six levels of refinement towards the domain boundary. We also discard finite cells that are entirely outside the physical domain. fig:RandomHolesForward shows the computational mesh with three solution snapshots of the wave spreading from the initial pulse on the left side. There were so many reflections until the third snapshot that no clear wavefront is visible anymore. The high polynomial degree allows the use of a coarse mesh for the wave field while still satisfying the sampling condition. For $p = 5$, the spatial discretization must resolve approximately any wavelength in the excitation signal with at least one element. Our mesh resolves the dominant wavelength with two cells and even wavelengths of the decaying, higher frequency content up to $2 f_\text{s}$ with at least one cell.

fig:fillRatiosCircles shows a histogram of the fill ratios of all cells. Approximately 80 of the cells are uncut, resulting in $n^\text{d} = 11\,936$ diagonal dofs and $n^\text{c} = 3\,473$ dofs supported by at least one cut cell. We consider cells below the fill ratio threshold $\epsilon = 10^{-10}$ empty. This example's smallest nonzero fill ratio is $\eta_{min} = \num{7.66e-4}$. While only a few cells are in the $[\epsilon, 0.1]$ bucket, they significantly reduce the global critical time step size, as discussed in sec:intro. fig:tCritsCircles shows the cell-wise critical time step sizes for the consistent SCM. The values are uniform across uncut cells but significantly lower for cells intersecting the domain boundary.
tab:criticaltimestepsizes compares the global critical time step, the cell-wise critical time step of uncut cells, and the minimum cell-wise critical time step of cut cells. While the minimum cell-wise value predicts the global critical time step well for the consistent and HRZ-lumped SCM, the Newmark IMEX method only depends on the critical time step of uncut cells. We expect this result since the Newmark IMEX method treats cut cells implicitly, and its critical time step only depends on the explicit, i.e., the uncut dofs. Therefore, the cell-wise critical time step for uncut dofs is a good and cheap estimate for the Newmark IMEX critical time step. Interestingly, the comparison of the consistent and HRZ-lumped SCM approaches shows a larger global critical time step for the consistent SCM despite the more restrictive cell-wise minimum that lies over 16 times below the uncut cell estimate. The eigenvalue stabilization approach presented in [40] may be able to lift the significant restrictions on the critical time step in an explicit treatment of cut cells in the SCM.

Comparison of the critical time step sizes for different SCM versions

global system
$\Delta t_\text{crit}^\text{global}$
cell-wise, uncut
$\displaystyle\Delta t_\text{crit}^\text{uncut}$
cell-wise, minimum
$\displaystyle\Delta t_\text{crit}^\text{cut, min}$
$\displaystyle\frac{\Delta t_\text{crit}^\text{uncut}}{\Delta t_\text{crit}^\text{cut, min}}$

CDM, consistent SCM       2.13     18.13     1.10     16.4
CDM, HRZ-lumped SCM       1.42     —”—                      1.30     13.9
Newmark IMEX SCM          18.14                —”—                      -                            -

Runtime versus time step size for the perforated plate

$L_2$ error versus time step size for the perforated plate

We now investigate the accuracy and runtime of several time integration methods by comparing them to a reference solution computed on a grid of $100 \times 40$ cells with $p = 5$ and integrated in time with $\Delta t = \SI{5e-5}{\second}$. Following [62], we approximate the relative error in the $L_2$-norm summing over randomly sampled coordinates $(x_i, y_i)$ inside the physical domain:

\begin{equation}
e_{L_2} = \sqrt{\frac{\sum_{i=1}^{n_p} (u(x_i, y_i) - u_{ref}(x_i, y_i))^2}{\sum_{i=1}^{n_p} u_{ref}(x_i, y_i)^2}}
\end{equation}

We first distribute random points in the extended, square domain but only keep the $n_\text{p} = 8572$ points inside the physical domain (i.e., those not inside one of the circular holes).

Our previous analysis of the critical time step sizes considered the CDM using the consistent SCM, the CDM using the HRZ-lumped SCM, and our Newmark IMEX approach that also uses the consistent SCM. We now include in our analysis a fully implicit scheme that integrates all dofs of the consistent SCM with an implicit trapezoidal Newmark method. We also consider a leap-frog algorithm since it is often used for similar problems, e.g., in  [37]. It subdivides each time step into $m$ smaller intervals to integrate the cut cells with higher accuracy and improved stability. We refer to [57] for details on the leap-frog formulation and its implementation. We measure the runtime of all simulations on a single Ryzen 7 2700x core. We use the Scipy function  where necessary to factorize the matrices and solve at each time step.

$L_2$-error versus runtime for the perforated plate

In the comparison of the runtime for a fixed number of time steps in fig:dtRuntimeCircles, the CDM using HRZ-lumped SCM and the CDM using the consistent SCM are most efficient, but they can only use small time step sizes due to their explicit treatment of cut cells. The plots exclude the unstable computations exceeding the critical time step size. Due to the unconditional stability, the implicit treatment of all cells in the trapezoidal Newmark approach is possible for the full range of time step sizes. However, the solution of a triangular matrix system with the factorized system matrix $\mathbf{S}$ and a different force vector $\mathbf{\hat{f}}_n$ (see Section <ref>) at each time step makes it significantly more expensive. In contrast, Newmark IMEX is almost as fast as the CDM while being impacted much less by its stability constraint. This result aligns with the analysis of the spring-coupled masses in the previous section. The leap-frog algorithm trades performance to extend its stability range since finer substeps for the cut cells increase the runtime. The two conditions, $\Delta t_\text{coarse} < \Delta t_\text{crit}^\text{uncut}$ and $\Delta t_\text{fine} = \Delta t_\text{coarse} / m < \Delta t_\text{crit}^\text{cut, min}$, determine the stability limit of the leap-frog algorithm. Therefore, $m = 17$ substeps are more expensive but allow a larger coarse time step than $m = 5$.

fig:dtL2Circles compares accuracy depending over the time step size. Except for the CDM applied to HRZ-lumped SCM, all methods converge to the error limit imposed by the chosen spatial discretization. The HRZ lumping of cut cells strongly deteriorates the overall accuracy (see also [39]).
The Newmark IMEX method is the most accurate for given time step sizes, followed by the leap-frog algorithms and the implicit trapezoidal Newmark method.
The two leap-frog algorithms behave similarly since increasing the number of substeps here impacts stability but not so much accuracy. The implicit trapezoidal Newmark method's error is comparatively high, particularly for large time step sizes. fig:RuntimeL2Circles combines the information from fig:dtRuntimeCircles and fig:dtL2Circles and is crucial for determining the most efficient of the investigated algorithms. Newmark IMEX significantly outperforms the other methods, confirming that our approach can leverage the advantages of implicitly integrating cut cells while explicitly integrating uncut cells without suffering from their disadvantages in this example.

§.§ Three-dimensional example

Whole structure

Pillar with source and receiver positions

3D example [63, 13]

Time step sizes and computation times

time step size $\Delta t$     critical time step $\Delta t_\text{crit}$     computation time

CDM, consistent SCM     $\SI{2e-7}{\second}$     $\SI{2.0001e-7}{\second}$     $\SI{510.99}{\second}$

CDM, HRZ-lumped SCM     $\SI{2e-7}{\second}$     $\SI{2.0015e-7}{\second}$     $\SI{326.17}{\second}$

Trapezoidal Newmark SCM     $\SI{1e-6}{\second}$     $\infty$     $\SI{278.4}{\second}$

Newmark IMEX SCM    $\SI{1e-6}{\second}$     $\SI{1.8551e-6}{\second}$     $\SI{116.78}{\second}$

Receiver $R0$ – left: whole signal; right: zoomed-in signal

Receiver $R1$ – left: whole signal; right: zoomed-in signal

Receiver $R2$ – left: whole signal; right: zoomed-in signal

Solutions of the 3D example: wave signals at the three receiver positions

Efficient methods for simulating wave propagation are essential for solving inverse problems in non-destructive testing, where many objective function and gradient evaluations with forward model evaluations are needed. We now analyze the geometry of the three-dimensional full waveform inversion performed in [13]. fig:3Dstructure illustrates the geometry and highlights the left pillar in blue that we extract and simulate in our performance analysis. The geometry is defined by a surface triangulation in the STL format. We define homogeneous Neumann (i.e., perfectly reflecting) boundary conditions on the entire surface. The left pillar of the model is embedded within an extended domain of size $\SI{2}{\meter} \times \SI{1.25}{\meter} \times \SI{1.25}{\meter}$. The material is concrete with a density of $\rho = \SI{2400}{\kg\per\meter\cubed}$ and a wave speed of $c = \SI{3000}{\meter\per\second}$. The density is scaled by $\alpha_f = 10^{-5}$ inside the fictitious domain. We position a source at $S0 = (\SI{1,04066}{\meter}, \SI{0.542926}{\meter}, \SI{1.0}{\meter})$ and three receivers at $R0 = (\SI{0.27966}{\meter}, \SI{0.39788}{\meter}, \SI{0.57757}{\meter})$, $R1 = (\SI{0.34066}{\meter}, \SI{0.18043}{\meter}, \SI{1.0}{\meter})$, and $R2 = (\SI{0.78096}{\meter}, \SI{0.47348}{\meter}, \SI{1.4967}{\meter})$. fig:3Dstructure_sr marks the source with a green cross and the receivers with red circles, where $R0$ is the lowest receiver, $R1$ is the middle receiver, and $R2$ is the highest receiver. The structure is excited at $S0$ with a two-cycle sine burst with a central frequency of $f_\text{s} = \SI{20}{\kilo\hertz}$

\begin{equation}
f_\text{t}(t) = \begin{cases}
\sin{(2 \pi f_\text{s} t)} \sin^2{\left(\frac{\pi f_\text{s} t}{n_\text{c}}\right)} &\qquad \text{, } t \leq \frac{n_\text{c}}{f_\text{s}} \\
0 &\qquad \text{, else}
\end{cases},
\end{equation}

where $n_\text{c} = 2$ denotes the number of cycles. Beyond $2 \times f_\text{s}$, the frequency content of the two-cycle sine burst decreases to almost zero. Thus, $2 \times f_\text{s}$ is considered as the maximum frequency of interest. We use $u, \dot{u} = 0$ as initial condition and simulate until $T = \SI{e-3}{\second}$.

Following sec:theory, we discretize the computational domain using $20 \times 20 \times 32$ cells of polynomial order $p = 4$ in the $x$-, $y$-, and $z$-directions. This mesh resolves the central wavelength with $2.4$ cells and the maximum wavelength with $1.2$ cells. We integrate cut cells on a quadtree with five levels of refinement towards the domain boundary. The SCM discretization of the wave field results in $n^\text{dof} = 325\,573$ basis functions. Given the thin structure of the pillar, only $n^\text{d} = 158\,053$ basis functions are supported exclusively by uncut cells, while $n^\text{c} = 167\,520$ basis functions are supported by at least one cut cell. We compare the performance of the explicit CDM, the trapezoidal Newmark method (treating all dofs implicitly), and our Newmark IMEX method, all three using the consistent SCM and the CDM using HRZ-lumped SCM. Tab:times3D shows the selected time step sizes, the critical time step sizes, and the computational times for the different time integration methods. The runtimes are measured from single-threaded computations on an Intel i9-9900K processor and include only the time integration with an initial matrix factorization, if necessary. The matrix factorization and triangular solution at each time step use  from Intel's oneMKL library. The relatively coarse time step for the Newmark IMEX and trapezoidal Newmark methods leads to $50$ time steps per period of the central frequency and $25$ time steps per period of the maximum frequency. The CDM methods with consistent SCM and HRZ-lumped SCM use a smaller time step size only slightly below their critical time step, resulting in $250$ time steps per period of the central frequency  and $125$ time steps per period of the maximum frequency.

fig:solutions3D compares the time signals at the receiver positions. The solutions of the CDM, the Newmark IMEX method, and the trapezoidal Newmark method are almost identical. The plots on the right side of fig:solutions3D magnify a time interval towards the end of the simulation. We use the CDM as a reference since it uses more time steps and is the most accurate.
While the Newmark IMEX looks identical to the CDM, the trapezoidal Newmark method shows minor deviations in amplitude and phase. The much higher number of time steps makes the CDM the most expensive method. In contrast, Newmark IMEX is the fastest of all methods despite the unfavorable ratio between uncut and cut dofs. The trapezoidal Newmark method requires more than twice the computational time of Newmark IMEX. When applying HRZ lumping to the cut cells, the accuracy of the simulation deteriorates significantly. The waves are similar at first, but the difference increases after several reflections because of the error introduced by the HRZ lumping of cells intersecting the boundary. Despite the fully diagonal mass matrix, the computation time is still almost three times higher than that of Newmark IMEX because it requires significantly more time steps to maintain stability.

§ CONCLUSION

We present an implicit-explicit (IMEX) time integration approach for the immersed scalar wave equation. The spatial discretization by the spectral cell method incorporates complex geometries by an extended domain. In the presented hybrid time integration, called the immersed Newmark IMEX method, we split the system obtained from the spatial discretization into two parts: one corresponding to the basis functions with support on at least one cell intersecting the domain boundary and another consisting of the remaining dofs whose basis function are supported only by uncut cells. The cut dofs are integrated in time with the implicit trapezoidal Newmark method to gain independence on their restrictive critical time step size. The uncut dofs are integrated with explicit second-order central differences. We use nodal lumping to obtain a diagonal mass matrix of the explicit subsystem.

We demonstrate that immersed Newmark IMEX eliminates the severe time step size restrictions imposed by the cut cells by determining the highest eigenfrequency of a cut cell when varying the cut ratios.
The analysis of a simplified system of ten spring-coupled masses shows that Newmark IMEX also preserves second-order convergence. Its stability depends on the explicit subsystem of the uncut dofs. Our two- and three-dimensional examples with significant geometric complexity confirm these results. Immersed Newmark IMEX outperforms other methods significantly when comparing runtime per accuracy since it benefits from the efficiency of the spectral element method  for uncut cells. The direct linear solver for the explicit subsystem costs less than using significantly more time steps. Even the leap-frog algorithm that limits the additional time steps to the explicit subsystem is more expensive per accuracy and time step size. Newmark IMEX performs best even in the unfavorable case of many cut cells, as our three-dimensional example shows. It is four times faster than the explicitly integrated spectral cell method with a consistent mass matrix and more than two times faster with a higher accuracy than the trapezoidal Newmark method.
We also observe that an SCM with an HRZ lumping of cut cells significantly reduces the accuracy while suffering from the same critical time step constraints.

While immersed Newmark IMEX appears viable for a broad range of problems, limitations exist. The implicit integration of cut cells requires solving a system of equations at each time step that we accelerate by computing a factorization in advance.
On the one hand, the memory consumption of direct solvers may limit the number of cut cells for computations on a single machine. The explicit schemes we consider do not store an inverse and can even be implemented matrix-free. On the other hand, larger problems may require a parallel computation on a distributed memory system to reduce the runtime. The involvement of a direct solver significantly complicates the necessary implementations, while methods like the HRZ-lumped spectral cell method are simple to parallelize. Finding ways to upscale the immersed Newmark IMEX method towards massively parallel computations remains a challenge for future investigation.

§ ACKNOWLEDGEMENTS

We gratefully thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for their support through the grants KO 4570/1-1 and RA 624/29-1. We also thank Harald Kloft and Robin Dörrie from the Institut für Tragwerksentwurf of the Technische Universität Braunschweig for providing the geometric model investigated in sec:bench3D which serves as a benchmark in the DFG project 414265976 TRR 277.

[1]
S. Pavlopoulou, C. Soutis, and G. Manson, “Non-destructive inspection of
adhesively bonded patch repairs using lamb waves,” Plastics, Rubber and
Composites, vol. 41, pp. 61–68, 03 2012.

[2]
P. Huthwaite and F. Simonetti, “High-resolution guided wave tomography,” Wave Motion, vol. 50, no. 5, pp. 979–993, 2013.

[3]
J. Rao, M. Ratassepp, and Z. Fan, “Guided wave tomography based on
full-waveform inversion,” IEEE Transactions on Ultrasonics,
Ferroelectrics, and Frequency Control, vol. 63, pp. 1–1, 02 2016.

[4]
A. Zimmermann, P. Huthwaite, and B. Pavlakovic, “High-resolution thickness
maps of corrosion using sh1 guided wave tomography,” Proceedings of the
Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 477,
p. 20200380, 01 2021.

[5]
P. Lailly, “The seismic inverse problem as a sequence of before stack
migrations,” Conference on Inverse Scattering—Theory and Application,
SIAM, pp. 206–220, 1983.

[6]
A. Tarantola, “Inversion of seismic reflection data in the acoustic
approximation,” GEOPHYSICS, vol. 49, no. 8, pp. 1259–1266, 1984.

[7]
A. Fichtner, Full Seismic Waveform Modelling and Inversion.
Advances in Geophysical and Environmental Mechanics and Mathematics,
Springer Berlin Heidelberg, 2011.

[8]
R. G. Pratt, L. Huang, N. Duric, and P. J. Littrup, “Sound-speed and
attenuation imaging of breast tissue using waveform tomography of
transmission ultrasound data,” in SPIE Medical Imaging, 2007.

[9]
G. Y. S. Sandhu, E. West, C. Li, O. Roy, and N. Duric, “3d frequency-domain
ultrasound waveform tomography breast imaging,” in Medical Imaging,

[10]
L. Guasch, O. Calderon Agudo, M.-X. Tang, P. Nachev, and M. Warner,
“Full-waveform inversion imaging of the human brain,” npj Digital
Medicine, vol. 3, 12 2020.

[11]
J. Rao, J. Yang, M. Ratassepp, and Z. Fan, “Multi-parameter reconstruction of
velocity and density using ultrasonic tomography based on full waveform
inversion,” Ultrasonics, vol. 101, p. 106004, 2020.

[12]
T. Bürchner, P. Kopp, S. Kollmannsberger, and E. Rank, “Immersed boundary
parametrizations for full waveform inversion,” Computer Methods in
Applied Mechanics and Engineering, vol. 406, p. 115893, 2023.

[13]
T. Bürchner, P. Kopp, S. Kollmannsberger, and E. Rank, “Isogeometric
multi-resolution full waveform inversion based on the finite cell method,”
Computer Methods in Applied Mechanics and Engineering, p. 116286, 2023.

[14]
K. R. Kelly, R. W. Ward, S. Treitel, and R. M. Alford, “Synthetic
seismograms; a finite-difference approach,” Geophysics, vol. 41,
pp. 2–27, 02 1976.

[15]
J. Virieux, “P-sv wave propagation in heterogeneous media: Velocity‐stress
finite‐difference method,” GEOPHYSICS, vol. 51, no. 4, pp. 889–901,

[16]
H. Bao, J. Bielak, O. Ghattas, L. F. Kallivokas, D. R. O'Hallaron, J. R.
Shewchuk, and J. Xu, “Large-scale simulation of elastic wave propagation in
heterogeneous media on parallel computers,” Computer Methods in Applied
Mechanics and Engineering, vol. 152, no. 1, pp. 85–102, 1998.
Containing papers presented at the Symposium on Advances in
Computational Mechanics.

[17]
T. J. Hughes, The finite element method: linear static and dynamic finite
element analysis.
Courier Corporation, 2012.

[18]
A. Düster, E. Rank, and B. Szabó, “The p-Version of the Finite Element and
Finite Cell Methods,” in Encyclopedia of Computational Mechanics
Second Edition (E. Stein, R. Borst, and T. Hughes, eds.), pp. 1–35, 2017.

[19]
J. A. Cottrell, A. Reali, Y. Bazilevs, and T. J. R. Hughes, “Isogeometric
analysis of structural vibrations,” Computer Methods in Applied
Mechanics and Engineering, vol. 195, no. 41, pp. 5257–5296, 2006.

[20]
T. J. Hughes, J. A. Evans, and A. Reali, “Finite element and nurbs
approximations of eigenvalue, boundary-value, and initial-value problems,”
Computer Methods in Applied Mechanics and Engineering, vol. 272,
pp. 290–320, 2014.

[21]
L. Demkowicz, Computing with hp-adaptive Finite Elements: Volume 1 One and
Two Dimensional Elliptic and Maxwell Problems.
Chapman & Hall/CRC Applied Mathematics & Nonlinear Science, CRC
Press, 2006.

[22]
T. J. R. Hughes, The finite element method: linear static and dynamic
finite element analysis.
Dover Publications, 2000.

[23]
N. M. Newmark, “A method of computation for structural dynamics,” Journal of the Engineering Mechanics Divison, vol. 85, 1959.

[24]
J. C. Butcher, “Implicit runge-kutta processes,” Mathematics of
Computation, vol. 18, pp. 50–64, 1964.

[25]
J. Dormand and P. Prince, “A family of embedded runge-kutta formulae,” Journal of Computational and Applied Mathematics, vol. 6, no. 1, pp. 19–26,

[26]
D. Komatitsch and J. Tromp, “Spectral-element simulations of global seismic
wave propagation—I. Validation,” Geophysical Journal International,
vol. 149, pp. 390–412, 05 2002.

[27]
G. C. Cohen, Higher-Order Numerical Methods for Transient Wave Equation.
Springer-Verlag: Berlin-Heidelberg, 2002.

[28]
T. J. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: Cad,
finite elements, nurbs, exact geometry and mesh refinement,” Computer
Methods in Applied Mechanics and Engineering, vol. 194, pp. 4135–4195, 10

[29]
J. Cotrell, T. Hughes, and Y. Bazilevs, Isogeometric Analysis: Toward
Integration of CAD and FEA.
Wiley, 2009.

[30]
T.-H. Nguyen, R. R. Hiemstra, S. Eisenträger, and D. Schillinger, “Towards
higher-order accurate mass lumping in explicit isogeometric analysis for
structural dynamics,” Computer Methods in Applied Mechanics and
Engineering, p. 116233, 2023.

[31]
S. Held, S. Eisenträger, and W. Dornisch, “An efficient mass lumping scheme
for isogeometric analysis based on approximate dual basis functions.”
http://arxiv.org/abs/2306.12257, 2023.

[32]
A. Düster, J. Parvizian, Z. Yang, and E. Rank, “The finite cell method for
three-dimensional problems of solid mechanics,” Computer Methods in
Applied Mechanics and Engineering, vol. 197, no. 45, pp. 3768–3782, 2008.

[33]
M. Joulaian, S. Hubrich, and A. Düster, “Numerical integration of
discontinuities on arbitrary domains based on moment fitting,” Comput.
Mech., vol. 57, p. 979–999, jun 2016.

[34]
S. Duczek, M. Joulaian, A. Düster, and U. Gabbert, “Numerical analysis of
lamb waves using the finite and spectral cell methods,” International
Journal for Numerical Methods in Engineering, vol. 99, no. 1, pp. 26–53,

[35]
M. Joulaian, S. Duczek, U. Gabbert, and A. Düster, “Finite and spectral cell
method for wave propagation in heterogeneous materials,” Computational
Mechanics, vol. 54, 2014.

[36]
S. Duczek, S. Liefold, and U. Gabbert, “The finite and spectral cell methods
for smart structure applications: Transient analysis,” Acta Mechanica,
vol. 226, no. 3, pp. 845–869, 2015.

[37]
S. Nicoli, K. Agathos, and E. Chatzi, “Moment fitted cut spectral elements for
explicit analysis of guided wave propagation,” Computer Methods in
Applied Mechanics and Engineering, vol. 398, 8 2022.

[38]
S. Nicoli, K. Agathos, P. Kudela, and E. Chatzi, “A moment-fitted extended
spectral cell method for structural health monitoring applications,” Applied Sciences, vol. 13, 09 2023.

[39]
M. Kelemen, A Review of Mass Lumping Schemes for the Spectral Cell
Master's thesis, Technical University of Munich, 2021.

[40]
S. Eisenträger, L. Radtke, W. Garhuom, S. Löhnert, A. Düster, D. Juhre, and
D. Schillinger, “An eigenvalue stabilization technique for immersed boundary
finite element methods in explicit dynamics,” arXiv:2310.11935v1,

[41]
D. Schillinger, L. Dedè, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, and
T. J. Hughes, “An isogeometric design-through-analysis methodology based on
adaptive hierarchical refinement of nurbs, immersed boundary methods, and
t-spline cad surfaces,” Computer Methods in Applied Mechanics and
Engineering, vol. 249-252, pp. 116–150, 2012.
Higher Order Finite Element and Isogeometric Methods.

[42]
E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, and A. Düster,
“Geometric modeling, isogeometric analysis and the finite cell method,”
Computer Methods in Applied Mechanics and Engineering, vol. 249-252,
pp. 104–115, 2012.
Higher Order Finite Element and Isogeometric Methods.

[43]
L. Leidinger, M. Breitenberger, A. Bauer, S. Hartmann, R. Wüchner, K.-U.
Bletzinger, F. Duddeck, and L. Song, “Explicit dynamic isogeometric b-rep
analysis of penalty-coupled trimmed nurbs shells,” Computer Methods in
Applied Mechanics and Engineering, vol. 351, pp. 891–927, 2019.

[44]
M. Meßmer, T. Teschemacher, L. F. Leidinger, R. Wüchner, and K.-U.
Bletzinger, “Efficient cad-integrated isogeometric analysis of trimmed
solids,” Computer Methods in Applied Mechanics and Engineering,
vol. 400, p. 115584, 2022.

[45]
Y. Voet, E. Sande, and A. Buffa, “A mathematical theory for mass lumping and
its generalization with applications to isogeometric analysis,” Computer Methods in Applied Mechanics and Engineering, vol. 410, p. 116033,

[46]
S. K. Stoter, S. C. Divi, E. H. van Brummelen, M. G. Larson, F. de Prenter,
and C. V. Verhoosel, “Critical time-step size analysis and mass scaling by
ghost-penalty for immersogeometric explicit dynamics,” Computer Methods
in Applied Mechanics and Engineering, vol. 412, p. 116074, 2023.

[47]
T. Belytschko, H.-J. Yen, and R. Mullen, “Mixed methods for time
integration,” Computer Methods in Applied Mechanics and Engineering,
vol. 17-18, pp. 259–275, 1979.

[48]
T. Belytschko and R. Mullen, “Stability of explicit-implicit mesh partitions
in time integration,” INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN
ENGINEERING, vol. 12, pp. 1575–1586, 1978.

[49]
T. J. R. Hughes and W. K. Liu, “Implicit-explicit finite elements in transient
analysis: Implementation and numerical examples,” Journal of Applied
Mechanics, Transactions ASME, vol. 45, p. 375 – 378, 1978.
Cited by: 176.

[50]
T. J. R. Hughes and W. K. Liu, “Implicit-explicit finite elements in transient
analysis: Stability theory.,” Journal of Applied Mechanics,
Transactions ASME, vol. 45, p. 371 – 374, 1978.
Cited by: 280.

[51]
K. J. Bathe, “Conserving energy and momentum in nonlinear dynamics: A simple
implicit time integration scheme,” Computers & Structures, vol. 85,
pp. 437–445, 4 2007.

[52]
D. Rabinovich, D. Givoli, and E. Turkel, “Single‐field identification of
inclusions and cavities in an elastic medium,” International Journal
for Numerical Methods in Engineering, 09 2023.

[53]
M. Petö, F. Duvigneau, and S. Eisenträger, “Enhanced numerical integration
scheme based on image-compression techniques: application to fictitious
domain methods,” Advanced Modeling and Simulation in Engineering
Sciences, vol. 7, 2020.

[54]
L. Kudela, N. Zander, S. Kollmannsberger, and E. Rank, “Smart octrees:
Accurately integrating discontinuous functions in 3d,” Computer Methods
in Applied Mechanics and Engineering, vol. 306, 04 2016.

[55]
M. Gérardin and D. Rixen, Mechanical Vibrations: Theory and Application
to Structural Dynamics.

[56]
W. K. Liu and T. Belytschko, “Mixed-time implicit-explicit finite elements for
transient analysis,” Computers & Structures, vol. 15, pp. 445–450, 1

[57]
C. Faßbender, Time Integration for the Spectral Cell Methods with
Application to the Full Waveform Inversion.
Master's thesis, Technical University of Munich, 2023.

[58]
F. Mossaiby, M. Joulaian, and A. Düster, “The spectral cell method for wave
propagation in heterogeneous materials simulated on multiple gpus and cpus,”
Computational Mechanics, vol. 63, pp. 805–819, 05 2019.

[59]
E. Hinton, T. Rock, and O. C. Zienkiewicz, “A note on mass lumping and related
processes in the finite element method,” Earthquake Engineering &
Structural Dynamics, vol. 4, pp. 245–249, 1976.

[60]
S. Duczek and H. Gravenkamp, “Critical assessment of different mass lumping
schemes for higher order serendipity finite elements,” Computer Methods
in Applied Mechanics and Engineering, vol. 350, pp. 836–897, 6 2019.

[61]
T. J. R. Hughes and R. A. Stephenson, “Convergence of implicit-explicit
algorithms in nonlinear transient analysis,” International Journal of
Engineering Science, vol. 19, p. 295 – 302, 1981.
Cited by: 11.

[62]
S. Duczek and H. Gravenkamp, “Mass lumping techniques in the spectral element
method: On the equivalence of the row-sum, nodal quadrature, and diagonal
scaling methods,” Computer Methods in Applied Mechanics and
Engineering, vol. 353, pp. 516–569, 8 2019.

[63]
H. Kloft and R. Dörrie, “Shelltonic benchmark.” Additive Manufacturing in
Construction, https://amc-trr277.de, 2023.

[1]
S. Pavlopoulou, C. Soutis, and G. Manson, “Non-destructive inspection of
adhesively bonded patch repairs using lamb waves,” Plastics, Rubber and
Composites, vol. 41, pp. 61–68, 03 2012.

[2]
P. Huthwaite and F. Simonetti, “High-resolution guided wave tomography,” Wave Motion, vol. 50, no. 5, pp. 979–993, 2013.

[3]
J. Rao, M. Ratassepp, and Z. Fan, “Guided wave tomography based on
full-waveform inversion,” IEEE Transactions on Ultrasonics,
Ferroelectrics, and Frequency Control, vol. 63, pp. 1–1, 02 2016.

[4]
A. Zimmermann, P. Huthwaite, and B. Pavlakovic, “High-resolution thickness
maps of corrosion using sh1 guided wave tomography,” Proceedings of the
Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 477,
p. 20200380, 01 2021.

[5]
P. Lailly, “The seismic inverse problem as a sequence of before stack
migrations,” Conference on Inverse Scattering—Theory and Application,
SIAM, pp. 206–220, 1983.

[6]
A. Tarantola, “Inversion of seismic reflection data in the acoustic
approximation,” GEOPHYSICS, vol. 49, no. 8, pp. 1259–1266, 1984.

[7]
A. Fichtner, Full Seismic Waveform Modelling and Inversion.
Advances in Geophysical and Environmental Mechanics and Mathematics,
Springer Berlin Heidelberg, 2011.

[8]
R. G. Pratt, L. Huang, N. Duric, and P. J. Littrup, “Sound-speed and
attenuation imaging of breast tissue using waveform tomography of
transmission ultrasound data,” in SPIE Medical Imaging, 2007.

[9]
G. Y. S. Sandhu, E. West, C. Li, O. Roy, and N. Duric, “3d frequency-domain
ultrasound waveform tomography breast imaging,” in Medical Imaging,

[10]
L. Guasch, O. Calderon Agudo, M.-X. Tang, P. Nachev, and M. Warner,
“Full-waveform inversion imaging of the human brain,” npj Digital
Medicine, vol. 3, 12 2020.

[11]
J. Rao, J. Yang, M. Ratassepp, and Z. Fan, “Multi-parameter reconstruction of
velocity and density using ultrasonic tomography based on full waveform
inversion,” Ultrasonics, vol. 101, p. 106004, 2020.

[12]
T. Bürchner, P. Kopp, S. Kollmannsberger, and E. Rank, “Immersed boundary
parametrizations for full waveform inversion,” Computer Methods in
Applied Mechanics and Engineering, vol. 406, p. 115893, 2023.

[13]
T. Bürchner, P. Kopp, S. Kollmannsberger, and E. Rank, “Isogeometric
multi-resolution full waveform inversion based on the finite cell method,”
Computer Methods in Applied Mechanics and Engineering, p. 116286, 2023.

[14]
K. R. Kelly, R. W. Ward, S. Treitel, and R. M. Alford, “Synthetic
seismograms; a finite-difference approach,” Geophysics, vol. 41,
pp. 2–27, 02 1976.

[15]
J. Virieux, “P-sv wave propagation in heterogeneous media: Velocity‐stress
finite‐difference method,” GEOPHYSICS, vol. 51, no. 4, pp. 889–901,

[16]
H. Bao, J. Bielak, O. Ghattas, L. F. Kallivokas, D. R. O'Hallaron, J. R.
Shewchuk, and J. Xu, “Large-scale simulation of elastic wave propagation in
heterogeneous media on parallel computers,” Computer Methods in Applied
Mechanics and Engineering, vol. 152, no. 1, pp. 85–102, 1998.
Containing papers presented at the Symposium on Advances in
Computational Mechanics.

[17]
T. J. Hughes, The finite element method: linear static and dynamic finite
element analysis.
Courier Corporation, 2012.

[18]
A. Düster, E. Rank, and B. Szabó, “The p-Version of the Finite Element and
Finite Cell Methods,” in Encyclopedia of Computational Mechanics
Second Edition (E. Stein, R. Borst, and T. Hughes, eds.), pp. 1–35, 2017.

[19]
J. A. Cottrell, A. Reali, Y. Bazilevs, and T. J. R. Hughes, “Isogeometric
analysis of structural vibrations,” Computer Methods in Applied
Mechanics and Engineering, vol. 195, no. 41, pp. 5257–5296, 2006.

[20]
T. J. Hughes, J. A. Evans, and A. Reali, “Finite element and nurbs
approximations of eigenvalue, boundary-value, and initial-value problems,”
Computer Methods in Applied Mechanics and Engineering, vol. 272,
pp. 290–320, 2014.

[21]
L. Demkowicz, Computing with hp-adaptive Finite Elements: Volume 1 One and
Two Dimensional Elliptic and Maxwell Problems.
Chapman & Hall/CRC Applied Mathematics & Nonlinear Science, CRC
Press, 2006.

[22]
T. J. R. Hughes, The finite element method: linear static and dynamic
finite element analysis.
Dover Publications, 2000.

[23]
N. M. Newmark, “A method of computation for structural dynamics,” Journal of the Engineering Mechanics Divison, vol. 85, 1959.

[24]
J. C. Butcher, “Implicit runge-kutta processes,” Mathematics of
Computation, vol. 18, pp. 50–64, 1964.

[25]
J. Dormand and P. Prince, “A family of embedded runge-kutta formulae,” Journal of Computational and Applied Mathematics, vol. 6, no. 1, pp. 19–26,

[26]
D. Komatitsch and J. Tromp, “Spectral-element simulations of global seismic
wave propagation—I. Validation,” Geophysical Journal International,
vol. 149, pp. 390–412, 05 2002.

[27]
G. C. Cohen, Higher-Order Numerical Methods for Transient Wave Equation.
Springer-Verlag: Berlin-Heidelberg, 2002.

[28]
T. J. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: Cad,
finite elements, nurbs, exact geometry and mesh refinement,” Computer
Methods in Applied Mechanics and Engineering, vol. 194, pp. 4135–4195, 10

[29]
J. Cotrell, T. Hughes, and Y. Bazilevs, Isogeometric Analysis: Toward
Integration of CAD and FEA.
Wiley, 2009.

[30]
T.-H. Nguyen, R. R. Hiemstra, S. Eisenträger, and D. Schillinger, “Towards
higher-order accurate mass lumping in explicit isogeometric analysis for
structural dynamics,” Computer Methods in Applied Mechanics and
Engineering, p. 116233, 2023.

[31]
S. Held, S. Eisenträger, and W. Dornisch, “An efficient mass lumping scheme
for isogeometric analysis based on approximate dual basis functions.”
http://arxiv.org/abs/2306.12257, 2023.

[32]
A. Düster, J. Parvizian, Z. Yang, and E. Rank, “The finite cell method for
three-dimensional problems of solid mechanics,” Computer Methods in
Applied Mechanics and Engineering, vol. 197, no. 45, pp. 3768–3782, 2008.

[33]
M. Joulaian, S. Hubrich, and A. Düster, “Numerical integration of
discontinuities on arbitrary domains based on moment fitting,” Comput.
Mech., vol. 57, p. 979–999, jun 2016.

[34]
S. Duczek, M. Joulaian, A. Düster, and U. Gabbert, “Numerical analysis of
lamb waves using the finite and spectral cell methods,” International
Journal for Numerical Methods in Engineering, vol. 99, no. 1, pp. 26–53,

[35]
M. Joulaian, S. Duczek, U. Gabbert, and A. Düster, “Finite and spectral cell
method for wave propagation in heterogeneous materials,” Computational
Mechanics, vol. 54, 2014.

[36]
S. Duczek, S. Liefold, and U. Gabbert, “The finite and spectral cell methods
for smart structure applications: Transient analysis,” Acta Mechanica,
vol. 226, no. 3, pp. 845–869, 2015.

[37]
S. Nicoli, K. Agathos, and E. Chatzi, “Moment fitted cut spectral elements for
explicit analysis of guided wave propagation,” Computer Methods in
Applied Mechanics and Engineering, vol. 398, 8 2022.

[38]
S. Nicoli, K. Agathos, P. Kudela, and E. Chatzi, “A moment-fitted extended
spectral cell method for structural health monitoring applications,” Applied Sciences, vol. 13, 09 2023.

[39]
M. Kelemen, A Review of Mass Lumping Schemes for the Spectral Cell
Master's thesis, Technical University of Munich, 2021.

[40]
S. Eisenträger, L. Radtke, W. Garhuom, S. Löhnert, A. Düster, D. Juhre, and
D. Schillinger, “An eigenvalue stabilization technique for immersed boundary
finite element methods in explicit dynamics,” arXiv:2310.11935v1,

[41]
D. Schillinger, L. Dedè, M. A. Scott, J. A. Evans, M. J. Borden, E. Rank, and
T. J. Hughes, “An isogeometric design-through-analysis methodology based on
adaptive hierarchical refinement of nurbs, immersed boundary methods, and
t-spline cad surfaces,” Computer Methods in Applied Mechanics and
Engineering, vol. 249-252, pp. 116–150, 2012.
Higher Order Finite Element and Isogeometric Methods.

[42]
E. Rank, M. Ruess, S. Kollmannsberger, D. Schillinger, and A. Düster,
“Geometric modeling, isogeometric analysis and the finite cell method,”
Computer Methods in Applied Mechanics and Engineering, vol. 249-252,
pp. 104–115, 2012.
Higher Order Finite Element and Isogeometric Methods.

[43]
L. Leidinger, M. Breitenberger, A. Bauer, S. Hartmann, R. Wüchner, K.-U.
Bletzinger, F. Duddeck, and L. Song, “Explicit dynamic isogeometric b-rep
analysis of penalty-coupled trimmed nurbs shells,” Computer Methods in
Applied Mechanics and Engineering, vol. 351, pp. 891–927, 2019.

[44]
M. Meßmer, T. Teschemacher, L. F. Leidinger, R. Wüchner, and K.-U.
Bletzinger, “Efficient cad-integrated isogeometric analysis of trimmed
solids,” Computer Methods in Applied Mechanics and Engineering,
vol. 400, p. 115584, 2022.

[45]
Y. Voet, E. Sande, and A. Buffa, “A mathematical theory for mass lumping and
its generalization with applications to isogeometric analysis,” Computer Methods in Applied Mechanics and Engineering, vol. 410, p. 116033,

[46]
S. K. Stoter, S. C. Divi, E. H. van Brummelen, M. G. Larson, F. de Prenter,
and C. V. Verhoosel, “Critical time-step size analysis and mass scaling by
ghost-penalty for immersogeometric explicit dynamics,” Computer Methods
in Applied Mechanics and Engineering, vol. 412, p. 116074, 2023.

[47]
T. Belytschko, H.-J. Yen, and R. Mullen, “Mixed methods for time
integration,” Computer Methods in Applied Mechanics and Engineering,
vol. 17-18, pp. 259–275, 1979.

[48]
T. Belytschko and R. Mullen, “Stability of explicit-implicit mesh partitions
in time integration,” INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN
ENGINEERING, vol. 12, pp. 1575–1586, 1978.

[49]
T. J. R. Hughes and W. K. Liu, “Implicit-explicit finite elements in transient
analysis: Implementation and numerical examples,” Journal of Applied
Mechanics, Transactions ASME, vol. 45, p. 375 – 378, 1978.
Cited by: 176.

[50]
T. J. R. Hughes and W. K. Liu, “Implicit-explicit finite elements in transient
analysis: Stability theory.,” Journal of Applied Mechanics,
Transactions ASME, vol. 45, p. 371 – 374, 1978.
Cited by: 280.

[51]
K. J. Bathe, “Conserving energy and momentum in nonlinear dynamics: A simple
implicit time integration scheme,” Computers & Structures, vol. 85,
pp. 437–445, 4 2007.

[52]
D. Rabinovich, D. Givoli, and E. Turkel, “Single‐field identification of
inclusions and cavities in an elastic medium,” International Journal
for Numerical Methods in Engineering, 09 2023.

[53]
M. Petö, F. Duvigneau, and S. Eisenträger, “Enhanced numerical integration
scheme based on image-compression techniques: application to fictitious
domain methods,” Advanced Modeling and Simulation in Engineering
Sciences, vol. 7, 2020.

[54]
L. Kudela, N. Zander, S. Kollmannsberger, and E. Rank, “Smart octrees:
Accurately integrating discontinuous functions in 3d,” Computer Methods
in Applied Mechanics and Engineering, vol. 306, 04 2016.

[55]
M. Gérardin and D. Rixen, Mechanical Vibrations: Theory and Application
to Structural Dynamics.

[56]
W. K. Liu and T. Belytschko, “Mixed-time implicit-explicit finite elements for
transient analysis,” Computers & Structures, vol. 15, pp. 445–450, 1

[57]
C. Faßbender, Time Integration for the Spectral Cell Methods with
Application to the Full Waveform Inversion.
Master's thesis, Technical University of Munich, 2023.

[58]
F. Mossaiby, M. Joulaian, and A. Düster, “The spectral cell method for wave
propagation in heterogeneous materials simulated on multiple gpus and cpus,”
Computational Mechanics, vol. 63, pp. 805–819, 05 2019.

[59]
E. Hinton, T. Rock, and O. C. Zienkiewicz, “A note on mass lumping and related
processes in the finite element method,” Earthquake Engineering &
Structural Dynamics, vol. 4, pp. 245–249, 1976.

[60]
S. Duczek and H. Gravenkamp, “Critical assessment of different mass lumping
schemes for higher order serendipity finite elements,” Computer Methods
in Applied Mechanics and Engineering, vol. 350, pp. 836–897, 6 2019.

[61]
T. J. R. Hughes and R. A. Stephenson, “Convergence of implicit-explicit
algorithms in nonlinear transient analysis,” International Journal of
Engineering Science, vol. 19, p. 295 – 302, 1981.
Cited by: 11.

[62]
S. Duczek and H. Gravenkamp, “Mass lumping techniques in the spectral element
method: On the equivalence of the row-sum, nodal quadrature, and diagonal
scaling methods,” Computer Methods in Applied Mechanics and
Engineering, vol. 353, pp. 516–569, 8 2019.

[63]
H. Kloft and R. Dörrie, “Shelltonic benchmark.” Additive Manufacturing in
Construction, https://amc-trr277.de, 2023.