A Wavelet-Free Approach for Multiresolution-Based Grid Adaptation for Conservation Laws

In recent years the concept of multiresolution-based adaptive discontinuous Galerkin (DG) schemes for hyperbolic conservation laws has been developed. The key idea is to perform a multiresolution analysis of the DG solution using multiwavelets defined on a hierarchy of nested grids. Typically this concept is applied to dyadic grid hierarchies where the explicit construction of the multiwavelets has to be performed only for one reference element. For non-uniform grid hierarchies multiwavelets have to be constructed for each element and, thus, becomes extremely expensive. To overcome this problem a multiresolution analysis is developed that avoids the explicit construction of multiwavelets.


Introduction
Solutions of hyperbolic conservation laws and convection-dominated problems often reveal a heterogeneous structure: in some regions solutions are smooth whereas in other regions strong gradients or even discontinuities occur. In smooth regions a coarse resolution is sufficient to realize a certain target accuracy whereas in non-smooth regions a fine resolution is required. Thus, considering uniform grids for non-smooth solutions results in unnecessary high computational cost. To reduce the computational cost, dynamical local grid adaptation can be performed. The majority of the available adaptive strategies can be categorized into the following three different paradigms: (i) Adaptive methods based on local error estimators. A posteriori error estimates have been developed in [1,14,30,38,[45][46][47]71]. Since the estimators are in general only bounded up to some constant by the error from below but not from above, these are not reliable and efficient. The existence of an upper bound requires a stable variational formulation of the underlying partial differential equation that is still lacking for the weak formulations of nonlinear systems of conservation laws. (ii) Sensor-based adaptive methods. A sensor is used to trigger local grid adaptation.
Examples for "sensors" are gradients, jumps or curvature of the solution, cf. [7,49,[60][61][62]. These methods do usually not provide any control of the error but are used frequently in practice due to their simplicity. (iii) Perturbation methods. The key idea here is to improve the efficiency of a given reference scheme on a uniformly refined reference grid by computing only on a locally refined adapted subgrid, while preserving the accuracy of the discretization on the full uniform grid, cf. [16,17,40,43,48]. This paradigm allows at least some control of the perturbation error between reference and adaptive scheme. In this context, the terms efficiency and reliability are used, however, with different meaning from in the context of a posteriori estimator based grid adaptation, see class (i). The term efficiency is interpreted here as the reduction of the computational cost in comparison to the cost of the reference scheme. The term reliability is used here in the sense of the capability of the adaptation process to maintain the accuracy of the reference scheme.
In this work, we focus on multiresolution-based grid adaptation. This concept belongs to the class of perturbation methods and does not rely on error estimates. To decide where to coarsen the reference grid, a reliable indicator to control local grid refinement is required.
To this end, a multiresolution analysis (MRA) is performed, where the data corresponding to the current solution are represented as data on a coarse level and difference information called details corresponding to data on successive refinement levels. This new representation of the data reveals insight into the local behavior of the solution. It can be shown that the details become locally small with increasing refinement level when the underlying function is locally smooth. As suggested by this so-called cancellation property, we may determine a locally refined grid performing data compression on local details using hard thresholding. This significantly reduces the amount of the data. Based on this thresholding, local grid adaptation is performed, where we refine an element whenever there exist significant details. The core issue of this strategy is to avoid using the fully refined grid at any point of the computation. Moreover, a repetition of time steps as frequently considered in estimator based adaptive schemes is not necessary. The concept of multiresolution-based grid adaptation has originally been developed for finite volume schemes, cf. [25,26,33,40,52,53,56,57,63,64]. It is motivated by the pioneering work of Harten [41][42][43][44] on the construction of cost-effective hybrid finite volume schemes for conservation laws using an MRA. The underlying idea of this adaptation strategy is to perform an MRA of the reference scheme and evolve only significant local contributions in time. Later on, the concept has been extended to discontinuous Galerkin (DG) schemes in [17,48,65]. In [48,65] an adaptive DG scheme for one-dimensional scalar conservation laws has been derived and analyzed. In [35][36][37] it has been generalized to nonlinear systems of conservation laws in multiple space dimensions.
In contrast to this approach, adaptive multiresolution DG methods in the spirit of the work of Alpert et al. [3,15] have been reported in [5,66]. The basic idea in these works is to design an adaptive DG method by representing the numerical solution as well as the 1 3 operators from the partial differential equation in a multiresolution representation. Thus, the evolution is carried out on the multiresolution coefficients. For that reason, their approaches strongly differ from the work of [17,48,65]. For the multiresolution representation multiwavelets are considered in both schemes. Both works can be considered as preliminary, since the presented adaptive strategy is not analyzed and only applied to very few numerical configurations.
A crucial part in the realization of the adaptive scheme is the construction of generators for the contributions of the orthogonal complements containing the difference information, i.e., the details between two successive refinement levels. Previous works, cf. [48,65], are based on the explicit construction of multiwavelets, i.e., orthogonal basis functions for these complement spaces. When considering an MRA on a hierarchy consisting of rectangular elements or triangles this construction can be performed on a single reference element and then the multiwavelets are transferred to the local element by an affine transformation. However, when considering the setting of more general non-uniform hierarchies, e.g., curvilinear grids or stretched grids, the local bases functions for the orthogonal complements have to be constructed independently on each cell in the hierarchy. In principle, this is possible, but computationally expensive.
To overcome this issue we propose an alternative approach to realize the grid adaptation which does not rely on the construction of basis functions for the complements. The key idea thereby is to represent the local contributions from the orthogonal complements in terms of the basis functions for the piecewise polynomial DG space on the next finer level. To this end, we formulate the MRA in terms of projectors to enable a realization independent of the local generators. In previous work the significance of local contributions was characterized by coefficients in the basis expansion of the detail information from orthogonal complements using multiwavelets. Here, we consider a local function norm for the characterization of significance. Thereby the construction of wavelets has become superfluous. This approach enables a very efficient realization of the adaptive concept to arbitrarily shaped elements since the explicit construction of generators for the complement spaces is avoided.
The outline of the paper is thus as follows. In Sect. 2 we briefly summarize our reference DG scheme for the discretization of hyperbolic conservation laws. This scheme is going to be accelerated by incorporating local grid adaptation. The indicator for local grid refinement is based on an MRA applied to the solution of the DG scheme. The underlying concept of the MRA is presented in Sect. 3. The novelty of this work is addressed in Sect. 4 where we discuss several options how to realize the MRA. In particular, we present a wavelet-free strategy that avoids the explicit construction of multiwavelets. The MRA is then applied to the reference DG scheme, see Sect. 5, resulting in the adaptive MR-DG scheme. To compare the influence of the MRA using either the classical approach or the wavelet-free approach we perform computations for well-known benchmark problems in one and two space dimensions. We conclude in Sect. 7 with a summary of our findings.

DG Discretization of Hyperbolic Conservation Laws
For convenience of the reader and to fix some notation we introduce the underlying problem of interest and briefly summarize its discretization by a DG scheme.
Conservation laws. We are interested in the numerical solution of the initial-boundary value problem for hyperbolic systems of conservation laws with the conserved quantities ∶ + × Ω → D defined on the domain Ω ⊂ d , the vector of fluxes Here D ⊂ m denotes the set of admissible states and Ω bc ⊂ Ω is the inflow boundary.
DG discretization. To discretize the problem (1) we apply a Runge-Kutta discontinuous Galerkin (RKDG) scheme following the ideas of Cockburn et al. [20][21][22][23]. The key ingredients are the discretization of the domain Ω , the introduction of a finite-dimensional DG space and the derivation of a weak formulation, Spatial discretization. We partition the domain Ω by G h ∶= {V } ∈I h with a finite number of cells V all of them having a Lipschitz boundary. The corresponding grid is characterized by the index set Following [32], we define a unit normal vector Γ for each face in Γ h such that it coincides with the outward pointing normal vector on Ω bc . In particular, we fix the orientation of the normal vector on the inner faces.
DG Space. On the partition G h of Ω we define a finite-dimensional DG space consisting of piecewise polynomials where Π p−1 (V ) is the space of all polynomials up to total degree p − 1 on the cell V . Thus, functions in S h may contain discontinuities and are not uniquely defined on the boundaries of the cells V , ∈ I h . For reasons of stability and efficiency, the DG space is assumed to be spanned by an orthogonal basis of locally supported functions S h = span ∈I h ,i∈P ,i , where the local degrees of freedom of Π p−1 are enumerated with the index set P.
Weak formulation. To determine an approximate solution in the finite-dimensional space S h instead of solving the infinite dimensional problem (1), we derive a weak formulation following [6,32]. This approach leads to a variational problem defining the semidiscrete DG solution of (1):  [6,32]. Note that the weak formulation (3) is stabilized (1a) t (t, ) + ∇ ⋅ ( (t, )) = 0, (t, ) ∈ + × Ω, Since the basis is not timedependent, the temporal evolution of the DG solution is described by the evolution of the coefficients. Choosing a single basis function ,i for v in (3) and using the orthogonality of the basis, we obtain a system of ordinary differential equations for the coefficients in the basis expansion of the solution (5) where all coefficient vectors ,i and right-hand sides L are comprised in the vectors U ∶= ( ,i ) ∈I h ,i∈P and L(U) ∶= (L( h , ,i )) ∈I h ,i∈P . Then we derive the fully-discrete DG scheme by discretizing and stabilizing the evolution equation for the coefficients (6). For that purpose, we discretize [0, T] by discrete time levels {t n } N n=0 . At each time level the semi-discrete DG solution h (t n , ⋅) is approximated by n h ∈ S m h applying a strong-stabilitypreserving Runge-Kutta (SSP-RK) scheme to (6), cf. [39]. Using an explicit time discretization requires a restriction of the time step size Δt . In case of a purely hyperbolic problem the largest possible time step size is limited by the well-known Courant-Friedrichs-Levy (CFL) condition [20,22]: Δt ⩽ h C CFL ∕C hyp , where the CFL number C CFL depends on the polynomial degree and C hyp is a bound for the spectral radius of the Jacobian of . In the nonlinear case this depends on the current solution and is recomputed in each time step.
To control oscillations near discontinuities, which typically arise in high-order schemes for hyperbolic conservation laws, we have to stabilize the fully discrete DG scheme locally by modifying high-order coefficients where we apply a projection limiter after each stage of the Runge-Kutta scheme. For our computations we choose the limiter by Cockburn et al. [21,23], because it is comparably simple, computationally efficient and effectively suppresses oscillations. It can be applied to quadrilateral grids as well as triangular grids, cf. [23]. Instead of using a limiter for stabilization, artificial viscosity might be used, cf. [8,9,59,76]. However, in numerical tests it turned out that typically this causes spurious oscillations that spoil the efficiency of our multiresolution-based grid adaptation.

Multiresolution Analysis
To accelerate the DG scheme presented in Sect. 2, we want to sparsify the grid G . To trigger grid refinement and grid coarsening we employ the MRA and perform data compression. In contrast to previous works, cf. [35][36][37]48], we outline the concept without specifying a wavelet basis.
Hierarchy of nested grids. The concept of the MRA introduced in [55] is based on a sequence S = {S l } l∈ 0 of nested spaces that are closed linear subspaces of L 2 (Ω) and the union of these spaces is dense in L 2 (Ω) . In the context of RKDG schemes, an appropriate multiresolution sequence can be set up by introducing a hierarchy of nested grids G ∶= {V } ∈I l characterized by the index sets I l .
These are all partitions of the domain Ω , i.e., Ω = ⋃ ∈G l V . The cells of two successive levels are assumed to be nested, i.e., each cell V , ∈ I l , is composed of cells V , ∈ M , where the refinement set M ⊂ I l+1 is characterized by V = ⋃ ∈M V . Obviously, the resolution becomes finer with the increasing refinement level l. For each of these grids we introduce a DG space By definition these spaces are closed subsets of L 2 (Ω) and due to the nestedness of the grids G l they are nested. The union of these spaces is dense in L 2 (Ω) whenever it holds lim We emphasize that the concept is not confined to uniform dyadic hierarchies but non-uniform hierarchies via grid mappings, cf. [34], such as triangulations are possible as well. Note that for a uniform dyadic hierarchy it holds |M | = 2 −d and |V | = 2 −d |V | , ∈ M , ∈ I l .
Multi-scale decomposition. To investigate the difference between two successive refinement levels, we consider the orthogonal complement space W of S with respect to S +1 defined by Since S +1 is finite dimensional, we decompose S +1 into the direct sum of S and its orthogonal complement W , i.e., S +1 = S ⊕ W , ∈ 0 .
Recursively applying this two-scale decomposition yields the multi-scale decomposition of the space S L , L ∈ , into the coarse discretization space S 0 and complement spaces W , 0 ⩽ ⩽ L − 1 : S L = S 0 ⊕ W 0 ⊕ ⋯ ⊕ W L−1 . Thus, each function u ∈ L 2 (Ω) can be represented by an (infinite) multi-scale decomposition u = u 0 + ∑ l∈ 0 d l with its contributions given by the orthogonal projections where P V ∶ L 2 (Ω) → V denotes the L 2 -projection to some closed linear subspace V ⊂ L 2 (Ω) . In particular, it holds u l+1 = u l + d l , l ∈ 0 . Applying this relation recursively we can express u L ≡ P S L (u) for L ∈ 0 in terms of the contributions of the orthogonal complement spaces and the coarsest discretization space S 0 , i.e.,

Separation of local contributions.
Since the spaces S l as well as W l are piecewise polynomials, i.e., discontinuities are present at element interfaces, the orthogonal projections (9) can be computed locally on each element. This allows to spatially separate the local contributions in the multiscale decomposition. For l ∈ 0 the compactly supported local contributions of the projections (9) are defined by where V denotes the characteristic function on cell V . Consequently, we may decompose u and d into the sum of their local contributions, i.e., Based on the infinite multi-scale decomposition (9) we formulate a localized multi-scale decomposition We emphasize that the local contributions d l may become small when the underlying function is locally smooth. To see this, we note that as a consequence of Whitney's theorem, cf. [31], we can estimate the local approximation error by for any u| V ∈ H p (V ) on a convex cell V , ∈ I l , l ∈ 0 . From this result we infer by the Cauchy-Schwartz inequality the cancellation property Thus, the local smoothness implies that the norm of the local contributions d decrease with order p for increasing . This fact allows for data compression and is the key of the multiresolution based grid adaptation. Significance, thresholding and grid adaptation. By neglecting small contributions from the orthogonal complement spaces, we may locally sparsify the multiscale decomposition of u L corresponding to an arbitrary but fixed refinement level L ∈ 0 and define an approximation u L, . For this purpose, we have to find a suitable norm to measure the local significance of d to decide if the contribution can be neglected. First of all, we define for ∈ I , ∈ 0 , a local complement space equipped with a local norm ‖ ⋅ ‖ ∶ W l, → defined by In principle, any equivalent norm can be used giving us some flexibility and also covering the results from previous works [37,48].
By introducing local threshold values ,L ⩾ 0 we distinguish between significant and nonsignificant contributions. To retain flexibility, we require only a weak constraint on the choice of the local thresholds assuming that there exists an max > 0 such that for all L ∈ there holds Then a local contribution For a uniform dyadic hierarchy a straightforward choice of the local thresholds is given by ,L = 2 −L L with L ⩾ 0 as considered in previous works, cf. [35][36][37]48]. If a local contribution d is not ‖ ⋅ ‖ -significant, we discard it and, thus, sparsify the multi-scale decomposition. For this purpose, we define the sparse approximation I is defined as the smallest set containing the indices of ‖ ⋅ ‖ -significant contributions, i.e., and being a tree, i.e., The set D can be determined by adding additional indices of non-significant contributions such that the tree condition (22) is fulfilled. Then, in analogy to classical wavelet analysis, cf. [24], the error introduced by the thresholding procedure can be estimated by for q ∈ {1, 2} using (18). Moreover, by definition of the L 2 -norm the thresholding proce- Since D has a tree structure, we can identify u L, with the orthogonal projection to a piecewise polynomial space corresponding to an adaptive grid. Hence, u L, can be written as where the adaptive grid G = ⋃ ∈I V is characterized by the index set Then, the sparsified representation u L, can be characterized equivalently by either the index set I or the index set D corresponding to the cells of the adaptive grid (25) and the significant local contributions of the orthogonal complements, respectively. The adaptive grid is illustrated in Fig. 1.

Realization of the Multiresolution Analysis
To implement the MRA basis functions or rather generators for the spaces S and W , ∈ 0 are required. For that purpose, we consider the general framework of stable biorthogonal bases for multiresolution representations that has been established in [18,27]. The choice of the basis determines how to compute the projections (9). Thus, it strongly influences the efficiency of the implementation.
Orthogonal bases. In the following, we refer to the basis functions for S as scaling functions and to the basis functions for W as multiwavelets, i.e., with index sets I S and I W corresponding to the global degrees of freedom of S and W , respectively. The local degrees of freedom of S and W on a cell V are enumerated with index sets P and P * , respectively. In particular, the local degrees of freedom of the orthogonal complement spaces can be enumerated by Due to the definition of the spaces there holds Then we encode the spatial and polynomial degrees of freedom, i.e., I S ∶= I × P and I W ∶= I × P * . For reasons of stability it is convenient to consider orthonormal bases, i.e., This setting is a special case of [18,27]. Moreover, for efficiency reasons we require that the basis functions are compactly supported, i.e., For technical reasons, we additionally require that the zeroth scaling function is constant, i.e., Since the multiwavelets form a basis of W , they are orthogonal to the scaling functions We expand the global functions u and d , ∈ 0 in the bases (26), i.e., Since u ≡ P S (u) and d ≡ P W (u) are defined by orthogonal projections of u ∈ L 2 (Ω) , the coefficients are given by Next, we make use of the locality of the bases (28) and expand the local contributions u , u +1 and d as For ease of notation, we comprise all local scaling functions and multiwavelets corresponding to a cell V in vectors, i.e., Additionally, we merge the vectors of scaling functions of all subcells in the refinement Thus, we may rewrite (33) in a more compact way as where the coefficients are comprised in vectors u ∈ |P| , d ∈ |P * | and u M ∈ |P| |M | , respectively. Due to (32) and (34), these vectors are given by where the inner products are applied component-wise.
Based on the nestedness of the spaces, we derive two-scale relations between the scaling functions and multiwavelets of successive refinement levels. To this end, we rewrite , ∈ I S , and , ∈ I W , in terms of scaling functions on level + 1 by Due to the locality of the bases (28), this reduces to where = ( , i) and = ( , j) . Then, we comprise the inner products in (38) in matrices ,0 ∈ |P|×|P| |M | and ,1 ∈ |P * |×|P| |M | , respectively. Thereby, we can write (38) equivalently in matrix-vector notation as Since the bases are assumed to be orthonormal, the matrix is an orthogonal matrix. Thus, is invertible and, consequently, scaling functions on a finer level can be represented by the sum of scaling functions and multiwavelets on the next coarser level, i.e., Furthermore, we deduce from the orthogonality of that where by |M |×|P| ∈ |M | |P|×|M | |P| we denote the identity matrix. Based on (39) and (41) we derive the following relations linking the coefficient vectors (37) of successive refinement levels corresponding to the multi-scale decomposition (36): for ∈ I , ∈ 0 with local transformation matrices ,0 and ,1 defined by (38) and (39). This provides an explicit representation of the local projections (11), i.e., Thus, in the discrete setting of coefficient vectors and basis expansions, cf. (36), these projections can be realized with matrix-vector multiplications. Due to the orthogonality of there holds ‖ ‖ 2 = ‖ T ‖ 2 = 1 . Hence, the two-scale transformation and its inverse, i.e., are well-conditioned.

Construction of multiwavelets.
To determine appropriate bases with the desired properties, i.e., orthonormal and compactly supported, where in particular the construction can be performed efficiently and stable, we first need to construct the scaling functions. These can be easily computed by orthogonalization and a following normalization of a local monomial basis. A well-known example in one spatial dimension is the shifted and normalized Legendre polynomials. Since the multiwavelets are not uniquely defined by the required conditions, cf. [29], several approaches for the construction of multiwavelets, i.e., the bases for the orthogonal complement spaces, are possible.
Algebraic approach. Following the ideas of [2,3,75] leads to an algebraic algorithm for an explicit construction of multiwavelets: in [2] Alpert proposed an algorithm for multiwavelets on uniform dyadic hierarchies in one spatial dimension. Yu et al. applied this idea in [75] for the construction of multiwavelets on a triangle with a uniform subdivision. The idea can easily be extended to the general multi-dimensional case, cf. [37]. The basic steps in this approach are: (i) construction of a local basis ,i , ∈ I and i ∈ P , for S +1 ⧵ S ensuring that W , ⊂ span i∈P * ,i , ∈ I , (ii) contribution of S in the basis functions are filtered by an orthogonalization with respect to the coarser scaling functions, (iii) orthogonalization of the basis functions among each other using a QR-decomposition, e.g., via Householder transformations, and (iv) normalization of the multiwavelets. Examples of multiwavelets are shown in Fig. 2.
Stable completion. Another possibility for the construction of multiwavelets comes from a discrete point of view. The basic idea is to find a completion of , i.e., for a given ,0 find ,1 ∈ |P * |×|P| |M | such that is an orthogonal matrix. Then, due to (39) the multiwavelets can be determined by This idea originates from [18], where a very general setting is considered. It has been carried out in the particular setting of DG spaces for the one-dimensional uniform dyadic hierarchy in [13] and for the one-dimensional non-uniform dyadic case in [4]. In both works an explicit formula for ,1 is derived. The drawback of this idea is that it cannot be easily realized in the multi-dimensional case. To find a completion for the general non-uniform multi-dimensional case is not trivial and it is open whether an explicit formula similar to the one-dimensional case can be derived. Furthermore, the computation of the completions from [4,13] for the one-dimensional case is still costly since it requires the inversion of a matrix.
Wavelet-free approach. With increasing polynomial order p and spatial dimension d the local degrees of W are increasing dramatically, e.g., in case of the dyadic hierarchy we have |P * | = (2 d − 1) |P| . Hence, the algebraic construction or the completion of the matrix separately for all cells in the hierarchy is computationally expensive. Therefore, we are interested in reducing the computational cost for the construction of multiwavelets. Whenever V and its subcells V , ∈ M , can be mapped by the same affine linear transformation to a reference cell V ref and its subcells V ref , respectively, the construction can be performed on the reference element V ref . Then the multiwavelets on V can be obtained by shifts and rescaling of the multiwavelets constructed on the reference element. Typical examples are dyadic grid hierarchies or regular triangulations, cf. [37,75].
For reasons of efficiency it is thus of major interest to avoid the explicit construction of a basis for the complement spaces. To this end, we propose an alternative approach avoiding the explicit use of multiwavelets: since W ⊂ S +1 , we may express d in terms of scaling functions on level + 1 using (39) by where due to (36) and (39) the coefficients are given by Applying (44) and (42) (50) is more costly than the classical approach using (44) with multiwavelets. Since u has to be computed from u M by (43) in the two-scale decomposition, we may compute d M in a more efficient way by Note that the MRA is based merely on the local contributions d and is not influenced by the choice of the local basis functions. Due to (48) and (51), we can represent d without an explicit use of multiwavelets M or the matrix ,1 . Thus, the construction of the multiwavelets can be avoided using (48) with (51). From the orthogonality of , in particular (40), (49) and (47), we conclude that this alternative representation of d is wellconditioned, i.e., ‖d M ‖ 2 ⩽ ‖u M ‖ 2 . Implementation aspects. The computational cost for the realization of the MRA essentially relies on the approach chosen for the representation of the local wavelet spaces. In the following we discuss the classical approach based on an explicit construction of the multiwavelets, cf. (36) and (44), and the alternative wavelet-free approach by (48) and (51) avoiding the costly construction of multiwavelets or rather the completion of the matrix . The starting point for an MRA is some given function u L ∈ S L ⊂ L 2 (Ω) on a fixed refinement level L ∈ 0 . This might come from a numerical scheme or from a projection of an L 2 -function. The task is to compute the multi-scale decomposition (10), i.e., u 0 , d 0 , ⋯ , d L−1 . Due to (9), the projections in the multiscale decomposition can be computed recursively from fine to coarse, cf. Fig. 3.
Consequently, the computation of the multi-scale decomposition consists of repeated two-scale transformations. The two-scale transformation, i.e., the projections to the next coarser spaces, can be performed cell-wise.
The implementation of the multi-scale transformation depends on the basis expansions for the different spaces. The difference between the classical and the waveletfree approach is the choice of the generating set for the orthogonal complement spaces W . In the classical approach the multiwavelets are used, whereas in the wavelet-free approach the scaling functions of S +1 are used.

Algorithm 1 (classic)
Thus, the implementation of the local two-scale transformation in Algorithms 1 and 2 and its inverse in Algorithms 3 and 4, respectively, is different for the two approaches.

Algorithm 4 (wavelet free)
To compare the computational cost and memory requirements of the two approaches we focus on a local two-scale transformation on a single element described by Algorithms 1-4 neglecting the expensive multiwavelet construction in the classical approach. By this comparison the computational cost and memory requirements of the full transformations can be estimated.
In Table 1 the memory requirement, i.e., the number of coefficients for the storage of d , u and u +1 for ∈ I and the number of matrix entries needed for the transformations, are compared. To quantify the overhead in storing more coefficients in the wavelet-free approach, we consider the ratio of the number of coefficients, The reduction of memory requirement resulting from storing less matrix entries in the wavelet-free approach can be quantified by analyzing the ratio of matrix entries, For dyadic hierarchies the ratios (52) and (53) are listed in Table 2 for d = 1, 2, 3.
On the one hand, avoiding the explicit construction of multiwavelets has the advantage that M ,1 is not needed anymore, cf. (51), and therefore |M |-times less matrix entries have to be stored. On the other hand, there is more memory needed for the storage of the coefficients, since d is expanded in the basis of the larger spaces S +1 .
Next, we compare the computational cost for the two-scale transformation and its inverse in Table 3. For that purpose, we compare the number of multiplications and additions for the computation of d , u from u +1 in Algorithms 1 and 2 and of its inverse in Algorithms 3 and 4.
To quantify the difference of the number of operations needed, we consider the ratio of number of total operations (52) ratio coeff ∶= # coefficients for wavelet-free approach # coefficients for classical approach = 1 + 1 2 |M | .
(53) ratio mat ∶= # matrix entries for wavelet-free approach # matrix entries for classical approach = 1 |M | .  In Fig. 4 the ratio of the total number of the operations ratio op is plotted for different spatial dimensions d and polynomial degrees p.
The realization of the multi-scale decomposition using the wavelet-free approach avoids the multiplication with ,1 and T ,1 , respectively. With increasing polynomial degree p and space dimension d the number of operations needed for the transformation is reduced significantly compared to the classical approach using multiwavelets.
In summary, we conclude that the wavelet-free approach not only avoids the costly construction of multiwavelets, but also is less expensive in terms of the computational cost in comparison with the classical approach based on multiwavelets. The only advantage of the classical approach over the wavelet-free approach is that less coefficients in the basis expansions have to be stored. However, in the classical approach more matrix entries, i.e., storage of submatrix ,1 , have to be stored. Consequently, one may benefit from this advantage whenever it is not necessary to store on each cell separately, e.g., in the case of a uniform dyadic hierarchy. Apart from that, the reduction of memory requirements for storing less matrix entries compensates for the overhead of storing more coefficients in the wavelet-free approach. Thus, even in terms of memory requirement, the wavelet-free approach is more efficient than the classical approach. Hence, it seems to be preferable to use the wavelet-free approach for the realization of the MRA.

Adaptive Multiresolution DG Scheme
We now intertwine the DG scheme and the MRA introduced in Sects. 2 and 3, respectively. The key idea is to apply the MRA to the evolution equations of the DG scheme. Then a multiscale decomposition is performed for the resulting evolution equations. We omit all equations for detail coefficients in a cell that is considered to be non-significant. The remaining evolution equations are projected back to an adaptive grid corresponding to the significant cells.
Reference scheme. The starting point for the grid adaptation is a hierarchy of nested grids. Corresponding to the grid on a fixed maximum refinement level L, we consider an MRA of the semi-discrete DG solution L defined by (3). The weak formulations defining the evolution of this solution can be written in a compact form as where the inner product is applied component-wise and L ∶ (S L ) m × S L → m is defined in (4). Note that L (⋅, ⋅) is linear in the second argument. Due to (10) and (11) we can rewrite the semi-discrete DG solutions in terms of the single-scale coefficients of the coarsest mesh and detail coefficients as In the following we refer to this scheme as well as to its fully-discrete counterpart as the reference scheme.
(54) ratio op ∶= total # operation wavelet-free approach total # operation classical approach Significant local contributions. Next, we derive weak formulations for the local contributions. By choosing test functions v ∈ S 0 and w ∈ W , 1 ⩽ ⩽ L − 1 , with compact support and using orthogonality, we conclude that the evolution of 0 , ∈ I 0 , and the evolution of the local contribution , ∈ I , 1 ⩽ ⩽ L − 1 , are specified by Here, the local contributions of the test functions are defined similarly to (11) To reduce the computational cost of the reference scheme (55), we only consider significant contributions in (58). While dealing with systems of equations, we have to account for the fact, that individual quantities might have different orders of magnitude. Thus, we augment the definition of significance (19) with a rescaling factor for vector-valued ∈ S L m : we call a local contribution ‖ ⋅ ‖ -significant iff For this purpose, we proceed analogously to (20) and define a sparsified representation of (t, ) by where evolution equations for the local contributions of the complements are partially neglected. Note that the significance of local contributions may change in time.
Fully discrete scheme. Note that the above ideas can be directly migrated to the fully discrete scheme where at each time step t n we consider only significant contributions in the local multi-scale decomposition (13) of the solution of the fully-discrete scheme n L , i.e., where the associated index set corresponding to significant local contributions D n ≈ D (t n ) approximates the index set of local contributions in the semi-discrete DG formulation. According to (24) and (25) this sparsified representation corresponds to a projection of n L to an adaptive grid. Here we omit the details and refer to [34] instead. Prediction of significant contributions. To determine which local contributions become significant in the next time step, in principle the solution of the reference scheme n+1 is required. However, the computation of the reference scheme has to be avoided. Thus, we have to predict the index set D n+1 of significant contributions for the next time step t n+1 by an index set D n+1 . The prediction of significant contributions at the new time step t n+1 can only be based on available information corresponding to the old time step t n , i.e., the set (57) To ensure the reliability of the prediction, we have to avoid that future significant contributions are lost. For that purpose we want to ensure where D n+1 is as small as possible. Otherwise the prediction results in excessive grid refinement causing unnecessary computational overhead. Thus, the "over-prediction" has significant influence on the efficiency of the adaptive scheme. For scalar hyperbolic conservation laws in one spatial dimension a prediction strategy has been proposed in [48,65]. This strategy is intertwined with the limiting process to ensure the reliability condition (62). For the proof of reliability of this strategy properties of the reference scheme that are only available for the scalar case in one spatial dimension are essential and therefore the extension to the multidimensional case is not trivial. Thus, we modify the idea of Harten's heuristic strategy [43] and define D n+1 as the smallest superset of D n fulfilling the following constraints: 1) significant contributions remain significant, i.e., D n ⊂ D n+1 ; 2) contributions in a local neighborhood of a significant contribution become significant, i.e., if ∈ D n , then {̃∶ Ṽneighbor of V } ⊂ D n+1 ; 3) new discontinuities may develop causing significant contributions on higher refinement levels, i.e., if max This type of strategy has been originally developed for finite volume schemes. Although the reliability of this heuristic strategy has never been proven to hold, it turned out that it gives satisfactory results in the context of finite volume schemes [16,56,57] as well as for DG schemes [19,[35][36][37].
Adaptive MR-DG scheme. By the prediction set D n an adaptive grid characterized by the index set I n is defined at each time step t n according to (25). Associated to this locally refined grid there exists a discretization space S n ⊂ S L . Analogously, we denote by I n the index set characterizing the adaptive grid defined by D n at time step t n . Furthermore, we denote by S n ⊂ S n ⊂ S L the associated discretization space. Thus, the adaptive DG solution n L, ∈ (S n ) m can be written equivalently in two different ways, i.e., A time step of the adaptive scheme consists of three steps.
1) Refinement. The grid is refined via prediction of significant contributions using the multi-scale representation. Consequently, n L, is represented in the basis of the richer space S n+1 ⊃ S n corresponding to the refined grid I n+1 of the next time step t n+1 , i.e., 2) Evolution. The evolution of the DG scheme with an RK scheme including local projection limiting at each Runge-Kutta stage is computed using the single-scale representation (64). The outcome of the evolution is denoted by ̃ n+1 L, . 3) Coarsening. Significant contributions might turn out to be non-significant, i.e., ̃ n+1 L, might contain non-significant contributions. Thus, we apply the thresholding procedure (20) using the rescaled definition of significance in (59) to compute corresponding to a coarser grid. For a detailed representation of the adaptive MR-DG scheme we refer to previous works [34,65], where it is discussed how to efficiently compute the integrals and how to initialize the adaptive grid avoiding the complexity of the fully refined grid. We emphasize that the MRA and the adaptive DG scheme are very much intertwined. In particular, performing the time evolution on the adaptive grid cannot be considered a DG scheme on an unstructured grid because the local numerical flux computation is connected to the numerical fluxes on the finest level by the multi-scale transformation. In the following we highlight two important conceptual issues.
Remark 1 (Choice of threshold value) Most crucial for the performance of the adaptive scheme is the choice of the local threshold values , . For reasons of efficiency and accuracy, it should be chosen neither too small nor too large to avoid over-refinement and under-refinement, respectively. Following previous works for uniform dyadic grid hierarchies, cf. [35][36][37]48], we choose where h and h L denote the (uniform) diameters of the cells on level and L, respectively. To determine the threshold value L > 0 we follow Harten's idea [42] of preserving the order of accuracy of the reference solution on the uniformly refined grid at the final time step t = t N for L → ∞ , i.e., For this purpose a heuristic strategy was developed in [50]: L = C thr h L with a constant C thr and an exponent . For a non-smooth solution we choose = 1 and = p otherwise. Note that a rigorous a-priori strategy was developed in [48] for scalar conservation laws that turns out to be too pessimistic in general. Furthermore, numerous computations verified that in most cases C thr = 1 is a suitable choice. However, if the solution exhibits weak discontinuities, then it is recommended to choose C thr in the order of the strength of this discontinuity. For a more detailed discussion on the choice of the threshold and the extension to non-uniform grid hierarchies we refer to [34].
Remark 2 (Adaptation and limiting) It is well-known that the MRA provides a reliable shock detector, cf. [41,42]. This motivates to combine adaptation and limiting: we add the MRA as an additional indicator for the detection of troubled cells, i.e., the minmod indicator and the limiter are only applied to cells on the finest refinement level, i.e., I n ∩ I L . In the one-dimensional adaptive scheme proposed in [48] limiting and prediction were intertwined: whenever a cell is indicated for limiting, the prediction enforces a refinement of this cell up to the finest level. Thereby, the stability of the scheme considering the additional troubled cell indicator can be ensured. Since our prediction strategy is not intertwined with the limiting process, we cannot prove the reliability of the additional troubled cell indicator. However, it turned out that if the threshold value is chosen carefully this strategy gives satisfactory results in practice, cf. [19,35,37].

Numerical Results
In [19,[35][36][37] numerical results are presented for numerous test cases using the adaptive MR-DG scheme. There, the classical approach (36) has been applied with the local norms based on the explicit use of multiwavelets. Here we focus on the validation of the newly developed wavelet-free approach (48) with the alternative choice for the local norms (17) and compare it to the classical approach using (67). For validation purposes we first compare the classical approach and the wavelet-free approach for classical one-dimensional problems, see Sect. 6.1. Then, for real applications we consider the well-known benchmark problems of a double Mach reflection and a shock-vortex interaction with a boundary layer in Sects. 6.2 and 6.3, respectively.

One-Dimensional Test Cases
To compare the classical approach and the wavelet-free approach we first consider classical one-dimensional test cases: the Sod Riemann problem [69], the Woodward interaction problem of two blast waves [73,74] and the Shu-Osher problem of the interaction between a moving shock and a density sine wave [67]. For each test case we determine the empirical order of accuracy to validate our reference scheme. For the "exact" solution we perform uniform computations on a very fine mesh corresponding to 12 refinement levels, i.e., 40 960 cells, instead of the analytical solution. Problem. The system of nonlinear conservation laws (1a) is determined by the onedimensional compressible Euler equations with the vector of conserved quantities = ( , v, E) T and the flux ( ) = ( v, v 2 + p, v(E + p∕ )) T . These equations express conservation of mass, momentum and total energy density. Here , v, E = e + 1 2 v 2 and e denote the density, velocity, specific total energy and specific internal energy, respectively. The system is closed by the equation of state for the pressure p = e( − 1) for a thermally and calorically perfect gas with the ratio of specific heats = 1.4 for air at standard conditions.
Configuration. The three test cases differ in the initial data, the boundary conditions and the computational domain.
1) Sod problem:  Computations. We perform adaptive computations starting with a grid on level 0 containing N 0 = 10 cells. We perform computations up to t = 0.25 , t = 0.038 and t = 0.8 for the Sod problem, the Woodward problem and the Shu-Osher problem, respectively, with varying maximum refinement levels 5 ⩽ L ⩽ 9 . Since this test case contains discontinuities, we consider the heuristic strategy with = 1 and C thr = 1 for the choice of the threshold value, i.e., L = h L . In Figs. 5, 6 and 7 the solutions (density ) and the adaptive grid after thresholding corresponding to the index set D n are exemplarily shown for   Fig. 10 Interaction of a moving shock with density sine waves at time t = 1.8 : comparison of classical approach and wavelet-free approach with the reference scheme using N 0 = 10 and 5 ⩽ L ⩽ 9 1 3 approaches the adaptive grids are almost the same, although the thresholding is performed using different local norms.
To confirm that the adaptive computations maintain the asymptotic behavior of the reference scheme we compare in Figs. 8a, 9a and 10a the L 1 -error in the density of the adaptive solution with the error of the reference solution. We observe that for Sod's Riemann problem and Woodward's interaction problem the errors of the adaptive solutions are very close to the discretization error of the reference scheme. For the Shu-Osher problem the error of the adaptive solutions is larger as the error of the reference scheme but the asymptotic behavior is the same. This confirms that by the choice of the threshold value the asymptotic convergence behavior of the reference scheme is maintained. This holds for both the classical and the wavelet-free approach.
To examine the efficiency of the classical approach and the wavelet-free approach, we compare the error in the density with the number of cells in the adaptive grids. Since the adaptive grids change dynamically, we consider in Figs. 8b, 9b and 10b the maximum number of cells over all time steps performing computations with different numbers of refinement levels L. We note that to realize a given error tolerance both the adaptive schemes using the classical and the wavelet-free approach require about the same number of grid cells whereas the reference scheme needs significantly more cells. This is also reflected in the CPU times recorded in Tables 4, 5 and 6. We note that for all test cases the adaptive computations outperform the non-adaptive computation. The computations using the wavelet-free approach are only slightly faster than those using the classic approach. Here the potential of the wavelet-free approach does not become evident because for the one-dimensional-dyadic grid hierarchy the multiwavelets can be computed on a reference element once and for all.

Double Mach Reflection
Problem. Again, we consider the compressible Euler equations. In the two-dimensional Configuration. We consider the reflection of a shock making a 60 • angle with a reflecting wall. This test case has been introduced by Woodward and Colella in [74]. Although the exact solution is not known, it is a well-known benchmark test case for numerical schemes for compressible  Fig. 11. This flow is challenging because the flow on the left-hand side of the separating shock is supersonic, whereas the flow on the right-hand side is subsonic. This change of the flow regime is known to cause problems in the numerical simulation, cf. [23,58].
At the lower boundary, i.e., x 2 = 0 we impose a reflecting wall between x 1 = 1∕6 and x 1 = 4 . Since the shock position in the free flow away from the wall is explicitly known, it is used to specify suitable boundary conditions for the flow field at all other boundaries, i.e., From ext we compute a suitable boundary condition bc . The overall configuration is shown in Fig. 11. Discretization. Simulation parameters are the following. We consider quadrilateral Cartesian grids with uniform cells on each level. The subdivision is dyadic, i.e., |M | = 2 d . The computations have been performed using a third-order scheme, i.e., p = 3 . For the time discretization an explicit third-order SSP-RK method with three stages and C CFL = 0.1 is used. To exactly compute the integrals in the scheme we use a tensor product based Gaussian quadrature formula with (p + 1) d points. To avoid limiting in regions, where the solution is smooth, we consider a Shu constant of M = 50 , see [21,23]. To reduce the Computations. We perform adaptive computations starting with a grid on level 0 containing 20 × 5 cells. We perform computations up to t = 0.2 with different maximum refinement levels 5 ⩽ L ⩽ 7 . Since this test case contains discontinuities, we consider the heuristic strategy with = 1 and C thr = 1 for the choice of the threshold value, i.e., L = h L . In Fig. 12 the adaptive solution (density ) and the corresponding adaptive grid are exemplarily shown for L = 6.
First, we find that all relevant structures of the solution are present in our results by comparing our adaptive results using L = 6 with the (uniform) results from other works, for instance [23,51]. Furthermore, we note that the grid is only refined near to discontinuities or in areas where the solution consists of small structures needing a high resolution. In regions with less fluctuations the adaptive grid is relatively coarse.
In Fig. 13 we compare the adaptive solutions and the corresponding grids for different maximum refinement levels L in the blow-up region on the right-hand side of the domain. Thereby we note that the structures are more sharply resolved with increasing resolution. Moreover, instabilities occur at the bottom of the "red triangle" with increasing resolution. It is well-known and observed by many others that these instabilities occur, cf. [23,58]. In fact, the absence of these instabilities is usually related to a too coarse resolution or too diffusive stabilization in the scheme, cf. [51,70]. Thus, it seems to be reasonable that we observe the instabilities only if the maximum refinement level is sufficiently large. Moreover, this shows that the adaptive scheme is capable of resolving these features even if they develop during the computation and are not present from the very beginning.
To assess the efficiency of the adaptive scheme and to bring out the difference between the classical and the wavelet-free approach we list the maximum number of cells in the grid N max ∶= max{|I n | ∶ n = 0, ⋯ , N} for the different computations in Table 7. For  instance, the adaptive grid for L = 6 consists of at most 40 768 and 33 469 cells using the wavelet-free approach and the classical approach, respectively. Again, we observe that the grids using the wavelet-free approach are slightly more refined than the grids using the classical approach. However, the uniform grid of the reference scheme on level L = 6 consists of 409 600 cells. Thus, the adaptive grids contain only approx 9.95% and 8.17% of the cells of the uniform grid of the reference scheme. For all considered maximum refinement levels the difference between both approaches is around 2%~3% points. Compared with the computational complexity of the reference scheme, i.e., N max = |I L | , this difference is of minor relevance, since the computational costs are reduced significantly with both approaches. The difference can be explained by the different local norms (17) and (67) that are equivalent but not identical. Note that due to Parseval's identity, (17) can be computed equivalently by either ‖d ‖ 2 ∕ √ �V � or ‖d M ‖ 2 ∕ √ �V � avoiding the costly computation of ,1 . Finally, we summarize in Table 8 the CPU times for the adaptive computations. We note that wavelet-free computations are faster than the classic approach on lower refinement levels but become more expensive with increasing number of refinement levels. This is caused by the fact that the adaptive grids for the wavelet-free computations become denser with increasing refinement levels in comparison to the classic computations. As can be concluded from Table 7 the ratio of the maximum number of cells N max for the waveletfree approach and the classic approach increases from 1.13% ( L = 4 ) to 1.35% ( L = 7 ). We emphasize that the difference is not significant as was already observed for the 1D test cases. As it was to be expected, on Cartesian grid hierarchies the classic approach is slightly preferable. To demonstrate the benefit of the wavelet-free approach we consider the following test case.

Shock-Vortex Interaction with a Boundary Layer
So far, all computations were performed on dyadic grid hierarchies. We emphasize that for a dyadic grid hierarchy the cost for the explicit construction of multiwavelets can be significantly reduced because it needs to be done only once for a reference cell. Then the multiwavelets are mapped to the local elements by an affine mapping. Thus, the benefit of the wavelet-free approach is limited. However, in case of grid hierarchies with a nonuniform subdivision, e.g., using stretched grids, there does not exist an affine mapping to a reference element and the multiwavelets have to be constructed element wise. Then it is mandatory to apply the wavelet-free approach because for this configuration the computational costs for the classical approach are prohibitively high. For an example we consider the test case of a shock-vortex interaction with a boundary layer using stretched Cartesian grids. For this configuration we do not perform the computation for the classical approach because the construction of the multiwavelets and the evaluation is too costly as can be concluded from Tables 2 and 3 1)) is based on the reference Mach number Ma and the heat capacity ratio . Due to the viscous terms we need to modify the reference DG scheme, see Sect. 5, where we stabilize the viscous fluxes applying the BR-2 scheme by Bassi et al. [10][11][12]. Again, the system is closed by the equation of state for the pressure for a thermally and calorically perfect gas with = 1.4.
Configuration. We investigate the interaction of a shock wave with the boundary layer of an adiabatic wall in two dimensions ( d = 2) which is a well-known test case for compressible Navier-Stokes equations, cf. [28,54,68,72]. For this purpose, we consider a quadratic shock tube [0, 1] 2 bounded by insulated adiabatic walls. The tube is filled with ideal gas at rest. Initially, a membrane is located at x 1 = 0.5 separating two reservoirs containing gas with different densities and pressures. After removing the membrane at t = 0 , a shock wave followed by a weak contact discontinuity and a weak right-moving rarefaction wave is moving to the right. Due to the presence of viscosity, boundary layers develop at the solid walls at x 2 = 0 and x 2 = 1 . When the shock reaches the wall at x 1 = 1 it is reflected. Then, the reflected wave interacts with the boundary layer resulting in complex wave structures.
Due to the symmetry of the problem, we compute the numerical solution only on the lower half of the tube, i.e., Ω = [0, 1] × [0, 0.5] and consider symmetry boundary conditions at x 2 = 0.5 . The initial configuration is shown in Fig. 14 Discretization. To adequately resolve the thin boundary layer developing at the wall y = 0 the grid hierarchy consists of stretched Cartesian meshes determined by applying a dyadic subdivision on the parameter domain Ω of the grid mapping ( ) = (x 1 , 0.5 + tanh( (2 x 2 − 1))∕(2 tanh( ))) T , ∈ Ω where the strength of the stretching is specified by . For tending to zero the stretching coincides with the identity resulting in a uniform dyadic hierarchy, whereas for increasing , the cells are more and more condensed toward x 2 = 0 . Here, we use = 1.8 . The time step is calculated by Δt = min(C visc h 2 , C CFL ∕C hyp h) with C visc = 0.01∕ and C CFL = 0.01 . The computations have been performed on a 16 × Intel Xeon Gold 5118 cluster with 2.3 GHz CPU frequency using 280 nodes.
Computations. We perform adaptive computations up to t = 1 with different maximum refinement levels 3 ⩽ L ⩽ 5 based on a grid consisting of 30 × 15 cells on level 0. Due to the stretching we apply a localized strategy for the threshold value choosing to [34], Sect. 4.5.2. For the choice of the threshold value we again apply the heuristic strategy using C thr = 1 and = 1 . To realize the MRA and grid adaptation, we only consider the wavelet-free approach to avoid the costly cell-wise construction of multiwavelets on the non-uniform hierarchy. The adaptive solutions and the corresponding grids at t = 1 are shown in Fig. 15. First of all, we note that with sufficient resolution, i.e., sufficiently large L, our solutions reveal the same structures as computed by other groups using solvers with uniform grids, cf. [28,68,72]. Hence, we observe that the grid adaptation is capable of refining the relevant structures in the solution. The solution using L = 5 reveals some indications that this flow might be unstable. This is backed up by observations in [28,68].
To access the efficiency of the grid adaptation, we again focus on the maximum number of cells during the computation in the adaptive grid N max provided in Table 9.
Here, we note that the adaptation is capable of reducing the computational cost significantly. For sake of completeness, we record in Table 10 the CPU times for the adaptive   computations. Since the full grids are at least three times the size of the adaptive grids, we refrain from performing the computations on the fully refined grids.

Conclusion
The wavelet-free approach enables a straight forward application of the adaptation concept to general hierarchies of nested grids. This approach requires only the construction of the single-scale basis functions which have to be constructed for the (adaptive) DG scheme anyhow.
The numerical results presented here lead to the conclusion that the wavelet-free approach is a suitable alternative to the classical wavelet-based approach to overcome the issues related to the construction of wavelets. However, the grids resulting from the use of the classical approach are slightly coarser than the ones using the more general waveletfree approach. Thereby, both approaches maintain the accuracy of the reference scheme. For that reason, the classical approach might be a good choice whenever the construction of the multiwavelets can be realized efficiently, e.g., when it is sufficient to construct them only on a single reference element. However, the additional refinement of the wavelet-free approach is comparably small in comparison to the overall gain from the grid adaptation. Thus, in most situations, especially when employing non-dyadic grids, this overhead is negligible. Moreover, the realization of the MRA and the adaptive scheme using the wavelet-free approach is computationally less costly than using the classical approach. Therefore, the wavelet-free approach is the preferable approach in the general case.
Funding Open Access funding enabled and organized by Projekt DEAL.

Compliance with Ethical Standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.