Higher Order Cut Finite Elements for the Wave Equation

The scalar wave equation is solved using higher order immersed finite elements. We demonstrate that higher order convergence can be obtained. Small cuts with the background mesh are stabilized by adding penalty terms to the weak formulation. This ensures that the condition numbers of the mass and stiffness matrix are independent of how the boundary cuts the mesh. The penalties consist of jumps in higher order derivatives integrated over the interior faces of the elements cut by the boundary. The dependence on the polynomial degree of the condition number of the stabilized mass matrix is estimated. We conclude that the condition number grows extremely fast when increasing the polynomial degree of the finite element space. The time step restriction of the resulting system is investigated numerically and is concluded not to be worse than for a standard (non-immersed) finite element method.

a very small intersection with the immersed domain. This will make some eigenvalues of the discrete system very small and in turn, result in poorly conditioned matrices. A suggested way to remedy this is by adding stabilizing terms to the weak formulation. A jump-stabilization was suggested in [4] for the case of piecewise linear elements, where the jump in the normal derivative is integrated over the faces of the elements intersected by the boundary. This stabilization makes it possible to prove that the condition numbers of the involved matrices are bounded independently of how the boundary cuts the elements. This form of stabilization has been used with good results in several recent papers, see for example [3,8,14,20], and has also been used for PDEs posed on surfaces in [9,5].
Thus, a lot of attention has been directed to the use of lower order elements. Higher order cut elements have received less attention so far. These are interesting in wave propagation problems. The reason for this is that the amount of work per dispersion error typically increases slower for higher order methods for this type of problems. In [14] it was suggested to stabilize higher order elements by integrating also jumps in higher derivatives over the faces. This is the intuitive higher order generalization of the stabilization first suggested in [4] and was also mentioned as a possibility in [2].
In this paper, we consider solving the scalar wave equation using higher order cut elements. Both the mass and stiffness matrix are stabilized using the higher order jump-stabilization. We present numerical results showing that the method results in higher order convergence rate. The time-step restriction of the resulting system is computed numerically and is concluded to be of the same size as for standard finite elements with aligned boundaries. Furthermore, we estimate how the condition number of the stabilized mass matrix depend on the polynomial degree of the basis functions. The estimate suggests that the condition number grows extremely fast with respect to the polynomial degree, which is supported by the numerical experiments. As a remedy for this behavior, we consider lowering the order of the elements close to the boundary. This results in a better condition number, but also in at least half an order lower convergence compared to having elements of full order everywhere. All numerical experiments are performed in two dimensions, but the generalization to three dimensions is immediate.
One reason why the considered stabilization is attractive is because it is quite easy to implement. Integrals over internal faces occur also in discontinuous Galerkin methods, thus making the implementation similar to what is already supported in many existing libraries.
The suggested jump-stabilization is one but not the only possibility for stabilizing an immersed method. In [10] a higher order discontinuous Galerkin method was suggested and proved to give optimal order of convergence. Here the problem of ill-conditioning was solved by associating elements that had small intersections with neighboring elements. Similar approaches has been used with higher order elements in for example [12,15], where elements with small intersection take their basis functions from an element inside the domain. One problem with these approaches is that it is not obvious how to choose which elements should merge with or associate to one another. A related alternative to these is the approach in [17], where individual basis functions were removed if they have a small support inside the domain. A different approach was used in [11] where streamline diffusion stabilization was added to the elements intersected by the boundary. This was proved to give up to fourth order convergence. However, this approach is restricted to in-terface problems. Another alternative is to use preconditioners to try to overcome problems with ill-conditioning, such as in [13]. However, only preconditioning does not solve the problem of severe time-step restrictions when using explicit timestepping. For this reason preconditioning alone is not sufficient in the context of wave-propagation.
This paper is organized in the following way. Notation and some basic problem setup is explained in Section 2.1, the stabilized weak formulation is described in Section 2.2, and the stability of the method is discussed in Section 2.3. Analysis of how fast the condition number increases when increasing the polynomial degree is presented in Section 2.4, and numerical experiments is presented in Section 3.

Notation and Setting
Consider the wave equation Let Ω ⊂ R d be immersed in a mesh, T , as in Figure 1. We assume that each element T ∈ T has some part which is inside Ω, that is: T ∩ Ω = ∅. Furthermore, let Ω T be the domain that corresponds to T , that is Let T Γ denote the set of elements that are intersected by ∂Ω: as in Figure 2. Let F Γ denote the faces seen in Figure 3. That is, the faces of the elements in T Γ , excluding the faces that make up ∂Ω T . To be precise, F Γ is defined as We assume that our background mesh is sufficiently fine, so that the immersed geometry is well resolved by the mesh. Furthermore, we shall restrict ourselves to meshes as the one in Figure 1, where we have a mesh consisting of hypercubes and our coordinate axes are aligned with the mesh faces. That is, the face normals have a nonzero component only in one of the coordinate directions. Denote the distance between two grid points in any coordinate direction by h.
Consider the situation in Figure 4, where two neighboring elements, T 1 and T 2 , are sharing a common face F . Denote by ∂ k n v the kth directional derivative in the direction of the face normal. That is, fix j ∈ {1, . . . , d} and let the normal of the face, n, be such that then define In the following we denote by (·, ·) X and ·, · Y the L 2 scalar products taken over the d and d − 1 dimensional domains X ⊂ R d and Y ⊂ R d−1 . Let · Z denote the corresponding norm, and let | · | H s (Z) denote the H s -semi-norm. By [v] we shall denote a jump over a face, F , that is: We shall assume that our basis functions are tensor products of one-dimensional polynomials of order p. In particular, we shall use Lagrange elements with Gauss-Lobatto nodes, in the following referred to as Q p -elements, p ∈ {1, 2, . . .}. Let V p h denote a continuous finite element space, consisting of Q p -elements on the mesh T : Define also the following semi-norm which is a norm on V p h in the case that Γ D = ∅.

The Stabilized Weak Formulation
Multiplying (1) by a test-function, integrating by parts, and applying boundary conditions by Nitsche's method [16] leads to a weak formulation of the following form: find u h such that for each fix t ∈ (0, t f ], u h ∈ V p h and where What makes this different from standard finite elements is that the integration on each element needs to be adapted to the part of the element that is inside the domain. As illustrated in Figure 5, some elements will have a very small intersection with the domain. Consider the mass-matrix in (5): Note that its smallest eigenvalue is smaller than each diagonal entry: Depending on the size of the cut with the background mesh some diagonal entries can become arbitrarily close to zero. Thus, both the mass and stiffness matrix can now be arbitrarily ill-conditioned depending on how the cut occurs. Because of this, one can not guarantee that the method is stable. One way to try to remedy this is by adding stabilizing terms, j, to the two bilinear forms where γ M , γ A > 0 are penalty parameters. This gives us the following weak formulation: find u h such that for each fix t ∈ (0, t f ], u h ∈ V p h and In [14] a stabilization term of the following form was suggested which in some sense is the intuitive extension of the stabilization which was suggested in [4]. The stabilization in (10) was also briefly mentioned as a possibility in [2]. As was discussed in [14] the bilinear form (8) can be shown to define a scalar product which is norm equivalent to the L 2 -norm on the whole background mesh: and a corresponding equivalence also holds for the gradient: The constants in (11) and (12) depend on the polynomial degree of our basis functions, but not on how the boundary cuts through the mesh. Let M denote the mass matrix with respect to the bilinear form M , and M T with respect to the scalar product on the background mesh, that is: (11) implies that the condition number, κ(M), of M is bounded by the condition number of M T : The property (12) is necessary in order to show that A(·, ·) is coercive in V p h with respect to the | · | -semi-norm on the background mesh: As we shall see in Section 2.3 this is needed in order to show that the method is stable with respect to time. The result in (14) follows by the same procedure as in [14], assuming that the following inverse inequality holds 1 This inequality follows the same scaling with respect to h and p as the corresponding standard inverse inequality, which relates the norm over a face to the norm over the whole element. See for example [22]. The stabilization in (10) is the basic form of stabilization that we shall consider. However, each time we differentiate we will introduce some dependence on the polynomial degree. It therefore seems reasonable that each term in the sum should be scaled in some way. Because of this, we consider a stabilization of the following form: where w j ∈ R + are some weights, which we are free to choose as we wish. The choice of weights will determine how large our constants C U , C L in (13) are, and in turn influence how well conditioned the mass matrix is.

Stability
The bilinear forms in (9) are symmetric. This is a quite important property, since this in the end will guarantee stability of the system. In order to show stability we want a bound over time on u Ω T . Define an energy, E, of the form Since both bilinear forms are at least positive semi-definite, this energy has the property E ≥ 0. The symmetry now allows us to show that for a homogeneous system, the energy is conserved: so that By the definition of the energy together with (11) and (14) this immediately implies that u h Ω T and ∇u h Ω T are both bounded. For the case Γ D = ∅ the semi-norm | · | is a norm for the space V p h and (14) implies that u h Ω T is also bounded. When Γ D = ∅ we can use that By integrating we obtain that u h Ω T is bounded since u h Ω T is bounded: Thus the system is stable. In total the system (9) discretizes to a system of the form with M, A ∈ R N ×N , and ξ ∈ R N , and where When solving this system in time we will have a restriction on the time-step, τ , of the form where α is a constant which depends on the time-stepping algorithm. If we for example use a classical 4th-order explicit Runge-Kutta α = 2 √ 2. The constant C F L is given by where λ max is the largest eigenvalue of the generalized eigenvalue problem: find One would expect that the added stabilization has some effect on the C F Lconstant. Because of this, we will investigate this constant experimentally in Section 3. It turns out that the C F L -constant is not worse than for the standard case (with boundaries aligned with the mesh).

Analysis of the Condition Number of the Mass Matrix
We would like to choose the weights in (16) in order to minimize the condition number of the mass matrix. This is particularly important when it comes to wavepropagation problems. For this application one typically uses an explicit timestepping method. When this is the case we need to solve a system involving the mass matrix in each time-step. In order to choose the weights we need to know how the condition number depends on the weights and the polynomial degree. To determine this, we follow essentially the same path as in [14] and keep track of the weights and the polynomial dependence of the involved inequalities. In the following, we denote by C various constants which do not depend on h or p, unless explicitly stated otherwise. We shall also by w denote the vector w = (w 1 , . . . , w p ), where w j are the weights in the stabilization term (16). We can now derive the following inequality, which is a weighted version of Lemma 5.1 in [14].
Lemma 1 Given two neighboring elements, T 1 and T 2 , sharing a face F (as in Figure 4), and v ∈ V p h , we have that: where Proof Denote by v i the restriction of v to T i and then extended by expression to the whole of T 1 ∪ T 2 . As in Figure 4, let x ∈ T 1 and denote by x F (x) the projection of x onto the face. Let n be the normal pointing towards T 1 and let j denote the only nonzero component, as in (2). We may now Taylor expand from the face: Using that Now introduce the following weighted l 2 (R p+1 )-norm: where α k > 0 and z ∈ R p+1 . If · 1 denotes the usual l 1 (R p+1 )-norm we have that: where Taking the L 2 (T 1 )-norm of (25) and using (26) now results in: Since v 2 lies in a finite dimensional polynomial space on T 1 ∪ T 2 the norms on T 1 and T 2 are equivalent: Using this in (28) and choosing gives us (23). Lemma 1 will now allow us to give a lower bound on the bilinear form M , which was defined in (8).
where L(w) is given by (24), N J is some sufficiently large integer and C l is a constant independent of h and p. Proof be a sequence of elements that need to be crossed in order to get to an element T N ∈ T \ T Γ as in Figure 6, and where we have used that L(w) ≥ 1 (since at least C 1 ≥ 1). Let now N J ≥ 1 denote some upper bound on the maximum number of jumps that needs to be made in the mesh. If our geometry is well resolved by our background mesh N J is a small integer. This gives us from which (29) follows. Fig. 6 A sequence of jumps from a boundary element T 0 ∈ T Γ to an inside element T N We proceed by estimating how a bound on the jumps depends on the polynomial degree.
Lemma 3 For the jumps in the normal derivative we have that: where T + F and T − F denotes the two elements sharing the face F .
Proof Note first that We shall need the following inequalities: which were discussed 2 in [19]. Although (33) holds for a whole element we shall use the corresponding inequality applied to a face: This is valid since a function v in the tensor product space over T will have a restriction v| F in the tensor product space over the face F . Note that the constants, C, in (33) and (34) are not necessarily the same. By combining (31), (32) and (34) we obtain (30).
Using Lemma 3 we can now bound the bilinear form M (·, ·) from above.

Lemma 4 An upper bound for
where and C g is a constant independent of h and p.
Proof Using the definition of j(·, ·) and applying Lemma 3 on each order of derivatives in the sum individually we have Let n F denote the number of faces that an element has in R d . We now have so we finally obtain: By using Lemma 2 and 4 we now have the following bound on the condition number.

Lemma 5 An upper bound for the condition number of the mass matrix is
where Proof Let λ(·) denote eigenvalues. From Lemma 2 and 4 we obtain which gives us (38).
Here, we would like to choose the weights in order to minimize the constant C M . However, we have the following unsatisfactory result, which shows that no matter how we choose the weights our analysis cannot yield a p-independent bound on the conditioning.

Lemma 6
The constant C M (w) in Lemma 5 fulfills C M (w) ≥ C 0 P (p), where C 0 does not depend on p or w. Here P (p) is the function which is independent of the choice of weights w.
Proof First note that

Now we have
where we first used that the l 1 (R p )-norm is greater than the l 2 (R p )-norm and finally Cauchy-Schwartz. From this the result follows.
The function P (p) increases incredibly fast when increasing the polynomial degree. This result could reflect either: 1. The analysis leading to Lemma 5 is not sharp. The bound C M is too generous, and a better bound exists. 2. The bound in Lemma 5 is not unnecessarily generous, so that the constant C M is in some sense "tight". This means that the condition number of the stabilized mass matrix (8) will grow faster than the function P (p), regardless of the choice of weights.
Alternative 2 is rather devastating from a time-stepping perspective, since in order to time-step (20) an inverse of the mass matrix needs to be available in each timestep. If this inversion is done with an iterative method the number of required iterations until convergence is going to be large. A combination of these two alternatives is, of course, possible. The estimate in Lemma 5 could be too pessimistic, but even the optimal bound increases incredibly fast. Given the results in Section 3 this appears to be the most plausible alternative.

Lowering the Order at the Boundary
To remedy the expected poor behavior of the condition number of the mass matrix we shall consider lowering the order of the elements close to the boundary. This will be done in the way illustrated in Figure 7. This idea is based on two observations: -In finite difference methods it is possible to lower the order close to the boundary and still get full convergence [6,21]. -By using lower order elements close to the boundary we only need to stabilize jumps in derivatives up to order p − 1.
Let N F (T ) denote the neighboring element of the element T sharing the face F with T . We now construct a new finite element spaceṼ p h in the following way. Elements which are intersected or have an intersected neighbor are lowered one order compared to the interior of the domain. More precisely: Otherwise .
(40) Using this space we can still guarantee stability.
The space in (40) will introduce hanging nodes between elements of different orders. This can be solved in several ways, but in the experiments in Section 3 we treat this by adding constraints that enforce continuity at the hanging nodes. Fig. 7 Illustration of the spaceṼ p h , that was constructed by lowering the order of the elements close to the boundary.

Choosing Weights in the Jump-Stabilization
In order to do a computation, we are forced to make some choice of the weights w i . The essence of Lemma 6 is that we can bound L(w)G(w) from below. So in order to choose weights let us assume that: From Lemma 4 it is seen that choosing w i 1 makes G(w) very large. In the same way, Lemma 2 tells us that choosing w i 1 for some i makes L(w) very large. From this observation it seems reasonable to try to enforce both bounds to be of about the same magnitude. In this way, we minimize L(w)G(w) with respect to w and enforce G(w) = L(w). This leaves us with where ∇ w denotes the gradient with respect to w. This now gives us the following choice of weights: There is no reason why this argumentation should lead to the optimal choice of weights, but it seems reasonable that this is not a particularly bad choice.

Numerical Experiments
In the following, we shall solve both an inner problem and an outer problem using finite element spaces of different orders. For the inner problem we use both the standard spaces V p h , defined by (4), andṼ p h , the corresponding spaces with lower polynomial order in the elements close to the boundary. For the outer problem we only use the space V p h . The weights from (41) are used, with p determined by the order of the polynomials at the boundary. In addition, the following parameters are used The scaling of γ D with respect to p follows from the inequality (15). When p = 1 these parameters coincide with the parameters used in [20]. There the effect of γ M on the condition number of the mass matrix was investigated numerically. For p = 1 this choice of γ A and γ D also coincides with the one in [4], where γ A was investigated numerically. The errors are computed in norms which are some quantities integrated over the domain Ω. It is worth noting that the geometry of Ω is represented by a level set function, ψ h . Both for the case when u ∈ V p h and when u ∈Ṽ p h the level set function is an element in the space where and where T B is a larger background mesh from which T was created. In order to perform the quadratures over the elements intersected by the boundary we have used the algorithm in [18], which generates the quadrature rules on the intersected elements with respect to ψ h . Thus, also the errors of the solution are calculated with respect to this approximation of the geometry. That is, the L 2 -norms are approximated as where ψ h is initialized by L 2 -projecting the analytic level set function onto the space W p h . Convergence-rates are estimated as where e i denotes an error corresponding to mesh size h i . Time-stepping is performed with a classical fourth order explicit Runge-Kutta, after rewriting the system (20) as a first order system in time. A time step, τ , of size is used. During the time-stepping we need to solve a system involving the mass matrix. When using higher order elements the condition number of the mass matrix is large, so an iterative method requires a lot of iterations. Because of this, a direct solver was used. On non-cut elements Gauss-Lobatto quadrature was used to assemble the mass matrix. This makes the mass matrix almost diagonal, which makes use of a direct solver cheaper. All off-diagonal entries in the mass matrix are related to degrees of freedom close to the immersed boundary. The library deal.II [1] was used to implement the method.

Standard Reference Problem with Aligned Boundary
It is relevant to compare some of the properties of the mass and stiffness matrix with standard (non-immersed) finite elements. The unstabilized mass and stiffness matrix were computed on a rectangular grid with size [−1.5, 1.5] × [−1.5, 1.5], with Neumann boundary conditions. As for the immersed case, quadrilateral Lagrange elements with Gauss-Lobatto nodes were used. The computed C F L is shown in Table 1. The C F L -number was computed according to (22). For a given p the value in Table 1 is the mean value when calculating the CFL-number over a number of grid sizes. The condition number of the mass matrix is shown in Figure 8 and the minimal and maximal eigenvalues of the mass matrix is shown in Figure 9. Since all eigenvalues should be proportional to h 2 , the eigenvalues have been scaled by h −2 for easier comparison.

An Immersed Inner Problem
Let Ω be a disk domain, centered at origo, with radius R = 1, and enforce homogeneous Dirichlet boundary condition along the boundary u| ∂Ω = 0.
Let J 0 denote the 0:th order Bessel-function and let α n denote its n:th zero. By starting from initial conditions: we can calculate the error in our numerical solution with respect to the analytic solution: u(x, t) = J 0 (α n x R ) cos(ω n t), ω n = α n R .
A few snapshots of the numerical solution are shown in Figure 10. The problem was solved with the given method until an end-time, t f , corresponding to a three periods: At the end-time the errors were computed. The calculated errors and estimated convergence rates for the different element orders are shown in Tables 2 to 4. One would expect that the L 2 (Ω)-and L 2 (∂Ω)errors are proportional to h p+1 and that the H 1 (Ω)-error is proportional to h p . The rates agree quite well with this. Computed C F L -constants for different element orders are shown in Table 1. The values were computed according to (22). We see that the C F L -constant is essentially the same as for the non-immersed case. In the same way as for the non-immersed case, Table 1 displays the mean value over a number of grid sizes, but the C F L -number only varied slightly when varying the grid size. By inserting the values in Table 1 into (21) one can see that it would have been possible to use a larger time-step than the one in (42).
How the condition number of the mass matrix depend on the grid size is shown in Figure 8, for the different orders of p. We see that the condition numbers are essentially constant when refining h, in agreement with (11). We also see that the condition numbers increase extremely rapidly when increasing the polynomial degree, as predicted by Lemma 6. It is also clear from Figure 8 that the condition number increase much faster than in the non-immersed case. The dashed lines in Figure 8 denote the function CP (p), where P is the function from (39). The constant C was chosen so that CP (1) agreed with the mean (with respect to h) of the condition numbers for V 1 h . The estimate from Lemma 6 appears to be too pessimistic.
The minimal and maximal eigenvalues for the different polynomial orders and refinements are seen in Figure 9. As can be seen, the scaled eigenvalues are essentially constant with respect to h. Thus the dependence on h is in agreement with the theoretical considerations in Section 2.4. We see that the minimal eigenvalues decrease quite fast when increasing the polynomial degree, and that they are substantially smaller than in the non-immersed case. The maximal eigenvalues also decrease but much slower than in the non-immersed case.   The errors and convergence rates when using the spaceṼ 2 h are shown in Table  5. For the errors in L 2 (Ω)-and H 1 (Ω)-norm it seems that we obtain the rate corresponding to the highest element (Q 2 ) in the space. Not unexpectedly we seem to get the lower order convergence for the L 2 (Ω)-error along the boundary. However, when using the spaceṼ 3 h the situation appears to be different. Here, it seems that one looses at least half an order for the rates of the L 2 (Ω)-and H 1 (Ω)-errors, which is rather unsatisfactory.
How the condition number of the mass-matrix depends on h for the two spaces V 2 h andṼ 3 h are shown in Figure 11. By comparing to Figure 8 we see that the condition number of the space V p h is essentially the same as for the corresponding space with the lowest order element everywhere. That is: which is not surprising since we expect that all ill-conditioning is due to the added penalty term, j(·, ·). The minimal and maximal eigenvalues of the mass matrix for the spaceṼ p h look essentially the same as for the space V p−1 h , while the CFLnumbers for the spaceṼ p h look essentially the same as for the space V p h .  Fig. 11 Condition number of the mass matrix when using the spaceṼ p h , for different h and p. The dashed lines denotes estimates according to the function P (p).
Here, g D is the function where we have chosen σ = 0.25, t c = 3. A few snapshots of the numerical solution is seen in Figure 13.
Here, we don't have an expression for the analytic solution. So when computing the errors we compare against a reference solution. The reference solution was computed on a grid twice as fine as the finest grid that we present errors for.
Given the less satisfying results using the spaceṼ p h in Section 3.2, we here only examine the convergence results for the space V p h . The computed errors after solving to the end time t f = 4 are shown in Table 7 to 9.
We see that the convergence is at least h p+1 for the L 2 (Ω)-error and h p for the H 1 (Ω)-error. In Table 7 to 9 the last column is the L 2 (Γ N )-error in the Neumann boundary condition which we see is close to the expected rate h p .

Discussion
The results in Section 3.2 and 3.3 show that it is possible to solve the wave equation and obtain up to 4th order convergence. In particular, it is also promising that the CFL-condition is not stricter than for the non-immersed case. However, both the theoretical results in Lemma 6 and the results in Section 3.2 show that there are problems with the conditioning of the mass matrix. It should be emphasized that even if the added stabilization creates some new problems it is by far better than using no stabilization at all. With the added stabilization the method can be proved to be stable, which is essential. It would, of course, be advantageous if one would be able to create a stabilization which does not lead to conditioning problems. However, the prospects for creating a good preconditioner for the mass matrix is rather good, since the stabilization maintains the symmetry of the mass matrix and since one obtains bounds on its spectrum from the analysis.
The choice of the weights in (41) were based on hand-waving arguments and can, therefore, be criticized. We have tried other choices of weights but have not presented the results here. This is mainly because they give similar results and we have no reason to believe that there exists a choice which makes the condition number significantly better.
The idea of lowering the order of the elements close to the boundary worked quite well for the spaceṼ 2 h . We obtained the convergence corresponding to the higher elements in the space, but the condition number corresponding to the lower order elements. Unfortunately, higher order convergence was not achieved when increasing the element order beyond 2. Thus, the procedure does not appear to be a plausible solution for going to higher orders.