Dispersion Properties of Explicit Finite Element Methods for Wave Propagation Modelling on Tetrahedral Meshes

We analyse the dispersion properties of two types of explicit finite element methods for modelling acoustic and elastic wave propagation on tetrahedral meshes, namely mass-lumped finite element methods and symmetric interior penalty discontinuous Galerkin methods, both combined with a suitable Lax--Wendroff time integration scheme. The dispersion properties are obtained semi-analytically using standard Fourier analysis. Based on the dispersion analysis, we give an indication of which method is the most efficient for a given accuracy, how many elements per wavelength are required for a given accuracy, and how sensitive the accuracy of the method is to poorly shaped elements.


Introduction
Realistic wave propagation problems often involve large three-dimensional domains consisting of heterogeneous materials with complex geometries and sharp interfaces. Solving such problems requires a numerical method that is efficient in terms of computation time and is flexible enough to capture the effect of a complex geometry.
Standard finite difference methods fall short, since they rely on Cartesian grids that cannot efficiently capture the effect of complex interfaces and boundary layers. Finite element methods can overcome this problem when the elements are aligned with those surfaces. However, the accuracy of the finite element method quickly deteriorates when the elements are poorly shaped or are poorly aligned with the geometry. Obtaining a high quality mesh is therefore quintessential. While both hexahedral and tetrahedral elements are commonly used for three-dimensional problems, tetrahedral elements have a big advantage in this respect, since they offer more geometric flexibility and since robust tetrahedral mesh generators based on the Delaunay criterion are available [29,35].
Apart from the construction of a high-quality mesh, finite element methods for wave propagation problems also require a (block)-diagonal mass matrix to enable explicit time-stepping. A diagonal mass matrix can be obtained with nodal basis functions and a quadrature rule, if the quadrature points coincide with the basis function nodes. This technique is known as mass-lumping. For quadrilaterals and hexahedra, masslumping is achieved by combining tensor-product basis functions with a Gauss-Lobatto quadrature rule, resulting in a scheme known as the spectral element method [30,33,20]. For triangles and tetrahedra, an efficient linear mass-lumped scheme is obtained by combining standard Lagrangian basis functions with a Newton-Cotes quadrature rule. For higher-degree triangles and tetrahedra, however, this approach results in an unstable, unsolvable, or inaccurate scheme. To remain accurate and stable, the space of the triangle or tetrahedron is enriched with higher-degree bubble functions. This approach has led to accurate mass-lumped triangles of degree 2 [15], 3 [7], 4 [25], 5 [5], 6 [27], 7-9 [24,9] and tetrahedra of degree 2 [25] and 3 [5].
Another way to obtain a (block)-diagonal mass matrix is by using discontinuous basis functions. The resulting schemes are known as Discontinuous Galerkin (DG) methods. The first DG methods for wave propagation problems were based on a first-order formulation of the wave equation [31,6]. In [32] and [17], DG methods were introduced that were based on the original second-order formulation of the wave problem. The advantage of finite element methods based on the second-order formulation is that they do not need to compute or store the auxiliary variables that appear in the first-order formulation. Moreover, they can be combined with a leap-frog or higher-order Lax-Wendroff time integration scheme that only requires K stages for a 2K-order accuracy. We focus on the symmetric interior penalty discontinuous Galerkin (SIPDG) method, presented and analysed in [17], which is based on the second-order formulation of the wave problem and which also remains energy-conservative on the discrete level. To remain accurate and stable, face integrals and interior penalty parameters are added to the discrete operator. We consider two choices for the penalty parameter: the penalty term derived in [28], based on the trace inequality of [36], and a recently developed sharper estimate [16], based on a more involved trace inequality.
To effectively apply these methods, it is crucial to know the required mesh resolution for a given accuracy. It is also useful to know which method is the most efficient for a given accuracy and how the mesh quality and material parameters, such as the P/Swave velocity ratio for elastic waves, affect the accuracy. A practical and common measure for the accuracy of these type of methods is the amount of numerical dispersion and dissipation. In this paper, we will focus mainly on the numerical dispersion, since the methods we consider are all energy-conservative and therefore do not suffer from numerical dissipation. We do, however, also investigate the spurious modes that appear when projecting a physical wave onto the discrete space.
The dispersion properties of DG methods based on the first-order formulation of the wave problem have already been analysed for Cartesian meshes [18,1], triangles [18,22], and tetrahedra [19]. For the SIPDG method, these properties have already been analysed for Cartesian meshes in [2,13] and triangles in [3] and for the mass-lumped finite element method this has already been analysed for quadrilaterals and hexahedra in [26,8,11] and for triangles in [23]. However, a dispersion analysis of the mass-lumped finite element and SIPDG methods for tetrahedra is, to the best of our knowledge, still missing, even though most realistic wave problems involve three-dimensional domains for which tetrahedral elements are particularly suitable. In this paper, we therefore present an extensive dispersion analysis of these methods for tetrahedra. This analysis is based on standard Fourier analysis. We use the analysis to obtain estimates for the required number of elements per wavelength and estimate the computational cost to obtain an indication of which method is the most efficient for a given accuracy. We consider both acoustic and elastic waves and also look at the effect of poorly shaped elements and high P/S-wave velocity ratios on the accuracy of the methods. This paper is organised as follows: in Section 2, we introduce the tensor notation used in this paper. The acoustic and elastic wave equations are presented in Section 3 and the mass-lumped and discontinuous Galerkin finite element methods are presented in Section 4. In Section 5, we explain how we analyse the dispersion properties of these methods. The results of this analysis are presented in Section 6 and the main conclusions are summarised in Section 7.

Some Tensor Notation
Before we present the acoustic and elastic wave equations, we explain the tensor notation that we use throughout this paper. We let the dot product of two tensors denote the summation over the last index of the left and first index of the right tensor. For the double dot product we also sum over the last-but-one index of the left and second index of the right tensor. A concatenation of two tensors denotes the standard tensor product. To give some examples, letn ∈ R d , u ∈ R m be two vectors, σ ∈ R d×m a second-order tensor, and C ∈ R d×m×m×d a fourth-order tensor. Then for all i = 1, . . . , d and j, q = 1, . . . , m.
In the next section we will use this tensor notation to present the acoustic and isotropic elastic wave equations.

The Acoustic and Isotropic Elastic Wave Equations
Let Ω ⊂ R 3 be a three-dimensional open domain with a Lipschitz boundary ∂Ω, and let (0, T ) be the time domain. Also, let {Γ d , Γ n } be a partition of ∂Ω, corresponding to Dirichlet and von Neumann boundary conditions, respectively. We define the following linear hyperbolic problem: where u : Ω × (0, T ) → R m is a vector of m variables that are to be solved, ∇ is the gradient operator, ρ : Ω → R + is a positive scalar field, C : Ω → R 3×m×m×3 a fourthorder tensor field, f : Ω × (0, T ) → R m the source field, andn : ∂Ω → R 3 the outward pointing normal unit vector.
By choosing the appropriate tensor and scalar field we can obtain the acoustic wave equation and the isotropic elastic wave equations.
These equations can be solved with the finite element methods described in the next section.

The Discontinuous Galerkin and Mass-Lumped Finite Element Method
4.1. The Classical Finite Element Method. Let T h be a tetrahedral tessellation of Ω, with h denoting the radius of the smallest sphere that can contain each element and let U h be the finite element space consisting of continuous element-wise polynomial basis functions satisfying boundary condition (1b). The classical conforming finite element formulation of (1) is finding where (·, ·) denotes the standard L 2 inner product, Π h : L 2 (Ω) → U h denotes the weighted L 2 -projection operator defined such that (ρΠ h u, w) = (ρu, w) for all w ∈ U h , and a : H 1 (Ω) m × H 1 (Ω) m → R is the (semi)-elliptic operator given by a(u, w) := Ω (∇u) t : C : ∇w dx.
Let {w (i) } n i=1 be the set of basis functions spanning U h , and let, for any u ∈ L 2 (Ω) m , the vector u ∈ R n be defined such that n i=1 u i w (i) = Π h u. Also, let M, A ∈ R n×n be the mass matrix and stiffness matrix, respectively, defined by M ij := (ρw (i) , w (j) ) and A ij := a(w (i) , w (j) ), and let f * : [0, T ] → R n be given by f * i := (f , w (i) ). The finite element method can then be formulated as finding u h : [0, T ] → R n such that The main drawback of the classical conforming finite element approach is that when an explicit time integration scheme is applied, a system of equations of the form M x = b needs to be solved at every time step, with M not (block)-diagonal. For large-scale problems, this results in a very inefficient time stepping scheme. This problem can be circumvented by lumping the mass matrix into a diagonal matrix or by using discontinuous basis functions.

4.2.
Mass-Lumping. When using nodal basis functions, the mass matrix can be lumped into a diagonal matrix by taking the sum over each row. This is equivalent to replacing the inner product (·, ·) by (·, ·) (L) h , in which the element integrals are approximated by a quadrature rule with quadrature points that coincide with the nodes of the basis functions. We can write (u, w) where Q e denotes the quadrature points on element e and ω e,x denote the quadrature weights. Let {x (i) } n i=1 denote the global set of integration points and define w (i) to be the nodal basis function corresponding to x (i) , so w (i) (x (j) ) = δ ij , with δ the Kronecker delta. Then the mass matrix becomes diagonal with entries M ii = e∈T x (i) ω e,x (i) ρ(x (i) ), where T x denotes the set of elements containing or adjacent to x.
For quadrilaterals and hexahedra, mass-lumping is achieved by using tensor-product basis functions and Gauss-Lobatto integration points. The resulting scheme is known as the spectral element method. For triangles and tetrahedra, mass-lumping is less straight-forward. Combining standard Lagrangian basis functions with a Newton-Cotes quadrature rule results in an efficient mass-lumped scheme for linear tetrahedra, but for higher-degree basis functions, this approach results either in an unstable scheme due to non-positive quadrature weights or in a scheme with a reduced order of convergence. This problem can be resolved by enriching the finite element space with higher-degree bubble functions and by adding integration points to the interior of the elements and faces. For example, by enriching the space of the quadratic tetrahedron with 3 degree-4 face bubble functions and 1 degree-4 interior bubble function, an enriched degree-2 mass-lumped tetrahedron that remains third-order accurate can be obtained [25].
In this paper we will analyse the standard linear mass-lumped finite element method, the mass-lumped finite element method of degree 2 derived in [25], and the 2 versions of degree 3 mass-lumped finite element methods derived in [5]. We will refer to these methods as ML1, ML2, ML3a and ML3b, respectively.

The Symmetric Interior Penalty Discontinuous Galerkin
Method. Another way to obtain a (block)-diagonal mass matrix is by allowing the finite element space U h to be discontinuous at the faces. When choosing basis functions that have support on only a single element, the mass matrix becomes block-diagonal with each block corresponding to a single element. When using orthogonal basis functions, the mass matrix even becomes strictly diagonal. In order to keep the finite element method stable and consistent with the analytic solution, the elliptic operator needs to be augmented. This can be accomplished with the symmetric interior penalty method [17], where a is replaced by the discrete (semi)-elliptic operator a where F h,in and F h,d are the internal faces and boundary faces on Γ d , respectively, α h ∈ e∈T L ∞ (∂e) is the penalty function, and { {·} }, [[·]] are the average trace operator and jump operator, respectively, defined as for all faces f ∈ F, where T f denotes the set of elements adjacent to face f , andn| ∂e denotes the outward pointing normal unit vector of element e. The bilinear form a (C) h is the same as the original elliptic operator a and is the part that remains when both input functions are continuous. The bilinear form a (D) h can be interpreted as the additional part that results from partial integration of the elliptic operator a when the first input function is discontinuous. Finally, the bilinear form a (IP ) h is the part that contains the interior penalty function needed to ensure stability of the scheme.
The penalty term can have a significant impact on the performance of the SIPDG method, since a larger penalty term results in a more restrictive bound on the time step size, but also because it can have a significant effect on the accuracy, as we will show in Section 6. Several lower bounds for the penalty term are based on the trace inequality of [36], including [34,14,28], among which we found the bound in [28] to be the sharpest. Recently, a sharper penalty term bound was presented in [16], which is based on a more involved trace inequality. In this paper we will consider both the penalty term of [16], given by (4a), and the one of [28], given by (4b): for all e ∈ T h , f ⊂ ∂e, where p denotes the degree of the polynomial basis functions, P p (e) denotes the space of polynomial functions of degree p or less in element e, ν h | ∂e∩f := |f |/|e| is a scaling function of order h −1 , with |e|, |f | the volume of e and area of f , respectively, c −1 n denotes the (pseudo)-inverse of the second-order tensor cn :=n · C ·n, wheren is the outward pointing normal unit vector, and d e denotes the diameter of the inscribed sphere of e. Although the first version requires more preprocessing time, it allows for an approximately 1.5 times larger time step [16].
We will refer to the SIPDG method with p = 1, 2, 3 using the penalty term as defined by (4a) as DG1a, DG2a, and DG3a, respectively, and to the same methods using the penalty term as defined by (4b) as DG1b, DG2b, and DG3b.

The Lax-Wendroff Time Integration Scheme.
To solve the resulting set of ODE's (3) in time, we use the Lax-Wendroff method [21,12], which is based on Taylor expansions in time and substitutes the time derivatives by matrix-vector operators using the original equations (3). For the second-order formulation, the resulting scheme is also known as Dablain's scheme [10]. The advantage of this scheme is that it is time-reversible, energy-conservative, and only requires K stages for a 2K-order of accuracy.
To introduce the scheme, let ∆t > 0 denote the time step size, and let U h (t i ) denote the approximation of u h at time t i := i∆t for i = 0, . . . , N T with N T the total number of time steps. The order-2K Lax-Wendroff method can be written as In case K = 1, this scheme reduces to the standard leap-frog or central difference scheme. When there is no source term, (5) simplifies to For the dispersion analysis, we will choose K equal to the polynomial degree p of the spatial discretization, since this will result in a 2p-order convergence rate of the dispersion error as shown in Section 6.

Dispersion Analysis
A common measure for the quality of a numerical method for wave propagation modelling is the amount of numerical dispersion and dissipation. Numerical dispersion refers in this context to the discrepancy between the numerical and physical wave propagation speed and numerical dissipation is the loss of energy in the numerical scheme. Since the schemes that we consider are all energy-conservative, they do not suffer from numerical dissipation. However, when projecting a physical wave onto the discrete space, this results in a superposition of a well-matching numerical wave and several numerical waves that have a completely different shape and frequency. We compute the number of these non-matching or spurious waves and refer to it as the eigenvector error, since it is related to the accuracy of the eigenvectors of M −1 A, while the dispersion error is related to the accuracy of the eigenvalues of M −1 A. We analyse the dispersion and eigenvector error using standard Fourier analysis, which is also known in this context as plane wave analysis. The main idea of this analysis is to compare physical plane waves with numerical plane waves on a homogeneous periodic domain, free from external forces, using a periodic mesh. To obtain a periodic tetrahedral mesh we subdivide a small cell into tetrahedra and repeat this pattern to fill the entire domain as illustrated in Figure 1. By using Fourier modes, we can then efficiently compute the numerical plane waves and their dispersion properties by solving eigenvalue problems on only a single cell.
Our analysis is similar to [11], but with the following extensions: • We extend the analysis to parallellepiped cells, since this allows for a more regular tetrahedral mesh. • We also compute the number of spurious modes that appear in the projection of the physical wave. • In the three-dimensional elastic case, there are two distinct secondary or shear waves with the same wave vector. To compute the dispersion and eigenvector error in this case, we consider the two best matching numerical plane waves.
We explain the dispersion analysis in more detail in the following subsections. First, we show how we can derive an analytical expression for the numerical plane waves using Fourier modes. After that, we show how we use this to compute the numerical dispersion and eigenvector error. In the last subsection we explain how we estimate the computational cost for each method.

5.1.
Analytic Expression for the Numerical Plane Waves. We first consider a periodic cubic domain of the form Ω := [0, N ) 3 , with N a positive integer, and later extend the results to parallelepiped domains which allow for more regular tetrahedral meshes. The physical plane wave has the following form: whereî := √ −1 is the imaginary number, κ ∈ R 3 is the wave vector, ω ∈ R is the angular velocity, and a ∈ R m is the amplitude vector. The wave vector must be of the form κ = κ z = 2π N z, with z ∈ Z 3 N , in order to satisfy the periodic boundary conditions. The numerical plane wave can be written in a similar form when using a periodic mesh. To obtain a periodic tetrahedral mesh, we subdivide the unit cell Ω 0 := [0, 1) 3 into tetrahedra and repeat this pattern N × N × N times to fill the entire domain as illustrated in Figure 1. We equip the mesh with a translation-invariant set of basis functions where each basis function has minimal support. In case of mass-lumping we use nodal basis functions and in case of DG we use basis functions that have support on only a single element. The numerical plane wave U h of the fully discrete scheme then has the form Here, U h (Ω k , t i ) denotes the coefficients of the basis functions corresponding to cell Ω k := k + Ω 0 at time t i . In case of mass-lumping, these basis functions are the nodal basis functions corresponding to the nodes on Ω k = k + [0, 1) 3 , while in case of DG, these are the basis functions that have support on one of the tetrahedra in Ω k . The vector U h,Ω 0 ∈ R n 0 denotes the basis function coefficients corresponding to cell Ω 0 at time 0 and x k = k are the coordinates of the front-left-bottom vertex of cell Ω k .
To show that this is indeed a numerical plane wave, let M (Ω k ,Ωm) , A (Ω k ,Ωm) ∈ R n 0 ×n 0 , for k, m ∈ Z 3 N , be submatrices of M and A, respectively, defined as follows: where {w (Ω k ,i) } n 0 i=0 denote the basis functions corresponding to cell Ω k and where a h = a (DG) h and (·, ·) h = (·, ·) in case of the DG method and a h = a and (·, ·) h = (·, ·) (L) h in case of the mass-lumped method.
Since the basis functions are translation invariant, the submatrices M (Ω k ,Ω k+∆k ) and A (Ω k ,Ω k+∆k ) are the same for any k ∈ Z 3 N with ∆k ∈ Z 3 N fixed. Furthermore, the submatrices M (Ω k ,Ωm) are only non-zero when k = m, since the mass matrix is diagonal in case of mass-lumping and block-diagonal, with each block corresponding to an element, in case of DG. The submatrices A (Ω k ,Ω k+∆k ) are only non-zero when ∆k ∈ {−1, 0, 1} 3 , since the nodal basis functions for mass-lumping and the local basis functions for DG do not interact when they are two or more cells apart. This implies that we only need to consider the submatrices M (Ω 0 ) := M (Ω 0 ,Ω 0 ) and A (Ω 0 ,Ω ∆k ) for ∆k = {−1, 0, 1} 3 . Now let κ = κ z := 2π N z, for some z ∈ Z 3 N , and let U h,0 ∈ R N 3 ×n 0 be the numerical wave at time t = 0: inv the inverse of M (Ω 0 ) and In other words, we can obtain eigenpairs of M −1 A by computing the eigenpairs of a small matrix S (κ) ∈ R n 0 ×n 0 . Note that since M (Ω 0 ) is symmetric positive definite, and A (κ) is Hermitian, S (κ) has n 0 distinct eigenpairs. Since there are N 3 choices for z ∈ Z 3 N and S (κz) has n 0 eigenpairs, we can obtain all of the N 3 × n 0 eigenpairs of M −1 A in this way. Now consider the numerical plane wave in (8) with (s h , U h,Ω 0 ) an eigenpair of S (κ) , so with (s h , U h,0 ) an eigenpair of M −1 A. We can rewrite U h as U h (t) = U h,0 e −î(ωt) . If we then substitute this wave into (6) we obtain It remains to determine the time step size ∆t. In the appendix we show that the numerical scheme is stable, if where σ max (M −1 A) denotes the spectral radius of M −1 A and c K is a constant, given by To obtain a bound on the spectral radius, recall that we can write every eigenpair of M −1 A in the form of (s h , U h,0 ), with U h,0 given in (9) and with (s h , U h,Ω 0 ) an eigenpair of S (κz) for some z ∈ Z 3 N . This implies that σ max (M −1 A) is equal to sup z∈Z 3 N σ max (S (κz) ). We can therefore bound σ max (M −1 A) as follows: the space of all distinct wave vectors κ. The constants c K can be computed numerically. For example, c K = 4, 12, 7.57 for K = 1, 2, 3, respectively. For higher values of K, see, for example, [27], where his σ t satisfies c K = 2σ t .
We can extend the results of this section to parallelepiped cells by applying a linear transformation x → T · x, with T ∈ R 3×3 a second-order tensor. The parallelepiped domain is then given by Ω = T · (0, N ) 3 and the cells are given by Ω 0 = T · [0, 1) 3 and Ω k = x k + Ω 0 , with x k = T · k the front-left-bottom vertex. The wave vectors κ z are of the form κ z = 2π N (T −t · z) and the wave vector space K 0 is given by K 0 := T −t · [0, 2π) 3 , with T −t the transposed inverse of T.

Computing the Dispersion and Eigenvector Error.
To explain how we compute the dispersion error, we first consider the acoustic wave equation. Let κ be a given wave vector and let u (κ) be the acoustic plane wave given by u (κ) (x, t) = eˆı (κ·x−ωt) . The angular velocity is given by ω = ±c|κ| with c the acoustic wave propagation speed. We compare this plane wave with the numerical plane waves.
To do this, we use the results from the previous subsection. There, we showed that for any eigenpair (s h , U h,Ω 0 ) of S (κ) we can obtain a numerical plane wave in the form of (8) with angular velocity ±ω h given by (10). Since S (κ) has n 0 eigenpairs, this means we can obtain n 0 discrete plane waves {U for a given wave vector κ. The corresponding wave propagation speeds {c |/|κ| and we can order the numerical plane waves such that We consider U (κ,1) h to be the matching numerical plane wave and U (κ,i) h , with i > 1, to be spurious modes. We then define the dispersion error as follows The complete procedure for computing e disp (κ) in the acoustic case is given by a. Compute all eigenpairs (s c. Compute the wave propagation speeds c h,0 , for i > 1, that may have a completely different shape and velocity. We can compute the number of these spurious waves by computing the projection error.

Now let u
To do this, we let U . We can then compute e vec (κ) by The complete procedure for computing e vec (κ) in the acoustic case is given by C. Compute the wave propagation speeds c For the isotropic elastic case, the procedure is very similar. Let κ be the wave vector and let u (κ) denote the elastic plane wave of the form u (κ) (x, t) = aeˆı (κ·x−ωt) , with a the amplitude vector, ω = ±c|κ| the angular velocity, and c the elastic wave propagation speed. In the elastic isotropic case, we have to distinguish between longitudinal or primary waves, where a is parallel with κ and the propagation speed is c = c P = (λ + 2µ)/ρ, and transversal, shear or secondary waves, where a is perpendicular to κ and the propagation speed is c = c S = µ/ρ.
For the analysis, we will only consider the secondary wave, since the wavelength λ = 2π/|κ| = 2πc/ω of this wave is shorter and therefore governs the required mesh resolution. In 3D, there are two linear independent amplitude vectors, a (κ,1) and a (κ,2) , that are perpendicular to κ and we will refer to the corresponding secondary plane waves as u (κ,1) and u (κ,2) . We will compare these physical plane waves with the numerical plane waves in a similar way as for the acoustic case.
Since, for a given κ and ω = ±c S |κ|, there are two linearly independent secondary waves, we compare the secondary wave velocity c = c S with the wave propagation speed of the two best matching numerical plane waves. In particular, we define the dispersion error as The procedures for computing this error is the same as for the acoustic case, with step d replaced by

The eigenvector is now computed by
onto the discrete space of cell Ω 0 , and U (κ) In other words, we compute the worst possible projection error for a linear combination of u The procedure for computing e vec (κ) is the same as for the acoustic case, with steps D-F replaced by We can use these errors to determine the required number of elements per wavelength. To compute these errors, we use a search algorithm, which requires the computation of e disp (κ) and e vec (κ) for a large number of wave vectors κ. The complete procedure for computing the dispersion and eigenvector error is given by: (1) Construct a cell Ω 0 and subdivide it into tetrahedra.
(3) Compute s h,max , given by (13). This is done with a search algorithm which requires the computation of σ max (S (κ) ), with S (κ) := M , for a large number of wave vectors κ. (4) Compute ∆t ≤ c K /s h,max , with c K given by (12). (5) For a given wavelength λ, compute the errors e disp (λ) and e vec (λ) given in (14).
For each λ, this requires the computation of e disp (κ) and e vec (κ), using steps a-d and A-F, for a large number of wave vectors κ.

Estimating the Computational Cost.
To compare the efficiency of the different methods, we also compute the number of degrees of freedom n vec , the number of nonzero entries of the stiffness matrix n mat , and the estimated computational cost n comp , for each wavelength λ. We define n vec to be the number of degrees of freedom per λ 3 -volume. This is computed by where n 0 is the number of basis functions corresponding to cell Ω 0 , and |Ω 0 | is the volume of Ω 0 . We define n mat to be the number of non-zero entries of the stiffness matrix per λ 3volume. In case of mass-lumping, we estimate this number by where |U q | is the number of degrees of freedom per node (|U q | = 1 in the acoustic and |U q | = 3 in the elastic case), Q Ω 0 is the set of nodes on Ω 0 , and N (q) are the neighbouring nodes of q that are connected with q through an element.
In case of the SIPDG method, we estimate this number by where |U e | is the number of basis functions with support on element e, T Ω 0 are the elements in Ω 0 , and N (e) are the neighbouring elements of e that are connected with e through a face.
To estimate the computational cost we look at the size of the matrix times the number of matrix-vector products. The resulting estimates gives a rough estimate of the relative CPU time of the different methods, since it estimates the number of computations when using a globally assembled matrix.
We define the computational cost n comp as the number of non-zero matrix entries per λ 3 -volume times the number of matrix-vector products during one oscillation in time. The duration of one oscillation is T 0 = λ/c, with c the wave propagation speed. The number of matrix-vector products during one oscillation is the number of stages of the Lax-Wendroff scheme K times the number of time steps N ∆t = T 0 /∆t = λ/(c∆t), where ∆t = c K /s h,max , with c K given by (12) and s h,max given by (13). We use this to compute n comp as follows:

Results and Comparisons
An overview of the different finite element methods that we analyse is given in Table  1. Each method is combined with an order-2p Lax-Wendroff time integration scheme, where p denotes the degree of the spatial discretization. Linear mass-lumped finite element method ML2 Degree-2 mass-lumped finite element method [25] ML3a, ML3b Degree-3 mass-lumped finite element methods [5] DGX Symmetric Interior Penalty Discontinous Galerkin method [17] of degree X = 1, 2, 3 DGXa DGX with penalty term derived in [16] and given by (4a) DGXb DGX with penalty term derived in [28] and given by (4b) To analyse the dispersion properties of these methods, we use standard Fourier analysis, as explained in Section 5. We consider a periodic mesh of congruent nearly-regular equifacial tetrahedra, known as the tetragonal disphenoid honeycomb. To obtain this mesh, we slice the unit cell Ω 0 := [0, 1) 3 into 6 tetrahedra with the planes x = y, x = z, and y = z and then apply the linear transformation x → T · x, with T := An illustration of this mesh is given in Figure 2. 6.1. Acoustic Waves on a Regular Mesh. We first consider the acoustic wave model with c = ρ = 1. Figure 3 illustrates the dispersion and eigenvector error with respect to the number of elements per wavelength N E := 3 λ 3 /|e| av , with λ the wavelength and |e| av the average element volume. The eigenvector error for ML1 is always zero, since it has only one degree of freedom per cell Ω k and therefore allows only one numerical plane wave for a given wave vector. From this figure we can obtain the order of convergence, which is 2p for the dispersion error and p+1 for the eigenvector error. These convergence rates are typical for symmetric finite element methods for eigenvalue problems, see, for example, [4] and the references therein. The 2p-order superconvergence of the dispersion error is also in accordance with the results of [26,2,13]. By extrapolating the results shown in Figure 3 we can obtain approximations of the errors of the form e = α(N E ) −β , where α is the leading constant and β is the order of convergence. The approximations are given in Table 2.  We can use these results to obtain estimates for the number of elements per wavelength required for a given accuracy, but we can also use them to obtain other properties, such as the number of time steps or the computational cost required for a given accuracy. An overview for a dispersion error of 0.01 and 0.001 is given in Table 3 and 4, respectively, and the relation between the accuracy and the computational cost is illustrated in Figure  4. Figure 4 shows that for linear elements, the mass-lumped method ML1 is significantly more efficient than the DG methods DG1a and DG1b, while for quadratic elements,   DG2a is significantly more efficient than ML2 and DG2b, and for cubic functions, DG3a is slightly more efficient than ML3b and significantly more efficient than DG3b and ML3a. In all cases, the DG methods using the sharper penalty term given by (4a) are significantly more efficient than those using the penalty term given by (4b). For a dispersion error of around 0.01 and higher, the linear mass-lumped method ML1 performs best in terms of computational cost, while for a dispersion error below 0.001 the best Tables 3 and 4 also show that for the case p = 1, the eigenvector error is always smaller than the dispersion error, but that for higher-order elements, the eigenvector error can become almost 5 times as large when the dispersion error is 0.01 and 10 times as large when the dispersion error is 0.001. This is due to the fact that the dispersion error converges with a faster rate (order 2p) than the eigenvector error (order p + 1) for higher-degree methods.
6.2. The Effect of Mesh Distortions. We also investigate the effect of the mesh quality on the dispersion error. To do this, we first create meshes of very flat elements by scaling the regular disphenoid mesh in the z-direction. After that, we create distorted meshes by displacing some of the vertices of the disphenoid honeycomb. To create flat elements, we scale the disphenoid mesh in the z-direction by a factor T z . The effect on the dispersion error is illustrated in Figure 5. For a mesh flattened by a factor 2, the dispersion error does not grow more than a factor 2.5, but flattening the mesh by a factor 10 increases the error by a factor between 10 and 100. In all cases, the mesh resolution remains the same and even becomes smaller in the z-direction. This means that the mesh quality can have a strong effect on the accuracy of the method and that using flat tetrahedra can significantly reduce the accuracy. The methods using lower-order elements are more sensitive to the mesh quality than the higher-order methods.  To create distorted meshes, we displace some of the vertices of the disphenoid mesh. In particular, we create a distorted mesh using the following steps: (1) Slice the cube [0, 0.5) 3 into 6 tetrahedra with the planes x = y, x = z, y = z.
(2) Repeat this pattern 2 × 2 × 2 times to pack the unit cell [0, 1) with 48 tetrahedra.  (4) Apply the transformation x → T · x, with T defined as in (15). In case of zero distortion, δ = 0, we obtain the original disphenoid honeycomb, scaled by a factor 0.5. When the distortion δ approaches 1, some of the elements become completely flat with zero volume.
An illustration of the mesh with distortion δ = 0.4 is given in Figure 6. In Figure  7, the dispersion and eigenvector error are plotted against the number of elements per wavelength for a heavily distorted mesh with δ = 0.9. These results show that the order of convergence remains 2p for the dispersion and p + 1 for the eigenvector error, even though the mesh is distorted. The distortion does, however, affect the leading constant of the errors. The effect of the mesh distortion on the dispersion error is illustrated in Figure 8. Again, the accuracy is not significantly affected by small distortions, but large distortions can reduce the accuracy by an order of magnitude.
6.3. Elastic Waves and the Effect of the P/S-wave Velocity Ratio. Besides the acoustic wave model, we also consider the isotropic elastic wave model. Figure 9 illustrates the dispersion and eigenvector error with respect to the number of elements per wavelength for the isotropic elastic wave model with µ = ρ = 1 and λ = 2, so with a P/S-wave velocity ratio of 2. Again, the order of convergence is 2p for the dispersion error and p + 1 for the eigenvector error.
By extrapolating these results we can again obtain approximations of the errors of the form e = α(N E ) −β , which are given in Table 5. Figure 10 illustrates the relation between the dispersion error and the computational cost, based on these results. The relative performance of the different methods is similar to the acoustic case.
We also look at the influence of the P/S-wave velocity ratio c P /c S on the dispersion error, where c S = √ µ denotes the S-wave velocity and c P = √ λ + 2µ denotes the P-wave velocity. This relation is illustrated in Figure 11. This figure shows that the DG methods are not really sensitive to the c P /c S ratio, since the dispersion error never grows more than a factor 1.5. The higher-order mass-lumped methods are slightly more sensitive, with a dispersion error becoming around 3 times as large for c P /c S = 10, compared to c P /c S = 2, while the linear mass-lumped method is very sensitive, with a dispersion error becoming almost 40 times as large in this case.

Conclusions
We analysed the dispersion properties of two types of explicit finite element methods for modelling wave propagation on tetrahedral meshes, namely mass-lumped finite elements methods and symmetric interior penalty discontinuous Galerkin (SIPDG) methods, both for degrees p = 1, 2, 3 and combined with an order-2p Lax-Wendroff time integration method. The analysed methods are listed in Table 1.
The dispersion properties are obtained semi-analytically using standard Fourier analysis. We used this to give an indication of which method is the most efficient for a given accuracy, how many elements per wavelength are required for a given accuracy, and how sensitive the accuracy of the method is to poorly shaped elements and high P/S-wave velocity ratios.
Based on the results we draw the following conclusions with regard to efficiency:   Figure 11. Relative dispersion error for the isotropic elastic wave model with different c P /c S ratios. Here, e disp,0 denotes the error for the original mesh with c P /c S = 2.
• The linear mass-lumped method is the most efficient method for a dispersion error of around 1% when using approximately regular tetrahedra. Heavily distorted elements, however, can significantly reduce its accuracy. • The degree-3 SIPDG method, with the penalty term derived in [16] and given by (4a), and the second degree-3 mass-lumped finite element method of [5] are the most efficient methods for a dispersion error of around 0.1% and less. • The SIPDG methods using the sharper penalty term bound derived in [16] are significantly more efficient than those using the penalty term of [28], which is based on the trace inequality of [36].
The required number of elements for a given accuracy can be obtained from the approximations given in Tables 2 and 5. We also draw the following conclusions with regard to accuracy: • Higher-order methods suffer more from spurious modes for the same dispersion error. This is due to the fact that for higher-order methods, the convergence rate of the dispersion error, 2p, is larger than the convergence rate of the eigenvector, p + 1. • All methods are significantly affected by a poor mesh quality, although lowerorder methods are more sensitive to this than higher-order methods. Flattening the tetrahedra by a factor 10 reduces the accuracy of the methods by 1-2 orders of magnitude, even though the mesh resolution remains the same and even improves in one direction. • The SIPDG methods are not really sensitive to high P/S-wave velocity ratios, while the accuracy of the higher-order mass-lumped methods reduces slightly when the P/S-wave velocity ratio is increased. The accuracy of the linear masslumped method, however, reduces by an order of magnitude when the P/S-wave velocity ratio is raised from 2 to 10.