Efficient quadrature rules for computing the stiffness matrices of mass-lumped tetrahedral elements for linear wave problems
S. Geevers, W. A. Mulder, J. J. W. van der Vegt
TL;DR
This paper introduces new efficient quadrature rules for computing stiffness matrices in mass-lumped tetrahedral finite elements, enabling faster wave propagation simulations with heterogeneous materials while maintaining accuracy.
Contribution
The paper develops novel quadrature rules for higher-degree mass-lumped tetrahedra that reduce computational cost and support heterogeneous materials without losing convergence rate.
Findings
New quadrature rules for degrees 3 and 4 tetrahedra with fewer points
Numerical results show improved efficiency over exact stiffness matrix computation
Optimal convergence order maintained with heterogeneous materials
Abstract
We present new and efficient quadrature rules for computing the stiffness matrices of mass-lumped tetrahedral elements for wave propagation modelling. These quadrature rules allow for a more efficient implementation of the mass-lumped finite element method and can handle materials that are heterogeneous within the element without loss of the convergence rate. The quadrature rules are designed for the specific function spaces of recently developed mass-lumped tetrahedra, which consist of standard polynomial function spaces enriched with higher-degree bubble functions. For the degree-2 mass-lumped tetrahedron, the most efficient quadrature rule seems to be an existing 14-point quadrature rule, but for tetrahedra of degrees 3 and 4, we construct new quadrature rules that require less integration points than those currently available in the literature. Several numerical examples confirm…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9| 2 | ||
| 3 | ||
| 4 | ||
| Type | Points | Description | |
|---|---|---|---|
| 1 | centre of tetrahedron | ||
| 4 | on line through vertex and centre | ||
| 6 | on line through edge-midpoint and centre | ||
| 12 | on plane through centre and two vertices | ||
| 24 | arbitrary position |
| Nodes | node parameters | ||
|---|---|---|---|
| Nodes | node parameters | ||
|---|---|---|---|
| - | |||
| Nodes | node parameters | ||
|---|---|---|---|
| - | |||
| Nodes | node parameters | ||
|---|---|---|---|
| method | method | method | |||
|---|---|---|---|---|---|
| 2n15 | 4n60 | 4n65 | |||
| 2n15q14 | 4n60q51 | 4n65q60 | |||
| 2n15q15 | 4n60q59 | 4n65q79 | |||
| 3n32 | 4n61 | ||||
| 3n32q21 | 4n61q60 | ||||
| 3n32q24 | 4n61q79 |
| method | method | method | |||
|---|---|---|---|---|---|
| 2n15 | 4n60 | 4n65 | |||
| 2n15q14 | 4n60q51 | 4n65q60 | |||
| 2n15q15 | 4n60q59 | 4n65q79 | |||
| 3n32 | 4n61 | ||||
| 3n32q21 | 4n61q60 | ||||
| 3n32q24 | 4n61q79 |
| time (s) | RMS error | |||||
|---|---|---|---|---|---|---|
| 2 | 15 | - | 366 | 988 | ||
| 14 | 371 | 750 | ||||
| 3 | 32 | - | 338 | 307 | ||
| 21 | 327 | 172 | ||||
| 24 | 336 | 204 | ||||
| 4 | 60 | - | 744 | 725 | ||
| 51 | 697 | 366 | ||||
| 59 | 723 | 407 | ||||
| 4 | 61 | - | 631 | 670 | ||
| 60 | 614 | 365 | ||||
| 79 | 614 | 471 | ||||
| 4 | 65 | - | 568 | 621 | ||
| 60 | 553 | 303 | ||||
| 79 | 553 | 388 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Efficient quadrature rules for computing the stiffness matrices of mass-lumped tetrahedral elements for linear wave problems*
S. Geevers1, W.A. Mulder2,3 and J.J.W. van der Vegt1
-
Department of Applied Mathematics, University of Twente, Enschede, the Netherlands
-
Shell Global Solutions International BV
-
Delft University of Technology
Abstract.
We present new and efficient quadrature rules for computing the stiffness matrices of mass-lumped tetrahedral elements for wave propagation modelling. These quadrature rules allow for a more efficient implementation of the mass-lumped finite element method and can handle materials that are heterogeneous within the element without loss of the convergence rate. The quadrature rules are designed for the specific function spaces of recently developed mass-lumped tetrahedra, which consist of standard polynomial function spaces enriched with higher-degree bubble functions. For the degree-2 mass-lumped tetrahedron, the most efficient quadrature rule seems to be an existing 14-point quadrature rule, but for tetrahedra of degrees 3 and 4, we construct new quadrature rules that require less integration points than those currently available in the literature. Several numerical examples confirm that this approach is more efficient than computing the stiffness matrix exactly and that an optimal order of convergence is maintained, even when material properties vary within the element.
*This work was funded by the Shell Global Solutions International B.V. under contract no. PT45999.
1. Introduction
Mass-lumped tetrahedral element methods are efficient methods for solving linear wave equations, such as the acoustic wave equation, the elastic wave equations, or Maxwell’s equations, on complex 3D domains with sharp material interfaces [25]. They offer the same convergence rate and geometric flexibility as standard continuous tetrahedral element methods, but also allow for explicit time-stepping due to a diagonal mass matrix.
To obtain mass-lumped elements, Lagrangian basis functions are combined with an inexact quadrature rule for computing the mass matrix. A lumped matrix is obtained when the quadrature points coincide with the basis function nodes. For quadrilateral and hexahedral elements, mass-lumping is achieved using tensor-product basis functions and Gauss–Lobatto points, resulting in the well-known spectral element method [20, 21, 15]. For linear triangular and tetrahedral elements, mass-lumping is achieved with standard Lagrangian basis functions and a Newton–Cotes integration rule. For higher-degree triangular and tetrahedral elements, however, the element space needs to be enriched with higher-degree bubble functions in order to maintain stability and an optimal order of convergence [6]. So far, mass-lumped triangular elements of degree 2 and 3 [5, 7], 4 [17], 5 [3], 6 [18], and 7-9 [16, 9] have been found. The first higher-degree tetrahedral elements were presented in [17] for degree 2 and in [3] for degree 3. Recently, we presented new mass-lumped tetrahedral elements of degrees 2 to 4 in [12]. The new degree-2 and degree-3 elements require significantly less nodes than the earlier versions, while mass-lumped tetrahedra of degree 4 had not been found before. Because of the reduced number of nodes, these new mass-lumped elements are also much more efficient than the earlier versions [12] and are therefore more suitable for large-scale 3D simulations.
A question that remains is how to efficiently compute the stiffness matrix for these elements. When the material parameters are piecewise constant, the stiffness matrix can be evaluated exactly [19]. Alternatively, we can use a quadrature rule to approximate the stiffness matrix. This latter approach can significantly reduce the number of computations as we will demonstrate in Section 6. Moreover, it also allows us to handle material parameters that vary within the element without loss of convergence rate as we will prove in Section 4.
Finding an efficient quadrature rule for the stiffness matrix for mass-lumped tetrahedra is not straightforward. For hexahedral elements, the stiffness matrix can be approximated with the same Gauss–Lobatto quadrature rule that is used for the mass matrix, but for mass-lumped tetrahedra, using the quadrature rule of the mass matrix to also evaluate the stiffness matrix turns out to be inefficient or inaccurate. In this paper, we therefore present new and efficient quadrature rules for computing the stiffness matrices for mass-lumped tetrahedra.
To obtain these quadrature rules, we show that the quadrature rule only needs to be exact for functions in the space , where denotes the degree of the element, denotes the space of polynomials up to degree , and denotes the space of partial derivatives of functions in the element space. Since the mass-lumped tetrahedra contain higher-degree bubble functions, so does the space that needs to be integrated exactly. Most quadrature rules in literature, however, are designed to be exact for spaces of the form . Such quadrature rules for tetrahedral domains can be found in, for example, [22, 13, 14, 8, 24, 23] and the references therein. We could choose equal to the highest polynomial degree that appears in , but the resulting number of quadrature points may then be suboptimal. Instead, we try to find quadrature rules that are exact for with a minimal number of quadrature points.
For the degree-2 tetrahedral element, the most efficient quadrature rule still seems to be the 14-point rule of [13] that is accurate for polynomials up to degree 5. For the degree-3 element and the three degree-4 elements of 60, 61, and 65 nodes, however, we present new quadrature rules that require 21, 51, 60, and 60 points, respectively, while using a quadrature rule from the current literature would require 24 [14], 59 [24], 79 [24], and 79 [24] points.
This paper is organised as follows. In Section 2, we introduce the mass-lumped finite element method, and in Section 3, we present our new quadrature rules for evaluating the stiffness matrix. We prove in Section 4 that the conditions used to obtain our quadrature rules result in optimal convergence rates. In Section 5, we analyze and compare the dispersion properties and resulting time step size for our new quadrature rules and several other rules available in the literature. In Section 6, we show numerical examples demonstrating that using our numerical quadrature rules for evaluating the stiffness matrix is more efficient than evaluating the integrals exactly and that the convergence rate is not lost when material parameters vary within the elements. Finally, we summarise our main conclusions in Section 7.
2. The mass-lumped finite element method
To present and analyze the mass-lumped finite element method, we consider the scalar wave equation given by
[TABLE]
where is the spatial domain, is the time domain, is the scalar field that needs to be solved, is the gradient operator, is the source term, are the initial values, and are positive spatial parameters. The spatial domain is assumed to be a bounded open domain with Lipschitz boundary , and the parameters and are assumed to be bounded by and , with strictly positive constants.
To solve the scalar wave equation with a finite element method, we consider the weak formulation. Let denote the standard Lesbesque space of square-integrable functions on , the standard Sobolev space of functions in that vanish on and have square-integrable weak derivatives, and , with a Banach space, the Bochner space of functions such that is square integrable on . Assume , , and . The weak formulation of (1) can then be written as finding u\in L^{2}\big{(}0,T;H^{1}_{0}(\Omega)\big{)}, with \partial_{t}u\in L^{2}\big{(}0,T;L^{2}(\Omega)\big{)} and \partial_{t}(\rho\partial_{t}u)\in L^{2}\big{(}0,T;H^{-1}(\Omega)\big{)}, such that , , and
[TABLE]
where denotes the standard inner-product and denotes the pairing between and .
This weak form of the wave equation can be solved with the mass-lumped finite element method, which consists of the following components:
- •
a tetrahedral mesh , where denotes the radius of the smallest sphere that can contain each element,
- •
a reference tetrahedron with reference space , where denotes the space of polynomials of degree or less and a space of higher-degree face and interior bubble functions,
- •
a set of reference nodes that can be used for both interpolation and quadrature on ,
- •
a set of quadrature weights .
Using these components, a finite element space can be constructed of the form
[TABLE]
where
[TABLE]
with the reference-to-physical element mapping. The interpolation points are given by , where
[TABLE]
and the inner-product is approximated by
[TABLE]
with and the volume of and , respectively, and , .
Assume and are all continuous. The finite element method can then be formulated as finding such that , , and
[TABLE]
where is the interpolation operator that interpolates a continuous function at the points by a function in .
Now let be the set of all interpolation points that do not lie on the boundary , and define nodal basis functions such that for all , with the Kronecker delta. Also define, for any continuous function , the interpolation vector such that for all . The finite element method can then be formulated as finding such that , , and
[TABLE]
Here, , with , is the mass matrix, and , with , is the stiffness matrix.
Since the interpolation points and quadrature points coincide, the mass matrix is diagonal with entries . Therefore, we can efficiently solve the system of ODE’s in (4) using an explicit time-stepping scheme. Standard conforming finite element methods do not result in a (block)-diagonal mass matrix and are therefore less suitable for solving wave equations on large three-dimensional meshes.
To remain accurate and stable, the mass-lumped finite element method needs to satisfy the following conditions [12]:
- C1
(Unisolvent). The space is unisolvent on the nodes . 2. C2
(Symmetry). The space and the set are invariant to affine mappings that map onto itself. 3. C3
(Face-conforming). If is zero at all nodes in , with a reference face, then is zero on . 4. C4
(Positivity). The weights are all strictly positive. 5. C5
(Accuracy). The quadrature rule is exact for functions in when .
The first three conditions are necessary to guarantee that the global basis functions are well-defined and continuous. The last two conditions are necessary for stability and for maintaining an optimal order of convergence.
When , these conditions can not all be met for standard polynomial spaces . Therefore, the element space needs to be enriched with higher-degree bubble functions. We will focus on the mass-lumped tetrahedral elements recently presented in [12]. An overview of these elements is given in Table 1. There, denotes the dimension of , which is equal to the number of nodes per element, denotes the span of the four face bubble functions, and denotes the span of the element bubble function, with , , , the four barycentric coordinates. We also used the notation for any two function spaces .
To apply these elements more efficiently, we also approximate the inner-product for the stiffness matrix, , with a quadrature rule. This also allows us to handle material parameters that vary within the element. It turns out that it is more efficient and sometimes even necessary to compute the stiffness matrix with a different quadrature rule than for the mass matrix. We will denote the quadrature points and weights for the stiffness matrix by and , respectively, and denote the corresponding approximated -product by .
The resulting finite element method remains stable and accurate if the following conditions are also satisfied:
- C6
(Positivity). The weights are all strictly positive. 2. C7
(Spurious-free). There is no function with zero gradient on all quadrature points except the constant function. In case of linear elasticity, there is no function with zero strain on all quadrature points except the six rigid motions. 3. C8
(Accuracy). If , the quadrature rule for the stiffness matrix is exact for functions in , where denotes the space of all partial derivatives of all functions in .
A proof that these three conditions are sufficient to maintain an optimal order of convergence is given in Section 4.
We constructed quadrature rules that satisfy these conditions for the specific function spaces of the higher-degree mass-lumped tetrahedra presented in Table 1. For the degree-2 element, the most efficient quadrature rule seems to be an existing 14-point rule that is fifth-order accurate, but for the higher-degree elements, we obtained new quadrature rules that require less points than those currently available in the literature. An overview of these rules is given in the next section.
3. Efficient quadrature rules for the stiffness matrix
To present the quadrature rules for the stiffness matrix, let be the reference tetrahedron with vertices at , , , and . The barycentric coordinates of this element are given by the three Cartesian coordinates , , , and the fourth coordinate . These coordinates are useful for describing , the space of affine mappings that map onto itself, since any function can be defined by a permutation of the barycentric coordinates. In particular, we can write for some , , for any .
Now, let denote point and all equivalent points , with . The quadrature rule will consist of several equivalence classes with quadrature weights that are the same within each equivalence class. To give an example of an equivalence class, consider the point . The barycentric coordinates of this point are given by , so the equivalence class consists of the four points , , , and when . An overview of the different types of points is given in Table 2. The configuration of a quadrature rule is given by the numbers , , , , and , which indicate that the quadrature rule has distinct points of type , points of type , points of type , etc.
To find a set of points and weights that satisfy accuracy condition C8, we construct a linear basis that spans . We describe this linear basis and linear span using the notation , which denotes the span of the functions and all equivalent functions , with and . To give an example, all equivalent versions of are given by the six functions , , , , , and , so denotes the span of these six functions.
After having constructed a basis for , we search for a quadrature rule that has a configuration with parameters. These parameters consist of location parameters and quadrature weights. Because of the symmetry, a quadrature rule that is exact for a function is exact for all its equivalent functions. Therefore, to satisfy C8, we end up with a nonlinear system of equations:
[TABLE]
We obtain solutions of this system using Newton’s method for a large number of different initial values and check for each solution if it satisfies C6 and C7. When we cannot find an admissible solution for configurations with parameters, we increase and try again. We continue this process until we find a suitable quadrature rule.
The complete algorithm for finding a quadrature rule of parameters can be summarised as follows:
Construct basis functions , such that 2. 2.
Choose a configuration for the quadrature rule such that
[TABLE] 3. 3.
Solve the system of equations given in (5) using Newton’s method with a random initial guess. 4. 4.
If the method converges within a maximum number of iterations, check if the solution satisfies conditions C6 and C7. 5. 5.
Repeat Steps 3 and 4 until a maximum number of trials is reached or an admissible solution is found.
The quadrature rules that were obtained in this way are given in Tables 3-6. There, denotes the number of nodes in each equivalence class and and denote the face bubble function and interior bubble function, respectively.
The quadrature rule with the least number of points we could find for the degree-2 15-node tetrahedron is the 14-point fifth-order accurate rule of [13]. We also found an accurate quadrature rule of 10 points with positive weights, but the resulting method was not accurate since condition C7 was not satisfied: it had one non-constant mode with gradient equal to zero at all 10 points. We also considered the 15-point quadrature rule used for the mass matrix, which also satisfies C6-C8. However, this rule significantly increases the condition number of the element matrix and therefore results in a considerably smaller time step size as shown in Section 5.
The quadrature rules for the degree-3 and degree-4 elements are new and require less quadrature points than rules currently available in the literature, since most quadrature rules in the literature are constructed to be exact for a function space of the form and not for the specific function spaces . To give an example, the highest polynomial degree of of the degree-4 61- or 65-node tetrahedron is 7, so contains a polynomial of degree 10. A quadrature rule that is order-10 accurate already requires 79 quadrature points [24], while our quadrature rule for these elements only requires 60 points. Similarly, our quadrature rules for the degree-3 32-node tetrahedron and the degree-4 60-node tetrahedron require 21 and 51 points, respectively, while the quadrature rules currently available in the literature for these elements require 24 and 59 points [24].
4. Error estimates
In this section, we prove that, when conditions C1-C8 are satisfied, the finite element method maintains an optimal order of convergence for a related elliptic problem. Convergence for the wave equation can then be derived in a way analogous to [12, Chapter 4.6].
Throughout this section, we will let denote the degree of the finite element space, by which we mean the largest degree such that . We will also let denote a positive constant that may depend on the domain , the regularity of the mesh, the parameters and , the reference space , and the reference quadrature rules and , but does not depend on the mesh resolution and the functions that appear in the inequalities.
4.1. Preliminary results
To obtain error bounds, we first define some norms and function spaces and list a few preliminary results. Let denote the Sobolev space of functions on with order- square-integrable weak derivatives and equip the space with norm
[TABLE]
where denotes the -norm, the partial derivative, and the order of the derivative. We also define the broken Sobolev spaces , equipped with norm
[TABLE]
Throughout this section, we will use the fact that for any three-dimensional domain .
We also define the semi-norms and for piecewise continuous functions and , and define to be the -projection operators projecting onto the discontinuous piecewise-polynomial spaces . Several useful properties of these spaces and operators are listed below.
Lemma 4.1**.**
Let . Then
[TABLE]
Proof.
These results follow immediately from the fact that the elements are affine equivalent with the reference element and that the reference space is finite dimensional. ∎
Lemma 4.2**.**
If conditions C1-C4 are satisfied, then becomes a full norm on and
[TABLE]
Furthermore, if conditions C1-C3, C6, and C7 are satisfied, then becomes a full norm on and
[TABLE]
Proof.
These inequalities follow immediately from the fact that the elements are affine equivalent with the reference element and that the reference element space is finite dimensional. ∎
Lemma 4.3**.**
Let and , with , and let . Then
[TABLE]
Furthermore, if , then
[TABLE]
Finally, if , with , then
[TABLE]
with the degree of the finite element space.
Proof.
The first, second, and last inequality follow from [4, Chapter 3.1]. For the fourth inequality, let , be a polynomial degree and a set of points such that is unisolvent on . Also, let be the interpolation operator that interpolates a function in at the nodes by a function in . We can then obtain the fourth inequality as follows:
[TABLE]
where we used Lemma 4.1 in the the second line and the triangle inequality in the third line. The last line follows from [4, Chapter 3.1].
The third inequality can be derived in a way analogous to the fourth inequality. ∎
4.2. Estimates on the integration error
Define the two integration errors and . In [12] we derived the following bounds on .
Lemma 4.4**.**
Let be the degree of the finite element space, , with , and . If conditions C1-C5 are satisfied, then
[TABLE]
We also derive bounds on the integration error for the stiffness matrix.
Lemma 4.5**.**
Let be the degree of the finite element space, with , and . If conditions C1-C3, C6, and C8 are satisfied, then
[TABLE]
Proof.
Using C8, we can write
[TABLE]
Inequality (6) then follows from the Cauchy–Schwarz inequality and Lemma 4.3.
Using C8 and the fact that for , we can also write
[TABLE]
Inequality (7) then follows from the Cauchy–Schwarz inequality and Lemma 4.3. ∎
4.3. Error estimates for a related elliptic problem
Let . The elliptic problem corresponding to (2) is finding such that
[TABLE]
The corresponding mass-lumped finite element method is finding such that
[TABLE]
In the next two theorems we prove optimal convergence in the -norm and -norm.
Theorem 4.6** (Optimal convergence in the -norm).**
Let be the solution of (8) and the solution of (9). Assume , , and , with . If conditions C1-C8 are satisfied, then
[TABLE]
with the degree of the finite element method.
Proof.
Define and . Using (8), we can write
[TABLE]
and using (9), we can obtain
[TABLE]
Subtracting these two equalities gives
[TABLE]
From Lemma 4.2, the positivity of , and Poincaré’s inequality, it follows that
[TABLE]
Using Lemma 4.5, Lemma 4.3, and the regularity of , we can obtain
[TABLE]
Using the Cauchy–Schwarz inequality, the boundedness of , and Lemma 4.3, we can also obtain
[TABLE]
Finally, using Lemma 4.4, we can obtain
[TABLE]
[TABLE]
Since , inequality (10) then follows from the above and Lemma 4.3. ∎
To prove optimal convergence in the -norm, we make the following regularity assumption: for any , the solution of (8) is in and satisfies
[TABLE]
This is certainly true when is and .
Theorem 4.7** (Optimal convergence in the -norm).**
Let be the solution of (8) and the solution of (9). Assume , , and , with , and assume the regularity condition (16) holds. If conditions C1-C8 are satisfied, then
[TABLE]
with the degree of the finite element method.
Proof.
Define to be the solution of the elliptic problem
[TABLE]
From the regularity assumption it follows that and
[TABLE]
From the definition of , it also follows that
[TABLE]
We can bound the term as follows:
[TABLE]
where we used the Cauchy–Schwarz inequality and the boundedness of in the first line, Theorem 4.6 and Lemma 4.3 in the second line, and the regularity assumption in the last line. It then remains to find a bound for .
To do this, use (8) to write
[TABLE]
and use (9) to write
[TABLE]
Subtracting these two equalities gives
[TABLE]
Now, set . We can write
[TABLE]
where we used C8 for the first and third equality. We can bound as follows:
[TABLE]
where we used the Cauchy–Schwarz inequality for the first inequality, Lemma 4.1, the regularity of , and Lemma 4.3 for the second inequality, the triangle inequality and the regularity assumption for the third inequality, and Theorem 4.6 and Lemma 4.3 for the last inequality. We can also bound as follows:
[TABLE]
Here, the first line follows from Lemma 4.5 and the fact that , the second line follows from the regularity of , the third line follows from Lemma 4.3, the fifth line follows from Lemma 4.3, and the last line follows from the regularity assumption. The fourth line follows from the fact that is piecewise polynomial of degree and therefore . By combining the bounds on and , we then obtain
[TABLE]
From Lemma 4.4, Lemma 4.3, and the regularity assumption, it also follows that
[TABLE]
Combining (20), (21), and (4.3) gives
[TABLE]
Combining this with (4.3) and (4.3) then gives (17). ∎
These results can be used to prove optimal convergence for the wave equation in a way analogous to [12, Chapter 4.6] by replacing by and by defining the projection operator of [12] such that for all .
4.4. Error estimates for the linear elastic case
So far, we only analyzed the scalar wave equation, but we can obtain error estimates for the elastic wave equations in an analogous way.
In the linear elastic case, the wave field is a vector field and (1a) becomes
[TABLE]
with , where is the elastic tensor field with symmetries and bounds
[TABLE]
with , strictly positive constants and .
The only part of the error analysis that requires some additional work in this case, is the second inequality of Lemma 4.2. Instead of , we need to show that, if conditions C1-C3, C6, and C7, are satisfied, then
[TABLE]
This result follows from the fact that is finite dimensional and from the relations
[TABLE]
where is the Jacobian of the element mapping, denotes the transposed of , , and .
Using the boundedness of , (23), and Korn’s inequality, we can show that the bilinear operator for the elastic case is still coercive. The other parts of the error analysis are analogous to the scalar case.
5. Dispersion analysis
To test the effect of the new quadrature rules on the accuracy and time step size, we first analyze the dispersion properties of the resulting mass-lumped finite element method along the lines of [11]. We consider a homogeneous unbounded domain with and consider plane waves of the form
[TABLE]
Here, denotes the imaginary number, denotes the wave vector, and denotes the angular velocity. We also let denote the wavelength and denote the wave propagation speed. The angular velocity and wave propagation speed satisfy the relation .
To analyze the numerical dispersion, we consider a translation-invariant mesh constructed from a repeated cell pattern as illustrated in Figure 1 and derive the propagation speeds of the numerical plane waves. The dispersion error is defined as the error in the numerical wave propagation speed: . Since the mesh is constructed from a repeated cell pattern, obtaining the numerical wave propagation speed for a given wave vector requires solving an eigenvalue problem related to only a single cell.
To construct the translation-invariant mesh, we subdivide the unit cell into tetrahedra and repeat this pattern to pack the entire 3D space. We also apply a linear transformation , with and , to the mesh in order to obtain more regular tetrahedra.
Let denote all the nodes on , let denote the translated nodes on the translated cell , and let denote the corresponding nodal basis functions. We define matrices as follows:
[TABLE]
For any wave vector , we then define the matrix as follows:
[TABLE]
When using an order- Dablain time integration scheme [10], with time step size , the angular velocities of the numerical plane waves are given by
[TABLE]
where are the eigenvalues of [11]. The numerical wave propagation speeds are given by . The dispersion error, for a given wavelength , is then given by
[TABLE]
For the dispersion analysis, we consider a mesh of congruent nearly-regular isofacial tetrahedra, known as the tetragonal disphenoid honeycomb. This mesh is obtained by slicing the unit cube into 6 tetrahedra with the planes , , and , and applying a linear transformation , with
[TABLE]
An illustration of this mesh is given in Figure 1.
We plot the dispersion error for different elements and quadrature rules against the number of elements per wavelength , where denotes the average element volume. We also compute the largest allowed time step size, given by
[TABLE]
with a constant depending on the order of the time integration scheme ( for , respectively) and
[TABLE]
the largest eigenvalue of the discrete spatial operator, with the space of distinct wave vectors. Details on the dispersion analysis can be found in [11].
We test the degree- mass-lumped finite element methods presented in [12] and given in Table 1 using exact stiffness matrix evaluation and using the quadrature rules presented in this paper and quadrature rules that are accurate up to degree from [24], with the highest polynomial degree of the enriched element space. For the degree-2 method, we also test using the quadrature rule of the mass matrix for evaluating the stiffness matrix. We combine each method with an order- Dablain time integration scheme.
The dispersion error versus the number of elements per wavelength is shown in Figure 2 and extrapolations of these graphs are given in Table 7. The figure and table show that the dispersion error of all methods converges with order , which is typical for eigenvalue and dispersion errors of symmetry-conserving finite element approximations, see for example, [2] and the references therein. The figure and table also show that there is hardly any difference in accuracy between the methods of the same degree and that the methods using a quadrature rule to evaluate the stiffness matrix have a nearly the same or even slightly smaller dispersion error than the methods using exact integration.
The largest allowed time step size for each method is given in Table 8. The table shows that the quadrature rules for the stiffness matrix tested here hardly affect the largest allowed time step size, except for the degree-2 15-point mass matrix quadrature rule, which reduces the allowed time step size by more than a factor 1.5. For the other quadrature rules, the largest allowed time step size remains nearly the same or becomes even slightly larger.
If the stiffness matrix is evaluated with a quadrature rule and the resulting accuracy and time step size remain nearly the same, then the number of computations to obtain a given accuracy mainly depends on the number of quadrature points. Since the quadrature rules presented in this paper satisfy these properties and require less quadrature points than those currently available in the literature, they can result in a reduction of the computational cost proportionate to the reduction in number of quadrature points.
6. Numerical tests
6.1. Algorithms for computing the element stiffness matrices
Before we present the numerical tests, we first briefly describe the algorithms for computing the element stiffness matrices. In particular, we show how we efficiently compute the element stiffness matrix-vector products on the fly. We do not store the matrices, since this requires storing and fetching significantly more data, and since it was shown in [19] that an on-the-fly approach is more efficient for higher-degree elements.
To describe the algorithms, let be an arbitrary element. We introduce the following notation.
- •
: nodes on reference element . Nodes of the different mass-lumped elements can be found in [12].
- •
: nodal basis function corresponding to .
- •
: nodal basis function of the physical element.
- •
: quadrature points for the stiffness matrix on reference element . Quadrature rules for the different elements are given in Section 3.
- •
: quadrature weight corresponding to .
- •
: the element stiffness matrix.
- •
: the wave field on .
- •
: the wave field at the nodes on .
When using exact integration, the stiffness matrix-vector product is given by
[TABLE]
for . After rewriting the integral as an integral over the reference element, this becomes
[TABLE]
where , is the gradient operator in reference coordinates, and is a tensor field, with the Jacobian of the element mapping and the transposed of . When is constant within each element, then is also constant and we can compute (26) using the algorithm of [19]:
[TABLE]
where are precomputed matrices, given by
[TABLE]
for , , with the derivative in reference coordinate . We can reduce the number of computations in (27) using the fact that is symmetric:
[TABLE]
where if and when . The complete algorithm can then be described as follows:
- A1.
Compute for , :
[TABLE] 2. A2.
Compute :
[TABLE]
The computational work is dominated by the first step where 6 matrix-vector products with matrices of size are computed.
Alternatively, we can compute by evaluating the integrals with a quadrature rule. Equation (26) then becomes
[TABLE]
where . We can compute this as follows:
[TABLE]
where are precomputed matrices, given by
[TABLE]
The complete algorithm can be described as follows:
- B1.
Compute for :
[TABLE] 2. B2.
Compute for :
[TABLE] 3. B3.
Compute :
[TABLE]
The computational work for this algorithm is dominated by the first and third step, which both require 3 matrix-vector products with matrices of size , so 6 of these matrix-vector products in total. Since for all the quadrature rules presented in this paper, this number of computations is smaller than for the previous algorithm, although only slightly. However, as we will show next, this quadrature-based algorithm is significantly more efficient than the exact-integral algorithm in case of linear elasticity. Moreover, this quadrature-based algorithm also works if varies within the element.
In case of linear elasticity, the wave field is a vector field and the term becomes the stress tensor , with the order-four elasticity tensor with symmetries and . The vector can in this case be written as a concatenation of three vectors , where is the wave field component at the nodes on . The parameter becomes , where , and becomes . The algorithm for computing the element stiffness matrix-vector product using exact integration then becomes
- A1*.
Compute for :
[TABLE] 2. A2*.
Define . Compute for :
[TABLE]
The computational work is again dominated by the first step, which now requires 27 matrix-vector products with matrices of size .
When using a quadrature rule, the algorithm becomes
- B1*.
Compute for :
[TABLE] 2. B2*.
Compute for :
[TABLE] 3. B3*.
Define . Compute for :
[TABLE]
The computational work for this algorithm is dominated again by the first and third step, which now both require 9 matrix-vector products with matrices of size , so 18 of these matrix-vector products in total. The number of computations is therefore reduced by more than a factor 1.5 when compared to the algorithm based on exact integration. Furthermore, the quadrature-based algorithm can also handle tensor fields that vary within the element.
Both algorithms can be slightly improved by exploiting the fact that the rows and columns of the matrices and the columns of matrices sum to zero. Furthermore, in case of isotropic elasticity, steps A2* and B2* can be computed more efficiently by exploiting the simple structure of the elasticity tensor .
In the next subsections, we demonstrate the superiority of the quadrature-based algorithm for the case of non-constant parameters and linear elasticity.
6.2. Acoustic wave on a heterogeneous domain
We first test the methods for an acoustic wave propagation problem with a heterogeneous domain. The acoustic wave equation is given by
[TABLE]
with spatial domain , time interval , pressure field , mass density , and acoustic wave speed . We choose and impose zero Neumann boundary conditions.
To construct an analytic solution, let , for , be distorted coordinates, with and , and define . Also let be the average mass density, the average wave speed, the wave vector, and the angular velocity, and let parameters and be given by
[TABLE]
Then the standing wave, given by
[TABLE]
is a solution of (31) that satisfies the zero Neumann boundary conditions.
Now, set km, , , for , and km/s, g/cm3. To test the numerical methods, we use and as initial conditions. We test on multiple unstructured meshes and simulate in time using a fourth-order time-stepping scheme [10] with time step size , where is the largest allowed time step size [11] and denotes the largest eigenvalue of the spatial operator, which is computed up to four decimals using power iteration. The root mean square (RMS) error is computed after two time oscillations, so at s.
Figure 3 shows the RMS error plotted against the cube root of the number of degrees of freedom and the wall-clock time for the mass-lumped tetrahedral element methods using the quadrature-based algorithm for the stiffness matrix as discussed in the previous subsection. The simulations shown here were performed with an OpenMP implementation on 24 cores of two Intel® Xeon® E5-2680 v3 CPUs running at 2.50GHz. Power-law fits of the left graph are also shown in Table 9. This graph shows optimal convergence rates of order and thereby confirms the error estimates of Section 4. In particular, it confirms that optimal convergence rates are maintained, even though the spatial parameters , vary within the element.
Figure 4 shows the same as Figure 3 for the methods using exact integration to evaluate the stiffness matrix and using a piecewise constant approximation of the mass density . Power-law fits of the left graph are again given in Table 9. The graph shows that, due to the piecewise constant approximation, only second-order convergence rates are obtained. The higher-degree elements now only result in more computations per element, without any significant gain in accuracy. When comparing with Figure 3, it follows that the quadrature-based approach is much more efficient than using exact integration with piecewise-constant parameter approximations.
6.3. Elastic wave on a homogeneous domain
We also test the methods for an elastic wave propagation problem on a homogeneous domain. The elastic wave equations are given by
[TABLE]
with the displacement field, the force field, the mass density, and the elasticity tensor. We consider an isotropic elastic medium, so , with the identity tensor, the transposed of , and the Lamé parameters.
We choose domain km3 with zero Neumann boundary conditions, and set the parameters with a constant mass density g/cm3, primary wave velocity km/s, and secondary/shear wave velocity km/s. A unit vertical force source with a 7-Hz Ricker-wavelet is placed at m and receivers are placed between m and m with a 50-m interval at m and m. The exact solution can be found in [1]. The simulation time is chosen such that reflections caused by the boundary conditions do not reach the receivers.
We tested the methods on multiple unstructured meshes and simulated over the time interval s, using the time-stepping algorithm as in the previous test case and omitting the initial time steps where the magnitude of the wavelet is smaller than . Simulations were also carried out in the same environment as in the previous test case. The RMS error is based on the errors at all receivers and for all directional components and is plotted against the cube root of the number of scalar degrees of freedom and elapsed time in Figure 5. Table 10 also lists the wall clock time and number of time steps for simulations with a RMS error of around and includes results of simulations where a degree- accurate quadrature rule taken from [24] is used, with the highest polynomial degree of the enriched element space. The left graph of Figure 5 confirms that the methods converge with optimal order. It also shows that there is hardly any difference in accuracy between using the quadrature-rule approach or exact integration algorithm for evaluating the stiffness matrix. Table 10 and the right graph of Figure 5 show that for the degree-3 and degree-4 elements, the quadrature-based algorithm reduces the computational cost by more than a factor 1.5, while for the degree-2 element, this algorithm also results in a moderate speed up. Furthermore, Table 10 also illustrates that the new quadrature rules presented in this paper are more efficient than those currently available in the literature.
7. Conclusion
We presented new and efficient quadrature rules for evaluating the stiffness matrices of mass-lumped tetrahedral elements for wave propagation modelling. These quadrature rules can significantly reduce the number of computations compared to algorithms that evaluate the stiffness matrix using exact integration, and can handle spatial parameters that vary within the element without loss of the optimal convergence rate. Obtaining these quadrature rules is not trivial, since degree- mass-lumped tetrahedral element spaces contain, apart from polynomials up to degree , numerous additional higher-degree bubble functions when . To obtain efficient quadrature rules, we therefore carefully analyzed the stability and accuracy requirements needed to maintain optimal convergence rates. The resulting conditions are presented in this paper, and we prove that, if these conditions are met, the resulting method can maintain an optimal order of convergence, even when the spatial parameters vary within the element. We found quadrature rules that satisfy these conditions for recently developed mass-lumped tetrahedral elements of degrees two to four.
For the degree-2 element, the quadrature rule with the least number of points we could find was the degree-5 accurate 14-point quadrature rule of [13], but for the degree-3 and degree-4 elements, we found new quadrature rules that require significantly less integration points than those currently available. A dispersion analysis shows that by using these quadrature rules, the accuracy and largest allowed time step size remain nearly the same. Several numerical examples also illustrate the accuracy and efficiency of the quadrature-based approach and its superiority to evaluating the integrals for the stiffness matrix exactly. In particular, the quadrature-based approach results in a computational speed-up of around a factor 1.5 in case of elastic waves. Furthermore, in case of a heterogeneous domain with spatial parameters that vary within the element, the quadrature-based approach results in optimal convergence rates, while exact integration combined with a piecewise constant approximation of the spatial parameters results in a convergence rate of at most order two.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] K. Aki and P. G. Richards. Quantitative Seismology, Theory and Methods, Vol. 1. New York , 1980.
- 2[2] D. Boffi. Finite element approximation of eigenvalue problems. Acta Numerica , 19:1–120, 2010.
- 3[3] M. J. S. Chin-Joe-Kong, W. A. Mulder, and M. Van Veldhuizen. Higher-order triangular and tetrahedral finite elements with mass lumping for solving the wave equation. Journal of Engineering Mathematics , 35(4):405–426, 1999.
- 4[4] P. G. Ciarlet. The finite element method for elliptic problems . North-Holland Publishing Company, Amsterdam, New York, Oxford, 1978.
- 5[5] G. Cohen, P. Joly, and N. Tordjman. Higher order triangular finite elements with mass lumping for the wave equation. In Proceedings of the Third International Conference on Mathematical and Numerical Aspects of Wave Propagation , pages 270–279. SIAM Philadelphia, 1995.
- 6[6] G. C. Cohen. Higher-order numerical methods for transient wave equations . Springer, 2002.
- 7[7] G. C. Cohen, P. Joly, J. E. Roberts, and N. Tordjman. Higher order triangular finite elements with mass lumping for the wave equation. SIAM Journal on Numerical Analysis , 38(6):2047–2078, 2001.
- 8[8] R. Cools. An encyclopaedia of cubature formulas. Journal of Complexity , 19(3):445–453, 2003.