Standard finite element methods approximate the solution of a given partial differential equation (PDE) by piecewise polynomial functions. A classic approach to the discretization of time dependent PDEs is to use finite element methods to discretize space and then use time stepping schemes to advance in time. We consider here a different approach based on finite element approximation simultaneously in space and time. This requires to mesh the full space-time domain. The use of approximation spaces based on piecewise ”total degree” polynomials in both space and time, leads to a higher number of degrees of freedom. On the upside,-refinement is made possible in space-time, allowing for straightforward higher order approximation. Furthermore, the space-time domain mesh is not forced to be a product mesh, as it is for time stepping schemes. Instead, we are allowed to use unstructured meshes. This gives us the possibility to devise suitable mesh design strategies in order to circumvent the CFL-condition, which usually limits the global time-step size by the size of the smallest spatial element in explicit time stepping schemes. In this paper we study a space-time discontinuous Galerkin (DG) method that combines tent pitched meshes [28, 8] and Trefftz basis functions . This will allow to advance in time elementwise and in parallel.
Space-time finite element methods for linear wave propagation go as far back as , and have been used with DG methods e.g. in [20, 9, 17, 5]. DG methods are based on discontinuous piecewise polynomial functions, and so-called numerical fluxes which impose continuity constraints at mesh inter-element boundaries.
On one hand, ideas of combining DG methods with so called tent pitched meshes can be found in e.g. [21, 20, 9, 18, 27]. Tent pitching techniques give a possible way of generating a space-time mesh, which complies with the causality properties of the hyperbolic PDE. The resulting mesh consists of tent shaped objects, each advancing locally optimal in time, with the PDE being explicitly solveable in each of them. Though the tent pitching strategy pairs well with DG methods, also other methods are applicable in combination with tent pitched meshes. In [11, 26], a scheme for semi-discretization for hyperbolic equations on tents is presented. On the other hand, in , Friedrichs theory is used in order to derive a conforming method and prove its convergence properties. We point out that tent pitching is not the only way to deal with the time step restiction of locally refined meshes. A stabilization for a conforming space-time finite element method on Cartesian (in time) meshes is presented in . Classical time-stepping schemes can still be applied successfully by splitting the domain into a coarse-mesh and a fine-mesh region, then explicit time stepping in the coarse-mesh region is combined with local implicit or explicit time stepping in the fine-mesh region. A fully explicit scheme can be found in [13, 12].
On the other hand, Trefftz methods, originating from , incorporate properties of the PDE into the test and trial spaces. This is done by choosing them as (local) solutions of the targeted differential equation. The use of Trefftz spaces allows to reduce the number of degrees of freedom, as compared to the total degree polynomial spaces, however keeping the same accuracy. Work on Trefftz-DG methods for different wave propagation problems includes [15, 3, 4, 7, 16, 19, 6]. In , a Trefftz-DG method in space-time for the second order wave equation is presented, proving -convergence in 1, 2 and 3 space dimensions, as well as -convergence, along with exponential convergence for analytic solutions in 1 space dimension. In 
, Trefftz-DG is applied to the coupled elasto-acoustic system, and well-posedness of the problem, as well as error estimates in mesh-dependent normes, are shown. Both[3, 4]
are formulated for meshes with tensor product structure in time. In a Trefftz-DG method for the acoustic wave equation in first order formulation and in arbitrary space dimension is presented. The formulation works for tensor product (in time) meshes, as well as for tent pitched meshes. Well-posedness and optimal -convergence are proven.
As Trefftz-DG formulations only contain interelement terms, they pose a natural choice to evolve the solution from the bottom to the top of tent elements. In this paper, we focus on the combination of Trefftz-DG formulation and tent pitched meshes, and on its efficient implementation.
This work proceeds as follows. First, we introduce the Trefftz-DG method in Section 2, starting by stating the model problem, defining the Trefftz spaces and finishing the section by formulating the method, as it was introduced in . We continue in Section 3 by reviewing different strategies of discretizing the Trefftz spaces. In Section 4, we discuss some numerical details on how to evolve the solution elementwise on a tent pitched mesh, and in Section 5, we show a way to recover the second order solution from the first order formulation. Finally, we present numerical results, which were obtained by implementation of the method in NGSolve [22, 23], in Section 6.
2 The Trefftz-DG method for the acoustic wave equation
2.1 The acoustic wave equation
Let be a space-time domain in , where
is a Lipschitz bounded domain with outward unit normal vector. Corresponding to Dirichlet and Neumann boundary conditions, we define such that they have disjoint interior, and one of which can be empty. We consider the acoustic wave equation in first order formulation, given by
where we assume that the wavespeed is piecewise constant on .
If the initial condition it the gradient of a scalar field , i.e. , then the first order system is equivalent to the second order system obtained by setting and :
The Laplacian , gradient and divergence are considered with respect to the space variable only.
2.2 Space-time meshes
The mesh of the space-time domain is assumed to consist of non-overlapping Lipschitz polytopes, where , with being the anisotropic diameter defined in (3.2). For each mesh face , for , we assume that it either lies below the characteristic speed , or is verical (parallel to the time axis). In more rigorous terms: Let be the normal vector to with , then either
Notice, however, that no CFL-condition or any other time step size restriction is imposed on the time-like faces.
A mesh with space-like faces only, as the one in Figure 1, is called a tent pitched mesh and can be obtained algorithmically as shown in [11, 26]. It is built by progressively advancing in time, stacking tent-shaped objects on top or each other, each of them union of -simplexes. The main idea is that the tent height is chosen such that the differential equation is explicitly solvable in each tent. Therefore, the local maximal time advance at a spacial point has to respect the causality constraint, which corresponds to a local CFL-condition. This allows to advance the solution tent by tent, not necessarily having to solve a global system. For independent tents, i.e. tents that are not on top of each other, the computations can be done in parallel.
The set of all faces of our mesh,
can be separated into two sets of internal faces
and sets of faces on the boundary of , split up according to their types of inital and boundary conditions
This classification of the faces is represented in Figure 1.
2.3 Trefftz spaces
By definition, Trefftz functions in the kernel of the considered differential operator. For the first order wave equation, we define the local and global Trefftz space as
respectively. Note that, by assuming that the solution is in , we require additional smoothness on the solution, as in general we only have that , for all .
We derive the Trefftz-DG method for any choice of discrete test and trial space with a Trefftz property, which we denote by . A possible choice for a polynomial is given in Section 3 below.
2.4 The Trefftz-DG method
Following , we derive the Trefftz-DG method for the IBVP in (2.1). The method is derived from a local weak formulation, obtained by multiplying the two equations in (2.1) by test and trial functions and , respectively, and integrating by parts on each element of the mesh . Then, adding the two equations gives
By choosing Trefftz test functions , the volume integrals over vanishes. We are left with:
Typical for DG methods, the continuity of the numeric solution on inter-element boundaries is enforced within the bilinear form of the method. To this end, the trace of the solution in the boundary integral has been replaced by the numeric fluxes , which we define below.
To do so, we need to introduce some standard DG notation. For a face shared by two elements , we define the average , jumps in space , and jumps in time , for scalar- and vector-valued functions as follows
where is the unit outer normal vector at split into its space and time components.
Across time-like faces, the information is passed by using centered fluxes, whereas, across space-like faces, the information is passed upward in time, resembling an up-wind scheme. More precisely, the fluxes on the inter-element faces are chosen as
where and are penalty parameters, which will be chosen constant (notice that they are needed on time-like and Dirichlet faces only). By and we denote the trace of the function on space-like faces from the adjacent element at higher and lower times, respectively.
Finally, we plug the definition of the fluxes into (2.4) and sum over all elements . Then the Trefftz-DG method for the wave equation reads:
On a tent pitched mesh, as the one in Figure 1, the method is semi-explicit, meaning that the solution on each tent only depends on the tents below, allowing to solve each tent explicitly, and tents independent from each other in parallel; details are given in Section 4.1 below. The situation where also vertical faces are present, is needed, for instants, in the case of piecewise constant wavespeed, is discussed in Section 4.2 below. Note that the method only includes integrals over element boundaries, thus only quadrature on dimensional simplices is needed.
3 Choice of discrete Trefftz spaces
So far, we have not specified what discretization of the Trefftz space to use. We introduce the straightforward choice, given by all polynomials in space-time that fulfill the first order wave equation. For an element in the mesh , we define the local polynomial Trefftz space as
where we denote by the space of polynomials on of degree . In general, it is possible to choose different polynomial degrees in different elements. Here, we choose a uniform , as this is consistent with the numerical examples below. The global Trefftz-DG space on the whole mesh is then given by The dimension of the elemental Trefftz space is given by
where we recall that for . Notice that, for the total degree polynomial space, one has .
Let us now assume that the first order problem is derived from a second order problem. Then it is natural to derive the vector valued Trefftz space for the first order problem from a scalar Trefftz space for the second order problem. We now detail this approach as it is the one we use for the numerical results presented in Section 6. Let us start by defining the polynomial Trefftz space for the second order problem:
We are able to construct a basis for this space using the recursion formula introduced in [19, Remark 13]. We recall it here, for completeness. We need some multi-index notation: for we denote and . Furthermore, let with 1 in the -th entry. Consider a space-time polynomial
We want to compute the coefficients such that the polynomial is Trefftz. This is done by inserting the polynomial into the second order wave equation and collecting terms of equal power to find that
has to hold for the polynomial to be Trefftz. To start the recursion, we need to choose polynomial bases (in the space variables only) for and , respectively. More precisely, we start by choosing polynomial basis functions for the space and for . Then we can introduce a basis for such that either and , or and for some . Hence, we can construct the basis for out of two sets of polynomial basis functions of and . This lets us determine the dimension as
Then, a Trefftz space for the first order system can be derived from
We have that
It is sufficient to compute the coefficients only once for and then fix the wavenumber by a coordinate transform. Furthermore, for numerical stability, it is convenient to shift the basis functions to the center of the element and scale them by its anisotropic diameter, which is defined by
for a mesh element . For reference coordinates , the coordinate transform given by
transforms the Trefftz basis of wavespeed 1 to Trefftz basis functions of arbitrary wavespeed . In the case of Trefftz functions for the first order system , we need to choose
4 Evolution within a tent
The tent pitched mesh allows to solve local tents explicitly. This is due to the fact that the slope of the mesh faces is below the characteristic speed , thus the local solution on a tent can be computed once the solution on its inflow boundary is known. In Section 4.1, we discuss how to evolve the solution within a tent with constant wavespeed inside the tent itself. The case where the wavespeed changes within a tent is considered in Section 4.2. Notice that, in the constant wavespeed case, tents coincide with mesh elements, while in the latter case tents on the interface contain more than one mesh element.
4.1 Constant wavespeed
Let us denote the bottom and top faces of the tent by and , respectively. Furthermore, tent faces on the boundary are denoted by for Dirichlet and for Neumann boundaries.
Since the solution is explicit on each tent, we only need to solve a local system of size . The system is derived from (2.5) and is given by the following equation
where, in the case , , and in the case , on a given face is the previously computed solution in the tent sharing that face in lower time.
For the numerical integration, we only need an integration rule for -simplices, in order to integrate over the boundary of the tent. We can define an integration rule on the spatial mesh once, which we can then map to the faces of the tent. This idea is visualized in Figure 2. After solving on the tent, we need to evaluate in the integration points on , and store these values for the next tent. On each spatial integration point, we only need to store the most recent results, leading to a total storage of: .
4.2 Piecewise constant wavespeed
Recall that we assume that the wavespeed is constant in time and piecewise constant in space. In this case, we always consider initial spatial meshes that are aligned with the discontinuities of the wavespeed. To treat such a jump within a space-time tent, we need to incorporate the jump terms from our DG formulation (2.5). This involves integrating on the time-like inter-element boundary contained inside the tent, denoted by . According to (2.5), one has to add to the left-hand side of (4.1) the term
Since the tent now includes two mesh elements, the system matrix is now of size . The extension to interfaces between more than two materials follows.
5 Recovery of the solution of the second order equation
In the case where the problem comes from a second order formulation we can substitute and to write the method in terms of test and trial functions from . Then the method (2.5) reads:
Clearly, the constant basis function does not contribute in this formulation. Thus, this formulation produces the same results as (2.5) with . In order to fix the constants and recover the solution to the second order wave equation, we modify the original formulation by adding the additional terms
to the bilinear form , and
to the right hand side , where and . Note that these terms preserve the consistency of the formulation.
Therefore, when evolving the solution inside a single tent, we need to add and to the left- and right-hand side, respectively, of the formulation discussed in Section 4.
6 Numerical results
In this section we present numerical test results in one-, two-, and three-spatial dimensions. The Trefftz-DG method was implemented in NGSolve [22, 23]. If not otherwise stated, we use the following settings for the numerical examples. We consider the problem (2.1) with initial and Dirichlet boundary conditions such that the analytical solution is , where is the standing wave
given here in 3+1 dimensions, and set the wavespeed . The penalty parameters are chosen as . We measure the error
at final time , which we choose at . All timings were performed on a server with two Intel(R) Xeon(R) CPU E5-2687W v4, with 12 cores each.
6.1 Approximation properties of Trefftz spaces
Independently of their combination with tent pitching, Trefftz polynomial spaces possess good approximation properties for wave solutions. In [19, Section 6] , -version approximation estimates for wave solutions in Trefftz polynomial spaces were proven. The derivation of -version approximation estimates in higher dimensions is, to the best of our knowledge, still open. Here, we compare the Trefftz space to a full discontinuous polynomial space. Recall that, thanks to the Trefftz property of the test functions, we were able to cancel the volume integral in the weak formulation (2.4). This does not hold for the full polynomial space. Therefore, we need to add the volume term to the left-hand side of the formulation (2.5), giving the new left-hand side:
We now need to solve for . Notice that, as opposed to the Trefftz-DG method (2.5) the (full) DG method requires the computation of integrals also in space-time volumes.
For the numerical results we have taken the unit square domain, uniformly meshed into space-time squares, and used the Trefftz space or the space of polynomials with maximal order in each component. In the latter case the number of degrees of freedom per element is hence given by .
The results are shown in Figure 3. Both choices exhibit similar, exponential, convergence speed in terms of polynomial degree, although the Trefftz space is only a subset of the polynomials of maciimal degree equal to . The benefits of the Trefftz space becomes clear when comparing errors versus number of degrees of freedom per element, as seen on the right in Figure 3.
6.2 Comparing space-time meshing strategies
In Section 4, we have seen how to advance the solution element wise on a tent pitched mesh. We now compare this approach to solving the full system on a Cartesian (in time) space-time slab. To solve the full system we use a block Jacobi solver. When comparing the timing of the two methods, we consider 4 different cases for the tent pitching approach, first solving the tents sequentially, and then solving them in parallel on 6, 12, and 24 threads. For this comparison, we choose a quasi-uniform mesh of the unit square in space and the final time equal to the mesh size, i.e. one CFL-conforming time step on the Cartesian mesh. For the -version comparison in Figure 4 on the top we fix , and in turn, for the -version comparison on the bottom we fix .
The results, in Figure 4 on the left, show that the error between the two mesh types differs only slightly. On the right in Figure 4, we compare the runtime of solving the full system on a Cartesian mesh, with solving the tents sequentially (on 1 thread), and solving them in parallel. Sequential tent pitching is about one magnitude faster. Moreover, we can investigate the effects of parallelising the computations. Using more threads only gives an advantage for small enough mesh sizes, as we are only able to solve independent tents in parallel.
6.3 Choice of spatial basis functions
As we have seen in Section 3, the recursion formula (3.1) for the derivation of the Trefftz basis functions, can be initialized with an arbitrary choice of polynomial basis functions in space variable only.
In the following, we compare three different choices for the initial polynomial basis functions: monomials, Legendre, and Chebychev polynomials. We compare them in 1+1 dimensions, on the space-time unit square. The mesh considered is the tent pitched mesh shown in Figure 5. The problem is solved globally using formulation (2.5).
The results in Figure 6 show that all choices behave the same for low degrees. However, for higher degrees, Legendre and Chebychev polynomials fail to approximate the solution, due to the bad conditioning of the system matrix, compared to the monomials. The good properties of the two sets of basis functions do not carry over when developed in the recursion.
6.4 Tent pitching in 2 and 3 space dimensions
As discussed in Section 4, we solve elementwise, and in parallel. For this example, we choose as a spatial domain the unit square and the unit cube. The initial quasi-uniform spatial mesh consists of triangles or tetrahedrons of maximal size . We then use tent pitching in and dimensions, until the algorithm stops at time , where we compute the error. The results of this are shown in Figure 7.
In Figure 7 on the left, we plot the error in terms of for different values of polynomial degree of . As typical for DG methods, in the case of a regular enough solution, we observe superconvergence, with the rate , outperforming the expected , see [19, Thm. 6.19]. We also consider convergence in terms of degree of the Trefftz space , and report the results in Figure 7, right plots. For our analytic solution, we can observe exponential convergence.
6.5 Dissipation of energy
For smooth enough functions the energy at a fixed time is given by
on the domain . We solve using the tent pitching algorithm.
The space mesh considered is a uniform partition of the interval into 5 elements. We measure the relative error in the energy given by
In Figure 8 on the left we can see that the error in energy increases over time. The results suggest that the energy of the numerical solution decreases linearly in time. In Figure 8 on the right, we compare the error in the energy at three different times , plotting it against the degree of Trefftz polynomials in a range from 2 to 8. We observe exponential convergence for increasing order. Furthermore, greater times seem to affect the error only by a multiplicative factor.
6.6 Non-uniformly refined spatial meshes
Now that we have verified the convergence of the Trefftz-DG method with tent pitching initialized on quasi-uniform spatial meshes, we test the advantageous of the method on a non-uniformly refined spatial mesh.
|mesh||total #dofs||L2-error||dof-rate||runtime [s]|
In this test, the refinement is applied to resolve a singular solution at the reentrant corner of an L-shaped domain, given by . The mesh refinement strategy used takes the diameter of a spatial mesh elements as
where is the distance of to the reentrant corner, fixing a minimal mesh size of . Motivated by the theoretical results in , we choose .
We consider a model problem with solution given, in polar coordinates, by
where denotes the Bessel function of the first kind. We consider , so that is singular at the origin. We solve up to time for . To avoid numerically integrate the singularity, we use the method to reconstruct the second order solution , introduced in Section 5, to measure the error given by .
The comparison between results obtained with uniform and non-uniform mesh refinement are shown in Figure 9 for Trefftz functions of degree . We compare the two different meshing strategies by plotting them against (global dof). For the uniformly refined meshes, the convergence rate is bounded by the smoothness of the solution , for . We observe a convergence rate of . Using the non-uniformly refined meshes, we are able to recover optimal convergence for the third order Trefftz polynomials, as seen in Figure 9.
Table 1, gives a closer look on some of the properties already visualized in Figure 9 and also shows the runtime (in seconds). For the computations we used 24 threads. In Figure 9 on the bottom right, we compare the run time with the degrees of freedom. We observe that the uniform and the non-uniform mesh take about the same time for comparable numbers of degrees of freedom. Thus, no significant locking, due to the spatial refinement, occurs.
6.7 Wave propagation in heterogeneous material
In the following example we investigate the reflection of a wave at an interface of two different materials. This experimental setup was also perfomed in [15, 2]. We consider the space-time domain , and problem (2.5) with homogeneous Dirichlet boundary conditions.
The wavespeed is the piecewise constant function given by
As initial condition, we take a Gaussian wave given by
where we choose and . The computations are performed with polynomial degree .
Snapshots of the solution are shown in Figure 10. In the Snapshots, the right part of the domain has spatial mesh sizes up to , whereas in the left part we choose as spatial mesh size of , in order to better capture the steeper wavefront in the slower traveling material. First, we see that the initial condition unfolds in the left homogeneous part of the medium. At , the wave crosses over into the material with higher wave velocity. In the next snapshot we can see that the wave splits into a part traveling to the right with a higher velocity and shallow wavefront, and a part reflected at the interface traveling backwards to the left. Finally, at the time , we can also observe the weaker Huygens wave, which traveled parallel to the interface, before traveling back towards the left.
In Figure 11 on the left we present a sketch of the actions described above, also indicating a region where we measured the output
The domain of measurement was chosen , with . The measurement over time is presented in Figure 11 on the right and shows that we are able to distinguish the three incoming waves. We can see the very weak Huygens wave arriving first, followed by the initial wave and the reflected one.
We have presented implementational aspects and numerical results for the Trefftz-DG method for the acoustic wave equation, originally presented in . The implementation in NGSolve was used to solve on Cartesian (in time) meshes and tent pitched meshes in up to 3+1 dimensions, with varying mesh sizes, polynomial degrees, and wavenumber.
The -convergence rates were shown to comply with the analytic results, showing superconvergence for analytic solutions and limited rates in the case of solutions with insufficient regularity. In the latter case, we were able to recover optimal convergence using non-uniform meshes. For analytic solutions we observed exponential convergence rates in the polynomial degree .
Possible developments include the extension to the case of electromagnetic waves (Maxwell’s equations).
This work has been supported by the Austrian Science Fund, grants no. W1245 and F65.
-  T. Apel. Anisotropic Finite Elements: Local Estimates and Applications. Advances in Numerical Mathematics. Teubner, Stuttgart, 1999.
-  W. Bangerth, M. Geiger, and R. Rannacher. Adaptive galerkin finite element methods for the wave equation. Comput. Meth. Appl. Math, pages 3–48, 2010.
-  L. Banjai, E. Georgoulis, and O. Lijoka. A Trefftz polynomial space-time discontinuous Galerkin method for the second order wave equation. SIAM J. Numer. Anal., 55(1):63–86, 2017.
-  H. Barucq, H. Calandra, J. Diaz, and E. Shishenina. Space-time Trefftz-DG approximation for elasto-acoustics. Applicable Analysis, 00:1 – 16, Aug. 2018.
-  W. Dörfler, S. Findeisen, and C. Wieners. Space-time discontinuous Galerkin discretizations for linear first-order hyperbolic evolution systems. Comput. Meth. in Appl. Math., 16:409–428, 2016.
-  H. Egger, F. Kretzschmar, S. M. Schnepp, I. Tsukerman, and T. Weiland. Transparent boundary conditions for a discontinuous galerkin trefftz method. Appl. Math. Comp., 267:42–55, 2015.
-  H. Egger, F. Kretzschmar, S. M. Schnepp, and T. Weiland. A space-time discontinuous Galerkin Trefftz method for time dependent Maxwell’s equations. SIAM J. Sci. Comp., 37(5):B689–B711, 2015.
-  J. Erickson, D. Guoy, J. Sullivan, and A. Üngör. Building space-time meshes over arbitrary spatial domains. Eng. Comput. (Lond.), 20:342–353, 08 2005.
-  R. S. Falk and G. R. Richter. Explicit finite element methods for symmetric hyperbolic systems. SIAM J. Numer. Anal., 36:935–952, 1999.
-  J. Gopalakrishnan, P. Monk, and P. Sepúlveda. A tent pitching scheme motivated by Friedrichs theory. Comp. Math. Appl., 70(5):1114 – 1135, 2015.
-  J. Gopalakrishnan, J. Schöberl, and C. Wintersteiger. Mapped tent pitching schemes for hyperbolic systems. SIAM J. Sci. Comp., 39, 04 2016.
-  M. Grote, M. Mehlin, and S. A. Sauter. Convergence analysis of energy conserving explicit local time-stepping methods for the wave equation. SIAM J. Numer. Anal., 56, 03 2017.
-  M. Grote and T. Mitkova. High-order explicit local time-stepping methods for damped wave equations. J. Comput. Appl. Math., 239, 09 2011.
-  T. J. R. Hughes and G. M. Hulbert. Space-time finite element methods for elastodynamics: Formulations and error estimates. Comput. Methods Appl. Mech. Eng., 66(3):339–363, Feb. 1988.
-  U. Köcher and M. Bause. Variational space—time methods for the wave equation. J. Sci. Comput., 61(2):424–453, Nov. 2014.
-  F. Kretzschmar, A. Moiola, I. Perugia, and S. M. Schnepp. A priori error analysis of space–time Trefftz discontinuous Galerkin methods for wave problems . SIMA J. Num. Anal., 36(4):1599–1635, 12 2015.
-  M. Lilienthal, S. M. Schnepp, and T. Weiland. Non-dissipative space–time hp-discontinuous Galerkin method for the time-dependent maxwell equations. J. Comp. Physics, 275:589 – 607, 2014.
-  R. Lowrie, P. Roe, and B. van Leer. A space-time discontinuous Galerkin method for the time-accurate numerical solution of hyperbolic conservation laws. in Proceedings of the 12th AIAA Computational Fluid Dynamics Conference, no. 95-1658, 06 1995.
-  A. Moiola and I. Perugia. A space–time Trefftz discontinuous Galerkin method for the acoustic wave equation in first-order formulation. Numer. Math., Jul 2017.
-  P. Monk and G. R. Richter. A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media. J. Sci. Comp., 22-23, 06 2005.
-  G. R. Richter. An explicit finite element method for the wave equation. Appl. Numer. Math., 16(1-2):65–80, Dec. 1994.
-  J. Schöberl. C++11 Implementation of Finite Elements in NGSolve. ASC Report 30/2014, Institute for Analysis and Scientific Computing, Vienna University of Technology, 2014.
-  J. Schöberl. NGSolve Finite Element Library. https://ngsolve.org. Accessed: 2019-07-24.
-  O. Steinbach and M. Zank. A Stabilized Space–Time Finite Element Method for the Wave Equation, pages 341–370. Springer International Publishing, Cham, 2019.
-  E. Trefftz. Ein Gegenstück zum Ritzschen Verfahren. Verhandl. 2er Internat. Kongress. Techn. Mechanik Zürich, 1926, 12–17 Sept., pages 131–137, 1926.
-  C. Wintersteiger. Mapped tent pitching method for hyperbolic conservation laws. Diplomarbeit, 2015.
-  L. Yin, A. Acharya, N. Sobh, R. B. Haber, and D. A. Tortorelli. A space-time discontinuous galerkin method for elastodynamic analysis. In B. Cockburn, G. E. Karniadakis, and C.-W. Shu, editors, Discontinuous Galerkin Methods, pages 459–464, Berlin, Heidelberg, 2000. Springer Berlin Heidelberg.
-  A. Üngör and A. Sheffer. Pitching tents in space-time: Mesh generation for discontinuous galerkin method. Int. J. Found. Comput. Sci., 13:201–221, 04 2002.