High Order Cut Finite Elements for the Elastic Wave Equation

A high order cut finite element method is formulated for solving the elastic wave equation. Both a single domain problem and an interface problem are treated. The boundary or interface are allowed to cut through the background mesh. To avoid problems with small cuts, stabilizing terms are added to the bilinear forms corresponding to the mass and stiffness matrix. The stabilizing terms penalize jumps in normal derivatives over the faces of the elements cut by the boundary/interface. This ensures a stable discretization independently of how the boundary/interface cuts the mesh. Nitsche's method is used to enforce boundary and interface conditions, resulting in symmetric bilinear forms. As a result of the symmetry, an energy estimate can be made and optimal order a priori error estimates are derived for the single domain problem. Finally, numerical experiments in two dimensions are presented that verify the order of accuracy and stability with respect to small cuts.


Introduction
The elastic wave equation is important in several applications. For example, materials in the earth's crust can be modeled as linear and elastic, and earthquakes give rise to seismic waves that propagate through the crust. Other examples include non-destructive testing and propagation of waves in beams and other solid structures. High order accurate methods are especially attractive when solving the elastic wave equation. The reason is that high order methods, in general, have lower work per dispersion error [11]. Seismic waves typically propagate over large distances and are therefore prone to dispersion error. Also, elastic waves often propagate in media with complicated geometries. So, it is of interest to have numerical methods with high order of accuracy that can handle complicated geometries. Examples of such methods are discontinuous Galerkin (dG) methods [6,16], and summation by parts (SBP) based finite difference methods [7,1]. The dG methods usually handle the complicated geometries by using an unstructured grid that conforms to the boundary, meanwhile the SBP-based finite difference methods use a curvilinear grid to handle the complicated geometries. While both of these methods work very well, it may at times be hard to find a good curvilinear mapping and it can be cumbersome to generate a good conforming grid.
In the present paper, we are interested in solving the elastic wave equation using the cut finite element method (Cut-FEM) with high order elements. Cut-FEM is an immersed method where boundaries and interfaces do not need to be aligned with the computational mesh. For details on Cut-FEM see for example the review paper [4]. Cut-FEM with high order elements has been studied in for example [8,9,18]. When using high order elements in Cut-FEM a few difficulties emerge. One problem is generating a high order quadrature on the elements that are cut by the boundary or interface. To accomplish this we use an algorithm by Saye [17]. Further, in the same way as for standard non-cut finite elements, the time step restriction becomes more severe when the element order is increased, which makes the time stepping more expensive. Finally, stabilization terms are added (introduced in [3,5,14]) in order to make the eigenvalues of the matrices bounded from above and below, independently of how the boundary/interface cuts the mesh. Unfortunately, this stabilization makes the condition number of the mass matrix increase fast when the element order increases. This can make time stepping the discrete system more expensive since we need to solve a system involving the mass matrix during time stepping. The present paper builds on the work in [8], where time-independent elasticity equations were solved using the Cut-FEM technique.
There are several reasons for using Cut-FEM to solve the elastic wave equation. One example is when a boundary or interface has a complicated geometry. Creating a computational mesh that conforms to this geometry can be expensive and timeconsuming. Using an immersed method could potentially be cheaper. Another example is when the geometry of a boundary or interface is not known a priori. This could, for example, be the case if the geometry of the interface is hard or impossible to measure. One way to get around this is to send waves from the surface that propagate toward the interface and get reflected from it. By measuring the reflected waves it is possible to solve an inverse problem to compute the shape of the interface. In order to solve the inverse problem, one would need to iterate over a lot of different interface geometries. Here, an immersed method would be useful since remeshing the interface geometry could be very time-consuming. Similar inversion problems have been of interest for some time and were introduced partly by Tarantola in the papers [20,21,22].
The present paper is organized as follows. In Section 2 the mathematical problems are stated. These are the elastic wave equation posed on a single domain and as an interface problem. This is followed by the explanation of the method in Section 3. In Section 4 we present a proof of convergence for the single domain problem, and in Section 5 we present numerical results on the order of conver-gence and robustness with respect to small cuts. Finally, we end with a discussion in Section 6.

Model of the Problem
We are interested in the elastic wave equation posed both on a single domain Ω ⊂ R d (Figure 1a), and as an interface problem on a composite domain Ω = Ω 1 ∪ Ω 2 ⊂ R d (Figure 1b). The interface problem is interesting when we have two materials in contact with each other, which occurs frequently in applications due to the layered structure of the earth's crust. On the other hand, the single domain problem is relevant if we have an inclusion of air or vacuum inside another material.  Let n denote the outward unit normal to ∂Ω, and assume that ∂Ω is partitioned such that ∂Ω = Γ N ∪ Γ D , with Γ N ∩ Γ D = ∅. The single domain problem reads: where u is the displacement vector, ρ is the density and σ is the stress tensor. We shall assume that Γ D and Γ N are sufficiently smooth. Furthermore, we assume that we are working with a linear, homogeneous and isotropic material. When this is the case the stress in the material is given by where δ ij is the Kronecker delta function and is the strain tensor defined as In (6) λ and µ are the Lamé-parameters, which are material dependent scalar constants.

Interface Problem
Consider now an interface problem on the domain illustrated in Figure 1b. We have a composite domain consisting of two elastic materials with material-parameters ρ i , λ i , µ i . In this case, the problem is given by where u i is the displacement vector in material i and the stress and strain tensors are defined analogously to (6) and (7). We assume that Γ D , Γ N and Γ I are sufficiently smooth. Here n i is the outward normal to Ω i and n is the normal pointing from Ω 2 to Ω 1 (n = n 2 ). · defines the jump over the interface: Since we have several normals defined (n 1 , n 2 and n), (12) can be interpreted in two different ways. To avoid any confusion we use the convention that the normal is fixed σ(u) · n = σ(u 2 ) · n − σ(u 1 ) · n, x ∈ Γ I .

Numerical Method
Let Ω be covered by a background mesh, T B , as in Figure 2a. We shall only consider the case when the mesh consists of quadrilaterals that are squares and of the same size. Let h denote their side length. Let the boundary or interface be partitioned as illustrated in Figure 2a. That is, for the single domain we assume where Γ A is aligned with the boundary of the mesh while Γ C cuts through it. Correspondingly for the interface problem, we assume that ∂Ω ∪ Γ I = Γ A ∪ Γ C . Let T C denote the elements that are intersected by Γ C : as illustrated in Figure 2b.  Let i ∈ {1, 2} denote an index indicating the domain, which will be omitted for the single domain problem. Let T i , denote the smallest set of elements in the background mesh covering Ω i , as illustrated in Figure 3. In particular, for the single domain problem, T is the smallest set of elements covering Ω. To be precise let Now introduce the spaces Where Q p (T ) denotes the p:th order Lagrange element with Gauss-Lobatto nodes over T . For high element orders, Gauss-Lobatto nodes result in a mass matrix with better properties than if equidistant nodes are used [10]. For the single domain problem, we solve for the solution u h ∈ V h , while for the interface problem we solve for the pair For the interface problem, this means that the degrees of freedom are doubled over the elements in the set T C .
Since the weak formulations for the single domain and the interface problem are very similar we discuss their derivation more or less simultaneously. We shall use the following standard inner products where the subscripts indicate over which domain the integration takes part. If u or v in (20) are tensors then contraction to a scalar is implied. Note that the angular brackets denote integration over a curve in 2D (or surface in 3D). By multiplying (1) or (8) by a test function, integrating by parts and simplifying (for details see for example [13]) we get Note that the Dirichlet boundary conditions are consistent with the following terms So in order to enforce the boundary conditions by Nitsche's method we add (22)-(24) to (21). Here, γ D is a constant controlling how strongly the Dirichlet boundary condition is enforced. We now have and In (25) the term B i corresponds to integration over the "bulk" and the terms D i and L D i enforce the Dirichlet boundary condition over Γ D i : Note that the terms (22)-(24) were added in a way so that a i in (26) is a symmetric bilinear form. Now (25) is the starting point for the weak formulations for both the single domain and the interface problem. Note also that for the single domain we have while for the interface problem

Stabilizing Small Cuts
A common problem for immersed methods is robustness with respect to small cuts. In order to understand this problem consider the single domain. Since Γ C intersects the mesh in an arbitrary way an element K may have an arbitrarily small intersection with the domain so that the size of K ∩ Ω h d . For each element we integrate over K ∩ Ω. For the mass matrix this means that the smallest eigenvalue can be arbitrarily small, and in turn that the condition number can be arbitrarily large. For the stiffness matrix, the problem is even worse. The term (22) that we add to enforce the boundary condition can make some eigenvalues of the stiffness matrix negative, which would make the method unstable.
A suggested way to remedy this problem is to add a stabilizing term, j i , both to the term that corresponds to the mass matrix and to the term that corresponds to the stiffness matrix: Here, γ i M and γ i A are scalar constants that control how much stabilization is added. In order to explain the definition of j i let F i denote the faces illustrated in Figure 4. That is, the faces of T C excluding the boundary faces of T i . To be precise let We now define the stabilization term as Here, ∂ k n v i denotes the k:th derivative in the direction of the face normal, n, and [·] defines the jump over a face F : Note that [·] is different from · in (15) since we have u i on both sides of F . The stabilization in (34) was suggested first in [3] and used first for the Poisson equation in [5]. For a nice explanation of why it works see [14]. With stabilization one can prove the following inequalities for the bilinear form where Ω i is defined as the domain that T i covers: In (36) C L and C U are positive constants that depend on the element order but not on h. From (36) we immediately get that the eigenvalues of the stabilized mass matrix are bounded independently of how the boundary/interface cuts the mesh. In turn, this bounds the condition number independently of the location of the boundary/interface. Unfortunately (as noted in both [18] and [8]) the constant in the bound increases very fast with the order of the elements. With stabilization, one can also show that the bilinear form A is continuous and coercive independently of how the boundary/interface cuts the mesh. This result and the one in (36) were proved for the time-independent elasticity equations in [8].

Weak Form for the Single Domain Problem
For the single domain we have that ∂Ω = Γ D ∪ Γ N , so by starting from (25) and adding the stabilizing terms we get the weak form for the single domain problem:

Weak Form for the Interface Problem
We now want to derive the weak formulation for the interface problem (8)- (14). First, let κ 1 > 0 and κ 2 > 0 fulfill κ 1 + κ 2 = 1 and let {·} to denote the following convex combination: By using that Ω i \ (Γ D i ∪ Γ N i ) = Γ I , n = n 2 = −n 1 and the condition (12) it is straightforward to verify that Note also that the interface condition (11) is consistent with the following terms Here, γ I is a positive constant which will control how strongly the interface condition is enforced. Now we add (25) for each domain, use (40) and add (41), (42) and stabilization to obtain the finite element method: (31), (27) and (28). The bilinear form I that enforce the interface conditions is given by The method contains a number of free parameters that need to be chosen. Clearly, the penalty parameters related to the stabilization should scale with the parameters of the materials. We choose to scale them as where We choose the constants related to the interface terms in the following way The scaling with respect to η i is analogous to the choice of parameters for the Poisson interface problem in [4]. The Nitsche parameter related to the Dirichlet boundary condition is chosen as Here, the scaling with p 2 of γ D and γ I follows from an inverse inequality. The numerical constants are chosen based on experience. We shall briefly discuss this in Section 6.

Imposition of Initial Conditions
In order to impose the initial conditions we first define the stabilised L 2 -projection, Π h u. For the single domain problem, Π h u is defined as the solution to the following problem: Given u, find Π h u ∈ V h such that For the interface problem, Π h u is defined analogously as the solution to: Given u, The initial conditions are now imposed as Note that, by setting the discrete initial conditions in this way, the initial conditions of the single domain problem, (4)- (5), only need to be defined on Ω and not on Ω .

Theory
In this section, we will present some theoretical results, in particular, a proof of convergence for the semi-discrete method for the single domain problem. The proof builds on the results presented in [8] where several time-independent problems were studied. During the analysis we will use the following norms: where we can note that | · | j is a semi-norm. Note that these norms only make sense if the argument is defined on Ω * . We will also use the -relation, which we define as where C is some constant that is independent of h. We also need a bounded extension operator, E : H s (Ω) → H s (Ω ). We shall assume that the solution is sufficiently smooth (s is sufficiently high) and that ∂Ω is sufficiently regular so that

Ritz Projection
In order to prove convergence we need a "Ritz-like" projection, which we define as the solution to the following problem: Given u, find R h u ∈ V h such that In this section, we will gather some results about the Ritz projection, which will be essential in the analysis to come. For brevity, we will from here on omit the "like" in the Ritz-like projection (56) and simply call it the Ritz projection. As shown in [8], given that γ D is sufficiently large, A is coercive and continuous with respect to |||·||| h . That is, there exists constants C r , C c > 0 such that For simplicity, we will assume that Γ D = ∅. When this holds, |||·||| h is indeed a norm (i.e. not only a semi-norm) and (56) has a unique solution. However, this assumption can likely be relaxed by looking for the solution in a constrained subspace of V h . One should note that this projection is nothing but the solution to the timeindependent elasticity problem. To see this, letû(x) = u(x, t f ), where t f is some fixed time, and definef so thatû is the solution to This means thatû will satisfy whereL is defined aŝ i.e. the same as L in (27) but using the right hand side data from (58). We can now formulate the finite element method to solve (58) as: Findû h ∈ V h such that Now, by subtracting (59) from (60) we can see that the solutionû, to the problem (58), in fact corresponds to the Ritz projection R h u in (56). So in principle the Ritz projection is obtained by solving a linear elasticity problem. This has been treated in detail in [8], where the results presented in Lemma 1 were derived.
Lemma 1 For the Ritz projection, R h u, in (56) the following error estimates hold Proof See Theorem 4.2 in [8].
We shall also need the following corollary.
Corollary 1 For the Ritz projection, R h u, in (56) the following holds Proof From (61) and the definition of |||·||| h in (53) we get from which (63) follows.

A priori Analysis
The analysis presented here is similar to the one presented in [19]. We wish to bound the error u h − u and in doing so we split the error in two parts, where e N = u h − R h u and e R = R h u − Eu. By Lemma 1 we directly get a bound for e R . To bound e N we first aim to find a bound on the "energy" of e N , which we define as To facilitate the proof, we will in this section assume that the discrete initial conditions are imposed using the Ritz-projection: Note that (66) is not the same initial conditions as in (51), which are used in the numerical experiments. The reason for this is that computing the Ritz-projection is more involved than computing the L 2 -projection. In practice, this most likely makes no difference since the result of both projections approximates the analytical solution with the same order of accuracy. However, the choice (66) makes the analysis simpler since it is equivalent to which by the definition of the energy in (65) gives us We are now ready to bound the energy.

Lemma 2 The following bound holds
Proof First, we have that When going to the third line we used the definition of the Ritz projection in (56). Finally we used (55) and the definition of e R . Now, choosing v h =ė N in (71) we can use the definition of the energy and that M is an inner product (so that Cauchy-Schwarz applies) to get By using we can divide both sides of (72) by 2 E e N and get where we in the last line used Lemma 1 and Corollary 1. Integrating and squaring (73) gives Finally, using (69) gives us the bound in (70).
We are now ready to state our a priori error estimates. They are summed up in Theorem 1.
Theorem 1 Let u be the solution to (1)-(5) and let u h be the solution to (38), then at any given time, t, the following a priori error estimates hold Proof Using the definition of E e N and Lemma 2 we get In order to bound e N and notė N note that Dividing (79) by 2 e N Ω and integrating over time gives by using (67). Combining (80) with (77) gives us Finally, we use the triangle inequality on (64) and combine (78) and (81) with the bounds on e R from Lemma 1 to get the estimates (75) and (76).

Time Step Restriction
Both of the weak forms (38) and (43) will discretize to a system of the form where M ∈ R N ×N is the mass-matrix, A ∈ R N ×N is the stiffness-matrix and L ∈ R N is the right-hand side vector. If we use explicit time stepping the largest time step, τ , we can take due to stability restrictions will be bounded by the C F L -number as: where α is a constant which depends on the chosen time stepping scheme. The C F L -number can be computed from the matrices in the discrete system. Let λ max be the largest eigenvalue of the generalized eigenvalue problem: find λ, x ∈ R N such that Ax − λMx = 0.
Then the C F L -number is given by It is important that the C F L -number does not decrease significantly when the smallest cut in the mesh approaches zero. Ideally, the time step restriction should not be more severe than for the standard non-cut finite element method.

Material Parameters
The problem for the single domain contains three material parameters, ρ, λ and µ. However, by rescaling (see [12]) one can show that the dimensionless equation only depends on the ratio, β, between the Lamé-parameters: Thus we can without loss of generality assume that the equation is already in dimensionless form and set ρ = µ = 1. Now we can obtain different physical behavior by varying λ. For the interface problem, we shall also assume that we are working in dimensionless form. By a corresponding analysis, it is possible to show that we can set ρ 1 = µ 2 = 1 and obtain different physical behavior by varying λ 1 , ρ 2 , λ 2 and µ 2 .

Numerical Experiments
In this section, we present some numerical examples. First, we investigate if the error converges with the expected order. This is done for the single domain problem in Section 5.1 and for the interface problem in Section 5.2. In Section 5.3 we investigate how the properties of the discretized matrices in (82) change when the smallest cut in the mesh approaches zero. To implement the method, we have used the finite element library deal.II [2]. A level set function has been used to represent the immersed boundary/interface. To generate high order quadrature rules on the intersected elements we have used the algorithm from [17].
In the experiments below the following material parameters have been used These parameters correspond to material 1 being sandstone and material 2 being granite, these are two of the most common rock types. Note that by using the present model we have assumed that the materials are linear, homogeneous and isotropic, which are possibly unrealistic for these types of rock. For waves in elastic materials, two different wave speeds are of importance. The pressure-, c p , and shear-wave speed, c s . These relate to the material parameters as The parameters in (87)  For time discretization we have used the explicit fourth order accurate classical Runge-Kutta. In the experiment below a time step, has been used. Since the condition number of the mass matrix is expected to be large a direct solver was used to invert M during the time stepping of (82).

Convergence for the Single Domain Problem
Assume that we have an elastic pressure wave traveling through R 2 in the xdirection: Here, ω is a constant which we choose as ω = π. Let this wave hit a circular inclusion (vacuum inside) with radius, R = 1. At the boundary of the inclusion, Γ N , a homogeneous Neumann boundary condition is enforced. If we consider this problem in all of R 2 the total solution,ũ, will be the sum of the incoming, u in , and reflected wave, The reflected wave can be computed analytically. The total analytical solution (given in [23]) is periodic in time and can be written as a series expansion in Bessel and Hankel functions. In this paper, we truncate the series and use it as our solutionũ. Since the solution is rather complicated we do not restate the series-expansion here, but merely refer the interested reader to [23]. Consider now the single domain problem in (1)-(3) posed on the finite domain as in Figure 1a. We have a finite square domain with side length L = 2π. As in Figure 2a, the outer boundary is aligned with the mesh but the inner boundary is not. We want to make the solution, u, on this truncated domain equal to the analytical solution,ũ, on R 2 . To achieve this, we set the initial conditions equal toũ: and impose a Dirichlet boundary condition on the outer boundary equal toũ: We solve this problem until the end time T = 2 (corresponding to one period) and compute the L 2 -error for decreasing mesh sizes. Snapshots of the solution at the initial time and a quarter of a period later are shown in Figure 5. The error in L 2 -norm as a function of element size is shown in Figure 6 for Q 1 -to Q 3 -elements. The straight lines in the figure denote the expected order of accuracy. We see that the order is a bit low for large h, but when going to finer h we get the expected order or even slightly higher order than expected.

Convergence for the Interface Problem
Consider now a similar setup as in Section 5.1. We have a plane wave of the form (89) traveling through a material in R 2 towards a disc. The material has properties ρ 1 , λ 1 , µ 1 , and the disc has radius equal to 1. However, instead of vacuum, we replace the material of the disc by another material with properties ρ 2 , λ 2 , µ 2 . In the same way as before, the reflected wave can be solved for analytically and the total solution,ũ, can be found in [23] in the form of a series expansion. We again truncate the series and use it as our solution. Now we solve the interface problem (8)-(12) posed on the finite domain in Figure 1b. Again we have a square domain with side length 2π. To make the solution of the problem equal to the analytical solution we again set the initial condition and the outer Dirichlet boundary condition equal toũ, as in (91)-(92). Snapshots of the solution at two different times are seen in Figure 7. We see that the displacement in the x-direction looks like the plane wave in (89), but since the wave-speed is lower in Ω 2 the plane wave gets distorted.
To verify the convergence we solve until the end time T = 2 (corresponding to one period) and then compute the error. The error in L 2 -norm as a function of element size is seen in Figure 8. We see that the order of accuracy is as expected for Q 1 -and Q 2 -elements. For Q 3 -elements the order is a bit low for large h, but eventually reaches the expected order when we go to finer h.

Matrix Properties with Decreasing Cut-Size
Consider the setup illustrated in Figure 9a for the single domain and in Figure 9b for the interface problem. For both setups, we have a rectangular domain on top of a square grid. For the single domain problem in Figure 9a the left, bottom and top boundary are aligned with the mesh, but the right domain boundary intersects the last column of elements with a cut of size h cut . For the interface problem, all boundaries are aligned with the mesh boundaries, but the immersed interface intersects the middle column of elements with a cut of size h cut . We are now interested in how the properties of the mass and stiffness matrix change when we vary the size of h cut . In the experiment, we use a background mesh containing 9 × 9 elements, which is slightly finer than what is illustrated in Figure 9.  How the condition number of the mass matrix changes is seen in Figure 10a for the single domain problem. We see that when the cut size is large (h cut /h ≈ 1) the condition number is small and initially grows when h cut is decreased. However, as the cut-size is decreased further the condition number becomes constant, as expected from the theory. We also see that the constant level increases very fast when we increase the order of the elements, which is consistent with results previously presented in [18,8].
In Figure 10b we see the condition number of the mass matrix for the interface problem. Note that we have f (h cut /h) on the x-axis, where f (x) = log 10 (x) − log 10 (1 − x) .
This makes the x-axis "almost logarithmic" as h cut /h approaches both 0 and 1, since f (x) is monotone on the interval (0, 1) and maps (0, 1) to (−∞, ∞). In Figure 10b we see that the behavior is analogous to the single domain problem as h cut /h approaches 0. We also see that the curve is almost mirrored in the point h cut = h/2. That the curve is not exactly mirrored can be explained by the difference in material parameters.
In the same way, the condition number of the stiffness matrix is seen in Figure 11a and 11b. We see that the dependence is similar as for the mass matrix in Figure 10a and 10b.
The C F L -number computed from (85) is shown in Figure 12a for the single domain problem and in Figure 12b for the interface problem. We see in the figures that the C F L -number is completely independent of the size of the cut. We also see that the C F L -number becomes smaller when we increase the order of the elements. This is also the case when using the standard (non-cut) finite element method.       The numerical experiments in Section 5.1 and 5.2 show that the method converges with the orders expected from Theorem 1. Furthermore, from the experiment in Section 5.3 we see that the method is robust when the size of the smallest cut in the mesh approaches zero.
The parameters (87) of the two materials used in the experiments for the interface problem are different but do not differ significantly. A future possibility would be to test how more extreme differences in material parameters affect the performance of the method. For the interface problem, the limit µ 2 → 0 is particularly important. For this case, material 2 stops being elastic and the problem on Ω 2 becomes equivalent to the acoustic wave equation [15]. One disadvantage of taking the limit µ 2 → 0 is that the problem on Ω 2 still is a system. Thus one future research direction would be to consider the problem of the elastic wave equation coupled directly with the acoustic wave equation.
The choice of numerical constants in front of γ i M , γ i A , γ I and γ D in (45), (47) and (48) is rather arbitrary. As far as we have seen the method is not particularly sensitive to the choice of constants. Still one can wonder what happens when they are chosen differently. If γ D and γ I are chosen too small coercivity is lost and the method becomes unstable, due to eigenvalues of the stiffness matrix becoming negative. This has nothing to do with the method being immersed. The same thing occurs also when symmetric Nitsche techniques are used in non-cut methods. Generally one wants to choose γ D and γ I close to the stability limit. If they are chosen larger than necessary the C F L -number becomes smaller. The influence of the stabilization parameters γ i A and γ i M on the condition numbers of the mass and stiffness matrix were discussed in [5,19], for linear P 1 -elements. There one could see that the condition numbers had a minimum when either stabilization parameter increased from 0. However, the condition number of either matrix increased rather slowly after passing the minimum. Thus, choosing γ i M or γ i A slightly larger than necessary does not have a severe effect.
As mentioned earlier, high order methods are typically attributed to being more efficient for hyperbolic problems. We have not investigated whether this is the case for the present method, but there are several aspects that would affect the efficiency. When increasing the order of elements the order of the quadrature must also be increased. Creating quadrature rules on the intersected elements is typically expensive, and using more quadrature points means more work. Whether it pays off to increase the order likely depends on what algorithm is being used to generate the quadrature. However, when solving wave propagation problems we are often interested in solving for an extended period of time. When this is the case the time spent on time integration is typically dominant. When time stepping (82) we need to be able to invert the mass matrix. If the number of degrees of freedoms is not too large we can afford to factorize it. Once factorized, inverting the mass matrix is very fast. However, if the number of degrees of freedom is very large we are forced to use an iterative method. This is potentially not efficient since we saw in Section 5.3 that the condition number is very large when the element order is high.