New Higher-Order Mass-Lumped Tetrahedral Elements for Wave Propagation Modelling
S. Geevers W. A. Mulder J. J. W. van der Vegt

Two-field Warm Inflation and Its Scalar Perturbations on Large Scales
    Jian-Yang Zhu
    December 30, 2023
=====================================================================




We present a new accuracy condition for the construction of continuous mass-lumped elements. This condition is less restrictive than the one currently used and enabled us to construct new mass-lumped tetrahedral elements of degrees 2 to 4. The new degree-2 and degree-3 tetrahedral elements require 15 and 32 nodes per element, respectively, while currently, these elements require 23 and 50 nodes, respectively. The new degree-4 elements require 60, 61 or 65 nodes per element. Tetrahedral elements of this degree had not been found yet. We prove that our accuracy condition results in a mass-lumped finite element method that converges with optimal order in the L^2-norm and energy-norm. A dispersion analysis and several numerical tests confirm that our elements maintain the optimal order of accuracy and show that the new mass-lumped tetrahedral elements are more efficient than the current ones.








§ INTRODUCTION

Wave propagation modelling has many applications in the fields of structural mechanics, electromagnetism and geosciences. In many of these applications, waves need to be modelled on a large and complex 3D geometry that requires a fast and robust numerical algorithm. 

The oldest and most popular algorithm is the finite difference method, which approximates the wave field on a uniform grid. This method is relatively easy to implement and is very efficient on simple geometries. However, its accuracy quickly deteriorates if the grid points are not aligned with sharp material interfaces and boundaries of the domain. A good alignment is often not possible with uniform grids. 

Unstructured meshes, on the other hand, offer more geometric flexibility and can be properly aligned with many complex geometries. Such meshes can be used with finite element methods. While more difficult to implement and requiring more computations, the finite element method can remain accurate on very complex geometries when using a proper mesh. When applied with mass lumping, the finite element method can in such cases become more efficient than the finite difference method <cit.>.

Mass lumping is important for applying the finite element method to wave propagation problems, since it allows for explicit time-stepping. When using an explicit time integration scheme, the finite element method requires the solution of a linear system Mx=b, with M the mass matrix, at every time step. When using the classical finite element method, the mass matrix is large and sparse, but not (block)-diagonal. This makes the numerical scheme very inefficient for large-scale simulations. Mass lumping avoids this problem by lumping the mass matrix M into a diagonal matrix. Usually, this is done with nodal basis functions and an inexact quadrature rule for M of which the quadrature points coincide with the basis functions nodes.

For quadrilaterals and hexahedra, mass lumping is relatively straightforward and is accomplished by using tensor product basis functions and Gauss–Lobatto quadrature points. The resulting method is known as the spectral element method. Quadrilaterals and hexahedra, however, offer less geometric flexibility than triangles and tetrahedra. 

For linear triangular and tetrahedral elements, mass lumping is done using standard Lagrangian basis functions and a Newton–Cotes integration rule. For higher-degree triangular and tetrahedral elements, however, this approach results in instabilities, a singular mass matrix, or a suboptimal convergence rate. The Newton–Cotes rule for quadratic triangular elements, for example, has zero weights at the vertices, resulting in a singular mass matrix. This can be resolved by enriching the quadratic element space with a cubic bubble function that vanishes on all edges and by adding an additional node at the centre of the triangle <cit.>. By enriching the element space with higher-degree bubble functions and combining it with a suitable quadrature rule, mass-lumped triangular elements were also obtained for degrees 3 <cit.>, 4 <cit.>, 5 <cit.>, 6 <cit.>, and 7 to 9 <cit.>. For tetrahedra, mass lumping can be accomplished in a similar way by adding higher-degree face and internal bubble functions to the element space. So far, this has resulted in mass-lumped tetrahedral elements of degrees 2 <cit.> and 3 <cit.>.

In this paper we show that the accuracy condition that was imposed on the quadrature rules of these higher-degree triangular and tetrahedral mass-lumped elements is too strong. This condition is that the quadrature rule of a degree-p element should be exact for polynomials up to degree p+p'-2 <cit.>, where p'>p is the highest polynomial degree of the functions in the enriched element space. Instead, we show that for p≥ 2 the quadrature rule only needs to be exact for functions in Ũ⊗𝒫_p-2, with Ũ the enriched element space and 𝒫_p-2 the set of polynomials up to degree p-2. We prove that by satisfying this condition, the finite element method can maintain an optimal order of convergence in the L^2-norm and energy-norm.

This new accuracy condition enabled us to develop several new mass-lumped tetrahedral elements of degrees 2 to 4. The new elements of degree 2 and 3 require 15 and 32 nodes per element, respectively, while the current versions require 23 and 50 nodes, respectively. Our degree-4 elements require 60, 61 or 65 nodes. Mass-lumped tetrahedral elements of this degree had not been found yet. A dispersion analysis and various numerical tests confirm the optimal order of convergence of these methods and show that the new mass-lumped tetrahedral elements are significantly more efficient than the current ones.

Although this paper focuses on wave propagation problems, more generally, mass lumping is useful for solving any type of evolution problem that requires explicit time-stepping. It is also useful for efficiently computing higher-order derivatives, which appear, for example, in the Korteweg–de Vries equation <cit.>.


This paper is constructed as follows: In Section <ref>, we present the scalar wave equation and the classical finite element method. In Section <ref>, we explain mass lumping. The stability is analyzed in Section <ref>. In Section <ref>, we present our new accuracy condition for the quadrature rule for the mass matrix and prove that, if this condition is satisfied, the mass-lumped finite element method can maintain an optimal order of convergence. This condition enabled us to derive several new mass-lumped tetrahedral elements of degrees 2 to 4, presented in Section <ref>. We analyze the dispersion properties of these new methods in Section <ref> and test the methods numerically in Section <ref>. In both sections we compare the new methods with existing finite element methods. Finally, we present our main conclusions in Section <ref>.






§ THE SCALAR WAVE EQUATION AND CLASSICAL FINITE ELEMENT METHOD


In this paper, we mainly focus on the scalar wave equation, which serves as a model problem for more complex wave problems such as the elastic wave equations and Maxwell's equations. Let Ω⊂ℝ^3 be a three-dimensional open bounded domain, with Lipschitz boundary ∂Ω, and let (0,T) be the time domain. The scalar wave equation can be written as


    ρ∂_t^2 u    = ∇· c∇ u + f       in Ω× (0,T), 
    
    u    =0       on ∂Ω, 
    
    u|_t=0   = u_0       in Ω, 
    ∂_tu|_t=0   = v_0       in Ω,



where u:Ω×(0,T)→ℝ is the unknown scalar field, ∇ is the gradient operator, ρ,c:Ω→ℝ^+ are positive scalar fields, and f:Ω×(0,T)→ℝ is the source term. We assume that the parameters ρ and c are bounded by ρ_0≤ρ≤ρ_1 and c_0≤ c≤ c_1 for some positive scalars ρ_0,ρ_1,c_0,c_1∈ℝ^+.

This equation can be solved with the finite element method, which is based on the weak formulation of (<ref>). Assume the initial conditions satisfy u_0∈ H^1_0(Ω) and v_0∈ L^2(Ω) and assume the source term satisfies f∈ L^2(0,T;L^2(Ω)). Here, L^2(Ω) denotes the space of square integrable functions on Ω, H^1_0 denotes the Sobolev space of functions on Ω that are zero on ∂Ω and have square-integrable weak derivatives, and L^2(0,T;U), with U a Banach space, denotes the Bochner space consisting of functions f:(0,T)→ U such that f(t)_U is square integrable in (0,T). The weak formulation of (<ref>) is finding u∈ L^2(0,T;H^1_0(Ω)), with ∂_tu∈ L^2(0,T;L^2(Ω)) and ∂_t(ρ∂_tu)∈ L^2(0,T;H^-1(Ω)), such that u|_t=0=u_0, ∂_tu|_t=0=v_0, and

    ⟨∂_t(ρ∂_t u), w⟩ + a(u,w)    = (f,w)        for all w∈ H^1_0(Ω),  a.e. t∈(0,T).

Here, ⟨·,·⟩ denotes the pairing between H^-1(Ω) and H_0^1(Ω), (·,·) denotes the L^2(Ω) inner product, and a(·,·):H^1_0(Ω)× H^1_0(Ω)→ℝ is the elliptic operator given by

    a(u,w):=∫_Ω c∇ u·∇ w  dx.

Because of the boundedness of ρ, it follows that the norm u_ρ^2:=(ρ u,u) is equivalent to the standard L^2(Ω)-norm. It can then be proven, in a way analogous to <cit.>, that (<ref>) is well-posed and has a unique solution.







The solution of (<ref>) can be approximated by the finite element method. Let 𝒯_h be a tetrahedral tesselation of Ω, with h the diameter of the smallest sphere that can contain each element in 𝒯_h, and let U_h  denote the finite element space consisting of continuous functions that are polynomial of degree at most p when restricted to a single element:

    U_h ={u∈ H^1_0(Ω)  |  u|_e∈𝒫_p(e)  for all e∈𝒯_h},

where 𝒫_p denotes the set of all polynomials of degree p or less. The classical conforming finite element method is finding u_h:[0,T]→ U_h, such that u_h|_t=0=Π_hu_0, ∂_tu_h|_t=0=Π_hv_0, and 

    (ρ∂_t^2 u_h,w) + a(u_h,w)    = (f,w)       for all w∈ U_h,  a.e. t∈(0,T),

where Π_h:L^2(Ω)→ U_h is the weighted L^2 projection operator defined such that (ρΠ_hu,w)=(ρ u,w) for all w∈ U_h. 

This can be rewritten as a set of ODE's using a linear basis {w_i}_i=1^n of U_h. For any function u∈ U_h we define u∈ℝ^n as the vector of coefficients such that u = ∑_i=1^n u_i w_i. The finite element method can then be formulated as solving _h:[0,T]→ℝ^n, such that _h|_t=0=Π_hu_0, ∂_t_h|_t=0=Π_hv_0, and

    M∂_t^2_h + A_h = f^*       for a.e. t∈(0,T),

where M,A∈ℝ^n× n are the mass matrix and stiffness matrix, respectively, given by M_ij := (ρ w_i,w_j), A_ij := a(w_i,w_j) for all i,j=1,…,n, and f^*∈ L^2(0,T;ℝ^n) is the source vector, given by f_i^*:=(f,w_i), for i=1,…,n, a.e. t∈(0,T).

When using an explicit time integration scheme, a system of the form Mx=b needs to be solved at every time step. Typically, the mass matrix M is large and sparse, but not (block)-diagonal, resulting in a very inefficient numerical scheme. A diagonal mass matrix can be obtained by a technique known as mass lumping. We will discuss this in the next section. 





§ MASS LUMPING


Mass lumping is usually done with nodal basis functions and an inexact quadrature rule for the mass matrix. A diagonal matrix is obtained when the integration points coincide with the nodes of the basis functions. However, when using elements of degree p≥ 2, this technique does not result in a stable and accurate finite element scheme. For example, for standard quadratic Lagrangian basis functions combined with a Newton–Cotes quadrature rule, the weights at the vertices of the quadratic tetrahedral element become negative, resulting in unstable modes. 

To overcome such problems, the elements are enriched with higher-degree face and interior bubble functions. These enriched elements are still affine-equivalent to a reference element ẽ. We can therefore write the discrete space in the form 

    U_h=H_0^1(Ω) ∩ U(𝒯_h, Ũ),

where

    U(𝒯_h,Ũ)    := {u∈ H^1(Ω)  |  u∘ϕ_e ∈Ũ, for all e∈𝒯_h },

with ϕ_e:ẽ→ e the reference-to-physical element mapping, and Ũ the reference space.  If Ũ=𝒫_p(ẽ) we obtain the standard elements of degree p. To obtain enriched elements, we set Ũ=𝒫_p(ẽ)⊕Ũ^+:={u  |  u=w+u^+  for some w∈𝒫_p(ẽ), u^+∈Ũ^+}, with Ũ^+ a space of higher-degree face and interior bubble functions.










A nodal basis and quadrature rule for U_h can be constructed from a nodal basis and quadrature rule for the reference space Ũ. In the next two subsections we will discuss this in more detail.




 §.§ Nodes and Nodal Basis Functions


A nodal basis for a space U_h consists of a set of nodes 𝒬_h and corresponding basis functions {w_}_∈𝒬_h, such that span{w_}_∈𝒬_h = U_h and w_() = δ_, for all ,∈𝒬_h, where δ_ denotes the Kronecker delta. This means that each basis function equals one at one particular node and zero at all the other nodes.

A common way to construct such a nodal basis for the space U(𝒯_h,Ũ) is using a nodal basis {w̃_}_∈𝒬̃ for the reference space Ũ. The element nodes 𝒬_e are obtained by mapping the reference nodes to the physical element: 𝒬_e:={ϕ_e()}_∈𝒬̃. The nodal basis functions of this element, {w_e,}_∈𝒬_e, are obtained by mapping the reference basis functions to the physical element. We can write these functions as w_e,:=w̃_ϕ_e^-1()∘ϕ_e^-1. The set of global nodes 𝒬_h is the union of all element nodes and the corresponding global basis functions are obtained by concatenating the corresponding element basis functions. Formally, we define the global nodal basis functions {w_}_∈𝒬_h as follows:

    w_|_e   := w̃_ϕ_e^-1()∘ϕ_e^-1,      e∈𝒯_, 
    
    0,    otherwise,

for all ∈𝒬_h, where 𝒯_ denotes the set of elements containing or adjacent to . To ensure that these global basis functions are well-defined and continuous, we need to impose the following additional conditions on 𝒬̃ and {w̃_}_∈𝒬̃:

    w̃_|_f̃   = 0       for all f̃∈ℱ̃, ∈𝒬̃∖f̃,

and


    𝒬̃   = s(𝒬̃)       for all s∈𝒮,  
    w̃_   = w̃_s()∘ s       for all s∈𝒮,



where ℱ̃ is the set of reference faces and 𝒮 is the set of all affine mappings that map ẽ onto itself. Condition (<ref>) implies that if a basis function is zero at the nodes on a face, then it should be zero on the entire face, and condition (<ref>) implies that the set of element nodes and basis functions are symmetric and do not depend on the choice of ϕ_e. A proof that {w_}_∈𝒬_h is indeed a set of well-defined and continuous nodal basis functions is given in Lemma <ref> and Theorem <ref>. 



It remains to incorporate the Dirichlet boundary condition u_h|_∂Ω=0. If u_h∈ U(𝒯_h,Ũ)=span{w_}_∈𝒬_h, then, because of (<ref>), this condition is satisfied when u_h=0 at all nodes on ∂Ω. A nodal basis for U_h therefore consists of all interior nodes 𝒬_h∖∂Ω and corresponding basis functions {w_}_∈𝒬_h∖∂Ω. 





 §.§ Quadrature Rule


To obtain a diagonal mass matrix, we approximate the integrals with an inexact quadrature rule of which the integration points coincide with the nodes of the nodal basis.

Let 𝒬_e:={ϕ_e()}_∈𝒬̃ be the set of nodes on e, and let {ω_e,}_∈𝒬_e be a set of corresponding weights. Together, the weights and nodes form a quadrature rule for the element. The quadrature rule is used to approximate the integrals of the mass matrix at the element as follows:

    (ρ u,w)_e = ∫_e ρ uw dx ≈∑_∈𝒬_eω_e,ρ_e()u()w()=:(ρ u,w)_𝒬_e,

where ρ_e:=ρ|_e denotes the scalar field ρ restricted to element e. We assume that ρ is continuous within each element, which implies that the approximation above is well defined. The global product (ρ u,w) is then approximated by

    (ρ u,w) ≈ (ρ u,w)_𝒬_h :=  ∑_e∈𝒯_h (ρ u,w)_𝒬_e.


Now let w_,w_, with ,∈𝒬_h, be nodal basis functions as described in the previous subsection. The corresponding mass matrix entry is given by

    (ρ w_,w_)_𝒬_h =δ_∑_ e∈𝒯_ω_e,ρ_e(),

This implies that the mass matrix is diagonal with entries of the form ∑_ e∈𝒯_ω_e,ρ_e().

The quadrature rules can be constructed from a reference quadrature rule. This rule consists of the reference nodes 𝒬̃ and a set of weights {ω̃_}_∈𝒬̃ and approximates integrals on the reference element as follows:

    ∫_ẽρ̃ũw̃ dx̃≈∑_∈𝒬̃ω_ρ̃()ũ()ṽ() =: (ρ̃ũ,w̃)_𝒬̃.

We can use this to approximate the integral of the physical element by

    (ρ u,w)_e = |e|/|ẽ|∫_ẽρ̃ũw̃ dx̃≈|e|/|ẽ|  (ρ̃ũ,w̃)_𝒬̃,

with |e| the volume of e, |ẽ| the volume of ẽ, and ρ̃:=ρ∘ϕ_e, ũ:=u∘ϕ_e, w̃:=w∘ϕ_e. This approximation is the same as (<ref>) when ω_e,=(|e|/|ẽ|)ω̃_ϕ_e^-1().

Now that we have introduced the quadrature rules for the mass matrix, we can present the mass-lumped finite element method.





 §.§ Mass-Lumped Finite Element Method

Assume ρ∈𝒞^0(𝒯_h), u_0∈ H^1_0(Ω)∩𝒞^0_0(Ω), v_0∈𝒞^0_0(Ω), and f∈ L^2(0,T;𝒞^0(Ω)). Here, 𝒞^0(𝒯_h) denotes the set of functions that are in 𝒞^0(e) when restricted to e. The mass-lumped finite element method is finding u_h:[0,T]→ U_h, such that u_h|_t=0=I_hu_0, ∂_tu_h|_t=0=I_hv_0, and 

    (ρ∂_t^2 u_h,w)_𝒬_h + a(u_h,w)    = (f,w)_𝒬_h      for all w∈ U_h,  a.e. t∈(0,T),

where I_h:𝒞^0(Ω)→ U(𝒯_h,Ũ) denotes the interpolation of a continuous function by a function in U(𝒯_h,Ũ) through the nodes of 𝒬_h.

To write this as a set of ODE's, let {^(i)}_i=1^n=𝒬_h∖∂Ω be a numbering of all interior nodes, and define w_i:=w_^(i) for all i=1,2,…,n. Then the mass-lumped finite element method can be formulated as solving _h:[0,T]→ℝ^n such that _h|_t=0=I_hu_0, ∂_t_h|_t=0=I_hv_0, and 

    M∂_t^2_h + A_h = f^*       for a.e. t∈(0,T),

where M_ij := (ρ w_i,w_j)_𝒬_h, A_ij := a(w_i,w_j) for all i,j=1,…,n, and f_i^*:=(f,w_i)_𝒬_h, for i=1,…,n, a.e. t∈(0,T). From (<ref>) it follows that M is now a diagonal matrix that can be written as

    M_ij   = δ_ij∑_e∈𝒯_^(i)ω_e,^(i)ρ_e(^(i)),       i,j=1,…,n.

This set of ODE's can be efficiently solved using an explicit time integration scheme such as the second-order leap-frog scheme or a higher-order Dablain scheme <cit.>, which is a type of Lax–Wendroff scheme <cit.> for second-order wave equations.

In the next sections we analyze the stability and accuracy of the mass-lumped finite element method and derive conditions for the quadrature rules.





 §.§ Stability of the Mass-Lumped Finite Element Method


To analyze the stability of the mass-lumped finite element method, we look at the behavior of the discrete energy. Consider the mass-lumped method given in (<ref>) and substitute w=∂_t u to obtain

    ∂_t E_h    =(f,∂_t u)_𝒬_h      for a.e. t∈(0,T),

where E_h:=1/2(ρ∂_tu,∂_tu)_𝒬_h+1/2a(u,u) is the discrete energy. This implies that the discrete energy remains bounded when the source term f is bounded and that the discrete energy is conserved when there is no source term. 

For stability it then remains to show that the discrete energy is a well-defined energy. This means that (ρ v,v)_𝒬_h+ a(u,u)>0 for all u,v∈ U_h, (u,v)≠ 0, which is the case when (ρ u,u)_𝒬_h>0 for any u∈ U_h, u≠ 0. Since we can write (ρ u,u)_𝒬_h=^tM, this is satisfied when M is positive definite. From (<ref>) it follows that this is the case when all weights of the quadrature rules are strictly positive, which is the case when the weights of the reference quadrature rule are strictly positive.





§ ACCURACY OF THE MASS-LUMPED FINITE ELEMENT METHOD




 §.§ A Less Restrictive Condition on the Accuracy of the Quadrature Rule


Let U_h=U(𝒯_h,Ũ), with Ũ=𝒫_p(ẽ)⊕Ũ^+, be the finite element space constructed as in Section <ref>, where p≥ 2 denotes the degree of the finite element method and Ũ^+⊂𝒫_p'(ẽ) is the space of higher-degree face and interior bubble functions.  Also, let the quadrature rule for the mass matrix be based on a reference element quadrature rule as described in Section <ref>. We will prove that an optimal convergence rate of the mass-lumped finite element method is obtained when all weights of the reference quadrature rule, {ω̃_}_∈𝒬̃, are strictly positive and

    ∫_ẽf̃ dx̃   =  ∑_∈𝒬̃ω̃_f̃()       for all f̃∈𝒫_p-2(ẽ)⊗Ũ,

where 𝒫_p-2(ẽ)⊗Ũ:={f  |  f=wu  for some w∈𝒫_p-2(ẽ), u∈Ũ}. This means that the quadrature rule of the reference element should be exact for products of the reference basis functions and polynomials of degree p-2. Until now, the condition used for the accuracy of the quadrature rule was

    ∫_ẽf̃ dx̃   =  ∑_∈𝒬̃ω̃_f̃(),       for all f̃∈𝒫_p+p'-2(ẽ),

see for example <cit.>, so it was imposed that the reference quadrature rule should be exact for functions in 𝒫_p+p'-2(ẽ), with p' the highest polynomial degree of the enriched space, which turns out to be significantly more restrictive for tetrahedral elements. By using (<ref>) instead of (<ref>) we are able to develop new mass-lumped elements that require significantly less nodes.

In the next subsections we will prove that the convergence rate of the mass-lumped finite element method remains optimal under the less severe condition (<ref>). The novel part of the proofs are the bounds on the integration error, derived in Section <ref>. This is the only part where we explicitly use condition (<ref>). Using these bounds we can prove optimal convergence in a rather standard way.





 §.§ Some Norms and Interpolation Properties

For the convergence analysis, we use multiple interpolation properties, which we will present in this subsection. Also, to make the analysis more readable, we will use C to denote some positive constant that may depend on the regularity of the mesh, the reference space Ũ, the reference quadrature rule, the domain Ω, and the parameters ρ,c, but does not depend on the mesh resolution h, the time interval (0,T), or the choice of the functions that appear in the inequality.

Let H^k(Ω), with k≥ 1, denote the Sobolev space, consisting of functions with square integrable order-k weak derivatives equipped with norm

    u_k^2    := ∑_||≤ kD^ u_0^2,       k≥ 1,

where ·_0 denotes the standard L^2(Ω)-norm, and D^:=∂_1^α_1∂_2^α_2∂_3^α_3 denotes a higher-order partial derivative of order ||:=α_1+α_2+α_3. Also let H^k(𝒯_h), with k≥ 1, denote the broken Sobolev space, consisting of functions that belong to H^k(e) when restricted to element e, for all e∈𝒯_h. We equip this space with the norm

    u_𝒯_h,k^2    :=  ∑_e∈𝒯_hu_e,k^2 := ∑_e∈𝒯_h(∑_||≤ kD^_e^2),       k≥ 1.

Now let I_h:𝒞^0(Ω)→ U(𝒯_h,Ũ) denote the interpolation by a function in U(𝒯_h,Ũ) through the nodes of 𝒬_h. This interpolation operator is well-defined for functions in H^2(𝒯_h)∩ H^1(Ω), since H^2(e) ⊂𝒞^0(e) when e is a three-dimensional element, and therefore H^2(𝒯_h)∩ H^1(Ω)⊂𝒞^0(Ω). For this interpolation operator, we can present the following approximation properties:


Let p≥ 2 be the degree of the finite element space and let u∈  H^1(Ω)∩ H^k(𝒯_h) with k≥ 2. Then

    u-I_hu_𝒯_h,l   ≤ Ch^min(p+1,k)-lu_𝒯_h,min(p+1,k),        l≤min(p+1,k).



This result follows from <cit.>.


Now assume that the weights for the reference quadrature rule are all strictly positive. For any function in H^1(Ω)∩ H^2(𝒯_h), we can then define the following discrete L^2 semi-norm:

    |u|_𝒬_h^2    := (u,u)_𝒬_h.

This discrete semi-norm is well defined, since H^1(Ω)∩ H^2(𝒯_h)⊂𝒞^0(Ω) as mentioned before. This becomes a full norm, ·_𝒬_h, that is equivalent to the L^2-norm, for functions in U_h:


If all the weights of the reference quadrature rule are strictly positive, then

    C^-1u_0    ≤u_𝒬_h≤ Cu_0      for all u∈ U_h.




Since the function space of the reference element Ũ:=span{w̃_}_∈𝒬̃ is finite-dimensional, and since all weights of the reference quadrature rule {ω̃_}_∈𝒬̃ are positive, there exists a constant C>0 depending on the reference quadrature rule and function space Ũ, such that

    C^-1u_ẽ   ≤ũ_𝒬̃≤ Cũ_ẽ      for all ũ∈Ũ.

where ũ_𝒬̃^2 := (ũ,ũ)_𝒬̃. Then (<ref>) follows from the relations

    u_0^2 = ∑_e∈𝒯_h|e|/|ẽ|ũ_e _ẽ^2,     u_𝒬_h^2 = ∑_e∈𝒯_h|e|/|ẽ|ũ_e _𝒬̃^2,

where ũ_e:=u∘ϕ_e.


Now let Π_h,q:L^2(Ω)→𝒫_q(𝒯_h) denote the L^2-projection onto the space of piecewise nonconforming polynomials of at most degree q: 

    𝒫_q(𝒯_h):={u∈ L^2(Ω)  |  u|_e∈𝒫_q(e)  for all e∈𝒯_h}.

We then present the following interpolation properties:


Let  u∈ H^k(𝒯_h) with k≥ 2, and let q≥ 0. Then

    u-Π_h,qu_0   ≤ Ch^min(q+1,k)u_𝒯_h,min(q+1,k).

Furthermore, if also u∈ H^1(Ω), if p≥max(q,2) is the degree of the finite element space, and if all the weights of the reference quadrature rule are strictly positive, then

    |u-Π_h,qu|_𝒬_h   ≤ Ch^min(q+1,k)u_𝒯_h,min(q+1,k).



The first inequality, (<ref>), follows from <cit.>. The second inequality can be derived as follows: 

    |u-Π_h,qu|_𝒬_h   =  |I_hu-Π_h,qu|_𝒬_h
       ≤ C I_hu-Π_h,qu_0
       ≤ C(I_hu-u_0 + u-Π_h,qu_0) 
       ≤  Ch^min(q+1,k)u_𝒯_h,min(q+1,k)

where we used Lemma <ref> in the second line, the triangle inequality in the third line, and Lemma <ref> and (<ref>) in the last line.






 §.§ Bounds on the Integration Error


In this section we will derive some useful bounds on the error of the quadrature rules for the mass matrix. The proofs of these bounds will be the only cases where we explicitly use the accuracy condition of the quadrature rule, given in (<ref>). Using these results we can prove optimal order of convergence of the mass-lumped finite element method in a rather standard way.

Let u,w∈ H^2(𝒯_h), and let r_h(u,w):=(u,w)-(u,w)_𝒬_h be the integration error of the mass matrix. We can derive the following bounds on r_h:


Let p≥ 2 be the degree of the finite element space, u∈ H^k(Ω) with k≥ 2, and w∈ U_h. If the reference quadrature rule satisfies (<ref>) and if all its weights are strictly positive, then

    |r_h(u,w)| ≤ Ch^min(p,k)u_min(p,k)w_1.

and

    |r_h(u,w)| ≤ Ch^min(p+1,k)u_min(p+1,k)w_𝒯_h,2.




Using (<ref>) and the fact that 𝒫_p-2(ẽ)⊗Ũ⊃𝒫_p(ẽ) for p≥ 2, we can write

    r_h(u,w)    = r_h((u-Π_h,p-1u) + (Π_h,p-1u-Π_h,p-2u) + Π_h,p-2u,  
       =    (w-Π_h,0w) + Π_h,0w  ) 
       = r_h(u-Π_h,p-1u, w) +  r_h(Π_h,p-1u-Π_h,p-2u, w-Π_h,0w).

From this, the Cauchy–Schwarz inequality, and Lemma <ref>, we can then obtain (<ref>). 

Using (<ref>), we can also write

    r_h(u,w)    = r_h([(u-Π_h,pu) + (Π_h,pu-Π_h,p-1u) + (Π_h,p-1u-Π_h,p-2u) + 
       =    Π_h,p-2u], [(w-Π_h,1w) + (Π_h,1w - Π_h,0w)+ Π_h,0w] ) 
       = r_h(u-Π_h,pu, w) + r_h(Π_h,pu-Π_h,p-1u, w-Π_h,0w) + 
       =     r_h(Π_h,p-1u-Π_h,p-2u, w-Π_h,1w).

From this, the Cauchy–Schwarz inequality, and Lemma <ref>, we can then obtain (<ref>). 






 §.§ Optimal Convergence for a Related Elliptic Problem

To prove optimal convergence of the mass-lumped finite element method, we first prove optimal convergence for a related elliptic problem.

Let v∈ H^2(𝒯_h). The elliptic problem related to (<ref>), is finding u∈ H_0^1(Ω) such that

    a(u,w)    = (v,w)       for all w∈ H_0^1(Ω).

This problem is well defined since a is coercive and bounded with respect to the H^1(Ω)-norm, which follows from the boundedness of c and Poincaré's inequality. 

The related mass-lumped method for solving this problem is finding u_h∈ U_h such that

    a(u_h,w)    = (v,w)_𝒬_h      for all w∈ U_h.


In the next theorems we prove optimal convergence of this method in the H^1-norm and L^2-norm.



Let u be the solution of (<ref>) and u_h the solution of (<ref>), with p≥ 2 the degree of the finite element space. Also, let k_u,k_v≥ 2, u∈ H^k_u(Ω), and v∈ H^k_v(Ω). If the reference quadrature rule satisfies (<ref>) and if all its weights are strictly positive, then

    u-u_h_1 ≤ Ch^min(p,k_u-1,k_v) ( u _min(p+1,k_u) +  v _min(p,k_v)).




By definition of u and u_h, we have

    a(u-u_h,w)    = r_h(v,w),       for all w∈ U_h.

By choosing w=I_hu-u_h we can then obtain

    a(I_hu-u_h,I_hu-u_h)    = -a(u-I_h u,I_hu-u_h) + r_h(v,I_hu-u_h),

From the coercivity of a it follows that

    I_hu-u_h_1^2 ≤ C a(I_hu-u_h,I_hu-u_h).

From the boundedness of a and Lemma <ref> it follows that

    |a( u-I_h u,I_hu-u_h)| ≤ Ch^min(p,k_u-1)u_min(p+1,k_u)I_hu-u_h_1.

Using Lemma <ref> we obtain

    |r_h(v,I_hu-u_h)| ≤ Ch^min(p,k_v)v_min(p,k_v)I_hu-u_h_1.

Combining (<ref>), (<ref>), (<ref>), and (<ref>) then gives

    I_hu-u_h_1 ≤ Ch^min(p,k_u-1,k_v) ( u _min(p+1,k_u) +  v _min(p,k_v)).

From Lemma <ref> it also follows that

    u-I_hu_1 ≤ Ch^min(p,k_u-1)u_min(p+1,k_u).

Combining (<ref>) and (<ref>) then results in (<ref>).


To prove optimal convergence in the L^2-norm, we make the following regularity assumption: for any v∈ L^2(Ω), the solution u of (<ref>) is in H^2(Ω) and satisfies

    u_2≤ Cv_0.

This is certainly true if ∂Ω is 𝒞^2 and c∈𝒞^1(Ω). 




Let u be the solution of (<ref>) and u_h the solution of (<ref>), with p≥ 2 the degree of the finite element space. Also, let k_u,k_v≥ 2, u∈ H^k_u(Ω), v∈ H^k_v(Ω), and assume that the regularity condition (<ref>) holds. If the reference quadrature rule satisfies (<ref>) and if all its weights are strictly positive, then

    u-u_h_0 ≤ Ch^min(p+1,k_u,k_v) ( u _min(p+1,k_u) +  v _min(p+1,k_v))

and

    |u-u_h|_𝒬_h≤ Ch^min(p+1,k_u,k_v) ( u _min(p+1,k_u) +  v _min(p+1,k_v)).




Let z∈ H^1_0(Ω) be the solution of

    a(z,w)    = (u-u_h,w),       for all w∈ H_0^1(Ω).

From the regularity assumption it follows that z∈ H^2(Ω) and z_2≤ Cu-u_h_0. Using the definition of z, u, and u_h, we can also write

    u-u_h_0^2    = a(u-u_h,z) 
       = a(u-u_h,z-I_hz) + a(u-u_h,I_hz) 
       = a(u-u_h,z-I_hz) + r_h(v,I_hz).

Using the boundedness of a, Theorem <ref>, Lemma <ref>, and the regularity assumption, we obtain

    |a(u-u_h,z-I_hz) |  
       ≤ Cu-u_h_1 z-I_hz_1 
       ≤ Ch^min(p+1,k_u,k_v+1) (u_min(p+1,k_u) + v_min(p,k_v))u-u_h_0.

From Lemma <ref>, Lemma <ref>, and the regularity assumption it also follows that

    |r_h(v,I_hz)| ≤ Ch^min(p+1,k_v)v_min(p+1,k_v)u-u_h_0.

Combining (<ref>), (<ref>), and (<ref>) results in (<ref>).

To derive (<ref>), we use Lemma <ref> to obtain

    |u-u_h|_𝒬_h   =  I_hu-u_h_𝒬_h≤ CI_hu-u_h_0.

Combining this inequality with (<ref>) and Lemma <ref> results in (<ref>).










 §.§ Some Additional Norms and Interpolation Properties

In order to analyze the convergence for the time dependent problem, we need to introduce an additional projection operator and some additional function spaces.

Let L denote the spatial operator L:=-∇· c∇, and let u∈ H_0^1(Ω) with Lu∈𝒞^0(Ω). We define the projection π_hu∈ U_h to be the solution of

    a(π_hu,w)    = (Lu,w)_𝒬_h,       for all w∈ U_h.

We can derive the following interpolation property of this projection operator:


Let p≥ 2 be the degree of the finite element space, and let c∈𝒞^k+1(Ω) and u∈ H^1_0(Ω)∩ H^k+2(Ω), with k≥ 2. If the reference quadrature rule satisfies (<ref>) and if all its weights are strictly positive, then

    u-π_hu_1    ≤ Ch^min(p,k) u _min(p+2,k+2) ,

Moreover, if regularity condition (<ref>) also holds, then

    u-π_hu_0    ≤ Ch^min(p+1,k) u _min(p+3,k+2), 
    
    |u-π_hu|_𝒬_h   ≤ Ch^min(p+1,k) u _min(p+3,k+2) .




From partial integration it follows that a(u,w)=(Lu,w) for all w∈ H_0^1(Ω). Also, by definition of the projection we have a(π_hu,w)=(Lu,w)_𝒬_h for all w∈ U_h. The inequalities then follow from Theorem <ref> and Theorem <ref> by taking v=Lu, k_u=k+2, k_v=k, and using the bounds Lu_q≤ Cu_q+2 for q≤ k.


We also extend the Sobolev spaces H^k(Ω) to Bochner spaces L^∞(0,T;H^k(Ω)), equipped with norm

    u_∞,k   := esssup_t∈(0,T)u_k.






 §.§ Optimal Convergence of the Mass-Lumped Element Method

In this section we prove the optimal convergence of the mass-lumped finite element method for the wave equation. We first derive an equation for the behavior of the numerical error and then prove optimal convergence in the energy-norm and L^2-norm.



Let u be the solution of (<ref>) and let u_h be the solution of (<ref>). If ρ∈𝒞^0(Ω),
∂_t^2u ∈ L^2(0,T; 𝒞^0_0(Ω)), and f ∈ L^2(0,T; 𝒞^0(Ω)), then Lu∈ L^2(0,T; 𝒞^0(Ω)), and 

    (ρ∂_t^2e_h,w)_𝒬_h+ a(e_h,w)    = -(ρ∂_t^2ϵ_h, w)_𝒬_h

for all w∈ U_h and almost every t∈(0,T), where e_h:=π_hu-u_h and ϵ_h:=u-π_hu.



Since ρ is bounded and continuous, it follows that ρ∂_t^2u∈ L^2(0,T;𝒞^0_0(Ω)). Since also f∈ L^2(0,T; 𝒞^0(Ω)), it follows that Lu∈ L^2(0,T; 𝒞^0(Ω)) and ρ∂_t^2u + Lu = f. This implies

    (ρ∂_t^2u,w)_𝒬_h + (Lu,w)_𝒬_h   = (f,w)_𝒬_h

for all w∈ U_h and almost every t∈(0,T). Using the definition of π_hu we can then obtain

    (ρ∂_t^2u,w)_𝒬_h + a(π_h u,w)    = (f,w)_𝒬_h

for all w∈ U_h and almost every t∈(0,T). By definition of u_h we have

    (ρ∂_t^2u_h,w)_𝒬_h + a(u_h,w)    = (f,w)_𝒬_h
 
for all w∈ U_h and almost every t∈(0,T). Subtracting this from the previous equality and reordering the terms results in (<ref>).




Let u be the solution of (<ref>) and let u_h be the solution of (<ref>), with p≥ 2 the degree of the finite element space. Let ρ∈𝒞(Ω), f∈ L^2(0,T; 𝒞^0(Ω)), and let c∈𝒞^k+1(Ω), u,∂_tu,∂_t^2u ∈ L^∞(0,T; H^k+2(Ω)) for some k≥ 2. Also, assume that regularity condition (<ref>) holds. If the reference quadrature rule satisfies (<ref>) and if all its weights are strictly positive, then

    u-u_h_∞,1 +  ∂_tu-∂_tu_h_∞,0≤
        Ch^min(p,k)( u_∞,min(p+3,k+2) + ∂_tu_∞,min(p+3,k+2) + T∂_t^2u_∞,min(p+3,k+2)) .




Define e_h:=π_hu-u_h and ϵ_h:=u-π_hu. From Lemma <ref>, it follows that

    (ρ∂_t^2e_h,w)_𝒬_h+ a(e_h,w)    = -(ρ∂_t^2ϵ_h, w)_𝒬_h,

for all w∈ U_h and almost every t∈(0,T). By substituting w=∂_te_h we can obtain

    ∂_tE_h    = -(ρ∂_t^2ϵ_h, ∂_te_h)_𝒬_h,

for almost every t∈(0,T), where E_h:=1/2(ρ∂_te_h,∂_te_h)_𝒬_h + 1/2a(e_h,e_h) is the discrete energy. Fix T'∈(0,T) and integrate (<ref>) over (0,T') to obtain

    E_h|_t=T' = E_h|_t=0 - ∫_0^T' (ρ∂_t^2ϵ_h, ∂_te_h)_𝒬_h dt .

Using the coercivity of a, the boundedness of ρ, and Lemma <ref>, we can derive

    e_h_1 +  ∂_te_h_0   ≤ CE_h^1/2,       a.e. t∈(0,T).

From the Cauchy–Schwarz inequality, the bounds of ρ, Lemma <ref>, and Lemma <ref>, we can also obtain

    |(ρ∂_t^2ϵ_h, ∂_te_h)_𝒬_h| ≤ Ch^min(p+1,k)∂_t^2u_min(p+3,k+2)∂_te_h_0,

for almost every t∈ (0,T). Finally, we can use Lemma <ref>, Lemma <ref>, and the boundedness of ρ and a to obtain

    E_h^1/2|_t=0≤ Ch^min(p,k)(u_∞,min(p+3,k+2) + ∂_tu_∞,min(p+3,k+2)).

By taking the supremum of (<ref>) for all T'∈(0,T) and using (<ref>), (<ref>), and (<ref>), we can obtain

    e_h_∞,1 +  ∂_te_h_∞,0≤
         Ch^min(p,k)( u_∞,min(p+3,k+2) + ∂_tu_∞,min(p+3,k+2) + T∂_t^2u_∞,min(p+3,k+2)) .

Using (<ref>) and Lemma <ref> we obtain (<ref>).










Let u be the solution of (<ref>) and let u_h be the solution of (<ref>), with p≥ 2 the degree of the finite element space. Let ρ∈𝒞(Ω), f∈ L^2(0,T; 𝒞^0(Ω)), ∂_t^2u∈ L^2(0,T; 𝒞^0_0(Ω)), and let c∈𝒞^k+1(Ω), u,∂_tu,∈ L^∞(0,T; H^k+2(Ω)) for some k≥ 2. Also, assume that regularity condition (<ref>) holds. If the reference quadrature rule satisfies (<ref>) and if all its weights are strictly positive, then 

    u-u_h_∞,0≤ Ch^min(p+1,k)(u_∞,min(p+3,k+2)  + T∂_tu_∞,min(p+3,k+2)).




Define e_h:=π_hu-u_h and ϵ_h:=u-π_hu. From Lemma <ref>, it follows that

    (ρ∂_t^2e_h,w)_𝒬_h + a(e_h,w)    = -(ρ∂_t^2ϵ_h, w)_𝒬_h,

for all w∈ U_h and almost every t∈(0,T). Fix T'∈(0,T) and choose w as

    w|_t=t'   := ∫_t'^T' e_h  dt.

This implies that w|_t=T'=0 and ∂_tw=-e_h. Using the relations

    (ρ∂_t^2e_h,w)_𝒬_h   = ∂_t(ρ∂_te_h,w)_𝒬_h + 1/2∂_t(ρ e_h,e_h)_𝒬_h, 
    
    a(e_h,w)    = -1/2∂_ta(w,w), 
    
    -(ρ∂_t^2ϵ_h, w)_𝒬_h   = -∂_t(ρ∂_tϵ_h, w)_𝒬_h - (ρ∂_tϵ_h, e_h)_𝒬_h,

we can rewrite (<ref>) as

    1/2∂_t(ρ e_h,e_h)_𝒬_h   = 1/2∂_ta(w,w) -∂_t(ρ∂_t(u-u_h),w)_𝒬_h - (ρ∂_tϵ_h,e_h)_𝒬_h,
 
for almost every t∈(0,T). Integrating (<ref>) over (0,T') and using the fact that w|_t=T'=0 and ∂_t(u-u_h)|_t=0,x∈𝒬_h=0, results in

    1/2(ρ e_h,e_h)_𝒬_h|_t=T' = 1/2(ρ e_h,e_h)_𝒬_h|_t=0 -1/2a(w,w)|_t=0 - ∫_0^T' (ρ∂_tϵ_h,e_h)_𝒬_h dt.

From the boundedness of ρ and Lemma <ref> it follows that

    e_h_0    ≤ Cρ^1/2 e_h_𝒬_h,       a.e. t∈(0,T).

Because of the coercivity of a we have

    -1/2a(w,w)|_t=0   < 0.

From the Cauchy–Schwarz inequality, the bounds of ρ, Lemma <ref>, and Lemma <ref>, we can also obtain

    |(ρ∂_tϵ_h, e_h)_𝒬_h| ≤ Ch^min(p+1,k)∂_tu_min(p+3,k+2)e_h_0,

for almost every t∈(0,T). Finally, we can use Lemma <ref>, Lemma <ref> and the boundedness of ρ to obtain

    ρ^1/2 e_h|_t=0_𝒬_h≤ Ch^min(p+1,k)u_∞, min(p+3,k+2) .

By taking the supremum of (<ref>) for all T'∈(0,T) and using (<ref>), (<ref>), (<ref>), and (<ref>), we can obtain

    e_h_∞,0≤ Ch^min(p+1,k)(u_∞,min(p+3,k+2)  + T∂_tu_∞,min(p+3,k+2)).

Using (<ref>) and Lemma <ref> we obtain (<ref>).











§ SEVERAL NEW MASS-LUMPED TETRAHEDRAL ELEMENTS OF DEGREES TWO TO FOUR



In this section, we present several novel mass-lumped tetrahedral elements for degree p=2,3,4. The new degree-2 and degree-3 elements use 15 and 32 nodes per element, respectively, while the current elements for these degrees require 23 and 50 nodes, respectively <cit.>. We also introduce several degree-4 elements, requiring 60, 61, and 65 nodes. Mass-lumped tetrahedral elements of degree four had not been found yet.







We present the mass-lumped tetrahedral elements using the reference tetrahedron with vertices at (0,0,0), (1,0,0), (0,1,0), and (0,0,1). In previous sections we used a tilde to denote coordinates and sets in the reference space, but since we only consider the reference space in this section, we will drop the tilde for readability. 

The nodes on the reference element are described using the notation {}, which denotes the node  and all equivalent nodes s(), with s∈𝒮. As shown in Lemma <ref>, any s∈𝒮 can be represented by a permutation of the barycentric coordinates. In this case, the barycentric coordinates are given by the three Cartesian coordinates x_1, x_2, x_3, and the additional coordinate x_4:=1-x_1-x_2-x_3, so any s∈𝒮 can be written as s(x_1,x_2,x_3)=(x_j,x_j,x_k), with i,j,k∈{1,2,3,4}, i≠ j, i≠ k, j≠ k. The barycentric coordinates of the node =(1/5,1/5,1/5), for example, are therefore given by (1/5,1/5,1/5,2/5), and the set of equivalent nodes {} consists of (1/5,1/5,1/5), (2/5,1/5,1/5), (1/5,2/5,1/5), and (1/5,1/5,2/5).

The reference function space, denoted by U, is the span of all nodal basis functions and is described in terms of {w}, which denotes the span of function w and all its equivalent functions w∘ s, with s∈𝒮. For example, all equivalent functions of w=x_1x_2 are x_1x_2, x_1x_3, x_1x_4, x_2x_3, x_2x_4, and x_3x_4, so {w} is the span of these six functions.

We assign the same weight to each equivalent node, so ω_=ω_s() for all s∈𝒮. From this and properties (<ref>) and (<ref>) it follows that if the quadrature rule is exact for a function w, then it is exact for all equivalent functions in {w}. If we can describe a function space in the form of {w_1,w_2,…,w_N}, by which we mean the span of w_1,w_2,…,w_N and all their equivalent versions, this means the quadrature rule is exact when it is exact for the N functions w_1,w_2,…,w_N. 

To give an example, the degree-3 element, given in Table <ref>, consists of the nodes (0,0,0), (a,0,0), (b,b,0), (c,c,c), and all equivalent nodes, and the function space for this element is given by U=𝒫_3⊕ℬ_f𝒫_1⊕ℬ_e𝒫_1, where ℬ_f:={β_f}:={x_1x_2x_3} are the face bubble functions and ℬ_e:={β_e}:={x_1x_2x_3x_4} is the internal bubble function and where we used the notation ℬ_f𝒫_k:=ℬ_f⊗𝒫_k, ℬ_e𝒫_k:=ℬ_e⊗𝒫_k. The quadrature rule should be exact for all functions in U⊗𝒫_1, which can be written as

    U⊗𝒫_1 = {x_1, x_1^2x_2, x_1^2x_2^2, β_fx_1,β_fx_1x_2, β_ex_1,β_ex_1x_2},

so as the span of 7 independent functions and all their equivalents. This means the quadrature rule should be exact for these 7 functions. Since this quadrature rule also has 7 parameters, namely 4 weights and three position parameters a,b,c, this results in a system of 7 equations with 7 unknowns. Solving this system results in the parameters given in Table <ref>. 

This approach has also been used to obtain the other elements presented in this paper. We have not yet found a systematic way to determine a suitable function space U⊃𝒫_p with a suitable configuration of the nodes. Instead, we just tried multiple configurations and checked if the resulting weights are all positive and the resulting nodes all lie on the reference triangle.

The degree-2 element with 15 nodes, the degree-3 element with 32 nodes, and the degree-4 element with 65 nodes are given in Tables <ref>, <ref>, and <ref>, respectively. In these tables, n denotes the number of nodes in the given equivalence class. Variants of the degree-4 element, requiring only 60 and 61 nodes, are given in Section <ref>.

In the next sections we test these new mass-lumped elements and compare them with the current mass-lumped elements and several discontinuous Galerkin approximations.





§ DISPERSION ANALYSIS



In this section we analyze the dispersion properties of the mass-lumped elements. The dispersion error is measured by the difference between the propagation speed of physical and numerical waves and is one of the main criteria to judge the quality of the finite elements for wave propagation modelling. We will use it to obtain an indication of the required mesh resolution for a given accuracy, and to compare different finite element methods in terms of accuracy and numerical cost.

For the analysis we will follow the same procedure as in <cit.>. We consider a homogeneous medium with ρ, c=1 and consider physical plane waves of the form 

    u=e^(·-ω t),

where :=√(-1) is the imaginary number,  is the wave vector, and ω is the angular velocity. Since ρ,c=1 we have a wave propagation speed c_P=1. For a given wave vector  we compute all corresponding numerical plane waves and determine the numerical wave with a propagation speed c_P,h=ω_h/|| closest to the physical wave velocity. The dispersion error is defined as the relative difference (c_P-c_P,h)/c_P. We then find the worst case among all possible wave directions for a fixed wave length λ=2π/||. We determine the dispersion error for different wavelengths and extrapolate the results to obtain a relation between the dispersion error and number of elements per wave length.



To obtain the numerical plane waves we construct a periodic tetrahedral mesh by packing a single parallelepiped cell with tetrahedra, and then repeating this pattern to fill the entire 3D-space. An illustration of such a mesh is given in Figure <ref>. Such a periodic mesh enables us to compute the numerical plane waves using Fourier modes and by solving an eigenvalue problem related to a single cell. 

To do this, let Ω_0 be the parallelepiped cell at the origin. We can write Ω_0:=T·[0,1)^3={y | y=T·, for some ∈[0,1)^3}, with T∈ℝ^3× 3 the second-order tensor whose columns are the vectors of the edges of Ω_0 connected to the origin. Let {^(Ω_0,i)}_i=1^n_0 be the set of nodes on Ω_0. For each k∈ℤ^3, we define the translated cell Ω_k:=T·k+Ω_0, and let {^(Ω_k,i)}_i=1^n_0 be the corresponding translated nodes. Then, for each node ^(Ω_k,i), we define w^(Ω_k,i) to be the corresponding nodal basis function. We can then define the following submatrices:

    M^(Ω_0)_ij   := (ρ w^(Ω_0,i),w^(Ω_0,j))_𝒬_h,       i,j=1,…,n_0, 
    
    A^(Ω_0,Ω_k)_ij   := a(w^(Ω_0,i),w^(Ω_k,j)),       k∈{-1,0,1}^3,  i,j=1,…,n_0.

For each wave vector  we then define the matrix

    S^()   := M_inv^(Ω_0)(∑_k∈{-1,0,1}^3 e^(·T·k)A^(Ω_0,Ω_k))

where M_inv^(Ω_0) denotes the inverse of M^(Ω_0). For an order-2K Dablain scheme, with time step size Δ t, the angular frequencies of the numerical plane waves {ω_h^(,i)}_i=1^n_0 are given by

    ω_h^(,i)   = ±1/Δ tarccos(∑_k=0^K1/(2k)!(-Δ t^2s_h^(,i))^k) ,

where {s_h^(,i)}_i=1^n_0 are the eigenvalues of σ(S^()) <cit.>. The numerical wave propagation speed is given by c_P,h^(,i)=|ω_h^(,i)|/||. The dispersion error, for a given wavelength λ, is then given by

    e_disp(λ) := sup_∈ℝ^3, ||=2π/λ( inf_i=1,…,n_0|c_P,h^(,i)-c_P|/c_P).


For our dispersion analysis, we will consider a congruent, nearly-regular, equifacial mesh, known as the tetragonal disphenoid honeycomb. This mesh can be obtained by a repeated pattern of cells, where a single cell can be obtained by slicing the unit cube into six tetrahedra with the planes x_1=x_2, x_1=x_3, and x_2=x_3, and then applying the transformation →T·, with 

    T := [       1    -1/3    -1/3;       0  √(8/9) -√(2/9);       0       0  √(2/3) ].
 

We will analyze the relation between the dispersion error and the number of elements per wavelength N_E:=(λ^3/|e|_av)^1/3, where |e|_av=2√(3)/27 denotes the average element volume. We will also look at the following quantities:

  
  * n_vec=n_0λ^3/|Ω_0|: the number of degrees of freedom per λ^3-volume. Here |Ω_0|=4√(3)/9 denotes the volume of Ω_0.
  
  * n_mat=λ^3/|Ω_0|∑_q∈𝒬_Ω_0 |𝒩(q)|: the number of non-zero entries of the stiffness matrix per λ^3-volume. Here, 𝒬_Ω_0 denotes the nodes on Ω_0 and |𝒩(q)| denotes the number of nodes connected with q through an element.
  
  * N_Δ t=T_0/Δ t: the number of time steps during one oscillation in time. Here T_0=λ/c_P denotes the duration of one oscillation and Δ t=√(c_K/s_h,max) is the largest allowed time step size for the order-2K Dablain scheme, with c_K a constant depending on the order of the time integration scheme (c_K=4,12,7.57,21.48 for K=1,2,3,4, respectively) and 
  
    s_h,max := sup_k∈𝒦_0max_i=1,…,n_0s^(,i)_h

the largest possible eigenvalue s_h, with 𝒦_0=T^-t·[0,2π) the space of distinct wave vectors. 
  
  * n_comp=n_matKN_Δ t: the estimated number of computations per λ^3-volume during one time oscillation, with K the number of stages of the order-2K Dablain scheme.

Details on the dispersion analysis and how the quantities listed above are computed can be found in <cit.>. 









We will refer to the standard linear mass-lumped finite element method as ML1. The higher-order mass-lumped methods will be referred to as ML[p]n[n], where p is the degree and n the number of nodes per element. In particular, the degree-2 method of <cit.> will be referred to as ML2n23 and the two versions of the degree-3 method in <cit.> will be referred to as ML3n50a and ML3n50b. The mass-lumped methods introduced in this paper will be referred to as ML2n15, ML3n32, ML4n60, ML4n61, and ML4n65.

We will also compare the mass-lumped methods with the symmetric interior penalty discontinuous Galerkin (SIPDG) method, introduced and analyzed in [Grote et al. 2006]. For the penalty term, we use the lower bound derived in <cit.>, since it was shown in <cit.> that this penalty term results in a significantly more efficient scheme than the penalty terms based on the more commonly used trace inequality of <cit.>. The quantities for the computational cost are computed in the same way as for the mass-lumped method, but now n_0 denotes the number of basis functions in Ω_0 and n_mat is computed as n_mat=λ^3/|Ω_0|∑_e∈𝒯_Ω_0 |U_e|^2|𝒩(e)|, where 𝒯_Ω_0 are the elements in Ω_0, |U_e| are the number of basis functions per element, and |𝒩(e)|=5 are the number of elements connected with e through a face, including e itself. We will refer to the SIPDG methods of degree 1, 2, 3, and 4, as DG1, DG2, DG3, and DG4, respectively. 

For the time integration, we combine each degree-p finite element method with an order-2p Dablain time integration scheme, since this results in order-2p convergence of the dispersion error. 
 
Figure <ref> illustrates the relation between the dispersion error and number of elements per wavelength. The dispersion error of the finite element methods converge with order 2p, which is typical for symmetric finite element methods for eigenvalue approximations, see, for example, <cit.> and the references therein. Using extrapolation, we obtain formulas for the dispersion error, given in Table <ref>. These formulas can be used to determine the required resolution of the mesh given the wavelength and desired accuracy. From the leading constants we can see that the new mass-lumped methods of degree 2 and 3 are more accurate for the same mesh resolution than the SIPDG and existing mass-lumped methods of these orders. The degree-4 mass-lumped method with 65 nodes is slightly more accurate than the versions using 60 and 61 nodes, but is slightly less accurate than the degree-4 discontinuous Galerkin method.

While some methods are more accurate for the same mesh resolution, this does not necessarily mean that these methods are more efficient, since the computational cost per element can greatly differ per method. To get an idea which method is most efficient for a given accuracy, we also look at the relation between the dispersion error and the estimated computational cost. This relation is illustrated in Figure <ref>. The required computational cost of each method for a dispersion error of 0.001 is also illustrated in Table <ref>. 

These results show that the new degree-2 mass-lumped method significantly outperforms the other degree-2 finite element methods, reducing the required computational cost for a given accuracy by one order of magnitude. The new degree-3 method is also significantly more efficient than the other degree-3 methods, reducing the required computational cost by more than a factor 5. These reductions in computational cost can be explained by the improved accuracy for the same mesh resolution, a reduction in the number of degrees of freedom and therefore reduction of the size of the stiffness matrix, and by a smaller number of time steps due to a larger allowed time step size. 

Among the degree-4 finite element methods, the mass-lumped method using 65 nodes performs best, mainly due to a smaller number of required time steps, although these differences are relatively small.

Figure <ref> also indicates that for a dispersion error between 1% and 0.1%, the new degree-2 mass lumped method performs best, while for smaller dispersion errors, the new degree-3 mass lumped method is most efficient. When we extrapolate the graphs, we find that the degree-4 mass-lumped method using 65 nodes will only outperform the degree-3 method for a dispersion error below 10^-5.

While the dispersion analysis provides useful information on the efficiency of the numerical methods, it does not include the effect of interpolation errors or inaccurate higher-frequency modes that may contaminate the numerical solution. Furthermore, the estimated computational cost is no perfect measure for the computation time, since the real computation time heavily depends on the implementation of the algorithm and the hardware that is used. In the next section we therefore also show the results of several numerical tests for the mass-lumped methods.







§ NUMERICAL TESTS





 §.§ Homogeneous Domain


We first tested and compared the old and new mass-lumped tetrahedral element methods on a homogeneous acoustic model using unstructured tetrahedral meshes. The domain is [-2,2]×[-1,1]×[0,2]km^3 and the acoustic wave propagation speed is c_P:=2km/s. A 3.5-Hz Ricker wavelet, starting from the peak, was placed at  _src:=(0,0,1000)m, and 56 receivers were placed on a line between x_r=-1375 and x_r=+1375m with a 50-m interval at y_r=0m and z_r=800m. Data were recorded for 0.6s, counting from the time at which the wavelet peaked, but the computations already started at the negative time -0.6s when the wavelet is approximately zero.






The exact solution, in case of an unbounded domain, is given by

    u(,t)    = w(t-r/c_P)/4π r,        t≤ 0.6 s,

where r:=|-_src| is the distance to the source and w(t):= (1-2π^2f^2t^2)e^-π^2f^2t^2 the Ricker-wavelet of peak frequency f=3.5Hz. This wavelet is zero up to machine precision for t≤-0.6s. For the bounded domain, on which we imposed zero Neumann boundary conditions, we add mirror sources to handle the reflections caused by the boundary conditions.


For the implementation of the mass-lumped methods we used the algorithm described in <cit.>. The time step size is based on the estimates of <cit.> multiplied by a factor 0.9. The simulations were carried out with OpenMP on 24 cores of two Intel® Xeon® E5-2680 v3 CPUs running at 2.50GHz. Figure <ref> shows the observed root mean square (RMS) errors of the receiver data for the various schemes against the number of degrees of freedom N and against wall clock time. The latter should not be taken too literal because it depends on code implementation, compiler and hardware, and even varies between runs. It can be further reduced by going to single precision, but then it becomes more difficult to measure the errors when they become small. Therefore, we ran a double-precision version of the code when preparing these figures.

Fourth-order time stepping was used for degrees higher than one <cit.>. For degree 4, we also considered 6th-order time stepping, but the errors were nearly the same as with 4th-order time stepping for the current example.

Power-law fits, given in Table <ref>, show that the RMS errors converge with approximately order p+1 and confirm that the new elements maintain an optimal order of convergence. Figure <ref> also shows that the new mass-lumped methods require less degrees of freedom and computation time for the same accuracy. For the degree-2 methods, the difference in wall clock time is up to one order of magnitude, while for the degree-3 methods this difference is up to a factor 2. The degree-4 methods become more efficient for errors below 10^-3.




 §.§ Elastic Salt Model


We also tested the methods on the more realistic salt model from <cit.>, made elastic by replacing the water layer at the top by rock. A 3-Hz Ricker wavelet vertical force source was placed on the surface at (2000,2200,0)m and 25 receivers were placed on a line between x_r=-1012.5 and x_r=7887.5m with a 25-m interval at y_r=2200m and z_r=0m. An illustration of this salt model is given in Figure <ref>. Figure <ref> displays vertical cross sections through the 3D vertical displacement wavefield of the 65-node degree-4 method, clipped at 25% of its maximum amplitude with red for positive and blue for negative values. Small amplitudes were replaced by the P-velocity to give an impression of the model. Figure <ref> also shows seismograms of this method for the displacement in the x- and z-direction clipped at 2% of the maximum amplitude. 










Simulations were carried out with the same implementation and in the same environment as for the homogeneous test case. The RMS errors are estimated by taking the traces for the 65-node degree-4 method as the `exact' solution. For the RMS error we use the data of receivers between x=2.1 and x=4km in the time interval [0,1.8]s. We selected this subset to exclude errors caused by the absorbing boundary layers. To compute the relative RMS error, we divide by the RMS of the data. 

An overview of the RMS errors and the wall-clock time is given in Table <ref>. The differences in the traces of the different degree-4 methods is of order 10^-5 and is much smaller then the estimated errors of the lower-degree elements. This indicates that the RMS errors of the degree-4 methods are of order 10^-5 and supports the idea that the degree-4 method can be used to estimate the accuracy of the lower-degree methods in this case. The Table illustrates again that the new degree-2 and degree-3 mass-lumped tetrahedral elements are more efficient than the current ones. The differences in computation time between the degree-4 methods is mainly due to the difference in number of time steps. The ML4n65 allows for a larger time step size than the other variants, which makes it slightly more efficient.







§ CONCLUSION


We developed a less restrictive accuracy condition for the construction of continuous mass-lumped elements, which enabled us to construct several new tetrahedral elements. The new degree-2 and degree-3 tetrahedral elements require 15 and 32 nodes, while the current versions require 23 and 50 nodes per element, respectively. These new elements require less degrees of freedom and allow larger time steps than the current versions. We also developed degree-4 tetrahedral elements with 60, 61, and 65 nodes per element. Mass-lumped tetrahedral elements of this degree had not been found yet. 

A dispersion analysis and numerical examples confirm that the new mass-lumped methods maintain an optimal order of accuracy and show that the new elements are significantly more efficient than the existing ones. In particular, the new degree-2 method is shown to be up to one order of magnitude faster than the current method, while the new degree-3 tetrahedral element results in a speed-up of up to a factor 2 for the same accuracy. The new degree-4 elements outperform the lower-degree elements for an accuracy below 10^-3. Among these degree-4 elements, the one with 65 nodes is  the most efficient, which is mainly due to a larger allowed time step size. The dispersion analysis also shows that the new degree-2 and degree-3 mass-lumped methods require significantly less degrees of freedom and number of time steps than the symmetric interior penalty discontinuous Galerkin methods of the same degree.

We have only considered tetrahedral elements in this paper, but the accuracy condition might also lead to more efficient triangular or higher-dimensional simplicial elements. Furthermore, although we focused only on linear wave propagation problems, mass lumping is useful for solving any type of evolution problem that requires explicit time-stepping.




abbrv




§ NODAL BASIS FUNCTIONS




Let e be a d-simplex in ℝ^d, with d≥ 0, and let s: ℝ^d→ℝ^d be an affine mapping that maps e onto itself. Then s can be represented by a permutation of the barycentric coordinates of e. In particular, there exists a permutation function P:ℝ^d+1→ℝ^d+1 such that 

    s()^*    = P(^*)       for all ∈ℝ^d,

where ^*, s()^*∈ℝ^d+1 denote the barycentric coordinates of  and s(), respectively.



Note that both → s()^* and → P(^*) are affine mappings from ℝ^d to ℝ^d+1. It therefore suffices to show that (<ref>) holds for all vertices of e. To do this, let v_i∈ℝ^d, for i=1,…,d+1, denote the vertices of e, and let p be the permutation of (1,…,d+1) such that s(v_i)=v_p(i). If we define P such that P(e_i)=e_p(i), with e_i∈ℝ^d+1 the unit vector in direction i, then

    s(v_i)^*    = v_p(i)^* 
       = e_p(i)
       = P(e_i) 
       = P(v_i^*)

for all i=1,…,d+1. 





Let 𝒬_h be defined as in Section <ref>, and let ∈e, for some ∈𝒬_h, e∈𝒯_h. If (<ref>) is satisfied, then ϕ_e^-1()∈𝒬̃.


By definition of 𝒬_h, there exists an element e^*∈𝒯_h and node ^*∈𝒬̃ such that =ϕ_e^*(^*), so ϕ_e^*^-1()=^*∈𝒬̃. Now construct a reference-to-reference element mapping s∈𝒮 such that s|_ϕ_e^*^-1(e^*∩e) = ϕ_e^-1∘ϕ_e^*|_ϕ_e^*^-1(e^*∩e). Then, using (<ref>), we can obtain ϕ_e^-1() = s∘ϕ_e^*^-1() = s(^*) ∈𝒬̃.




Let 𝒬_h be defined as in Section <ref>, and let ∈𝒬_h. If the conditions in (<ref>) and (<ref>) are satisfied, then the basis function w_, given in (<ref>), is well-defined and continuous and satisfies

    w_()   =δ_      for all ∈𝒬_h,

where δ_ denotes the Kronecker delta function.


The fact that w_ is well defined follows immediately from Lemma <ref>.

To prove that w_ is continuous, let f=∂ e^-∩∂ e^+ be any face adjacent to the elements e^-,e^+∈𝒬_h. It is sufficient to show that w_|_∂ e^+∩ f = w_|_∂ e^-∩ f, where w|_∂ e∩ f denotes the trace of function w restricted to e on face f. Suppose that ∉e^-∪e^+. Then w_|_∂ e^-∩ f = 0 = w_|_∂ e^+∩ f.

Now suppose that ∈e^+∖ f. Then ∉e^- and ^+∉f̃^+, where ^+:=ϕ_e^+^-1() and f̃^+:=ϕ_e^+^-1(f). From (<ref>) it then follows that w_|_∂ e^+∩ f = w̃_^+∘ϕ_e^+^-1|_f = 0 = w_|_∂ e^-∩ f.

Finally, suppose that ∈ f. Let s∈𝒮 be a reference-to-reference element mapping such that s|_f̃^+ = ϕ_e^-^-1∘ϕ_e^+|_f̃^+. We can then derive

    w_|_∂ e^+∩ f   = w̃_^+∘ϕ_e^+^-1|_f 
       = w̃_s(^+)∘ s ∘ϕ_e^+^-1|_f 
       = w̃_ϕ_e^-^-1()∘ϕ_e^-^-1|_f 
       = w_|_∂ e^-∩ f,

where the second line follows from (<ref>).

Since w_ is continuous, w_() is well-defined for any ∈Ω. To prove property (<ref>), suppose there is an element e∈𝒯_h such that ,∈e. Define :=ϕ_e^-1() and :=ϕ_e^-1(). Then w_()=w̃_() = δ_=δ_.

Now suppose that there is no element e such that ,∈e. Then ≠ and there exists an element e such that ∈e and ∉e. By definition of w_ it then follows that w_()=0=δ_.





§ VARIANTS OF THE DEGREE-FOUR MASS-LUMPED TETRAHEDRAL ELEMENT









