Triangular Spectral Element simulation of two-dimensional elastic wave propagation using unstructured triangular grids

A Triangular Spectral Element Method (TSEM) is presented to simulate elastic wave propagation using an unstructured triangulation of the physical domain. TSEM makes use of a variational formulation of elastodynamics based on unstructured straight-sided triangles that allow enhanced ﬂexibility when dealing with complex geometries and velocity structures. The displacement ﬁeld is expanded into a high-order polynomial spectral approximation over each triangular subdomain. Continuity between the subdomains of the triangulation is enforced using a multidimensional Lagrangian interpolation built on a set of Fekete points of the triangle. High-order accuracy is achieved by resorting to an analytical computation of the associated internal product and bilinear forms leading to a non-diagonal mass matrix formulation. Therefore, the time stepping involves the solution of a sparse linear algebraic system even in the explicit case. In this paper the accuracy and the geometrical ﬂexibility of the TSEM is explored. Comparison with classical spectral elements on quadrangular grids shows similar results in terms of accuracy and stability even for long simulations. Surface and interface waves are shown to be accurately modelled even in the case of complex topography with the TSEM. Numerical results are presented for 2-D canonical examples as well as more speciﬁc problems, such as 2-D elastic wave scattering by a cylinder embedded in an homogeneous half-space. They all illustrate the enhanced geometrical ﬂexibility of the TSEM. This clearly suggests the need of further investigations in computational seismology speciﬁcally targeted towards efﬁcient implementations of the TSEM both in the time and the frequency domain.


I N T RO D U C T I O N
Many problems in geophysics need to infer the physical and chemical parameter distributions of the Earth's interior from information provided by seismic wave propagation through complex media. Moreover, numerical simulations of earthquake-induced wave propagation within heterogeneous geological structures have important implications in terms of earthquake risk assessment and of strong motion predictions. Continuous improvements during the last decades, both in terms of acquisition techniques and seismic network density, lead today the possibility of exploiting new information provided by high-resolution multicomponent seismograms, over a wide range of frequencies. The development of new seismic interpretation methods, which take advantage of these high-resolution first-order spatial derivatives. The weak formulation has the nice property to handle naturally both interface continuity and freeboundary conditions, allowing very accurate resolution of evanescent interface and surface waves that are of major interest in seismology. Finally, spectral elements are much easier to implement than global spectral methods on today's generation parallel computers, which involve non-uniform distributed access to the memory.
Previous work with spectral element method (SEM) in computational seismology typically used explicit time stepping and quadrilateral elements, with a clever choice of the polynomial basis for the function expansion based on the Legendre family of Jacobi polynomials, and with a set of non-uniform collocation points defined as the Gauss-Lobatto-Legendre (GLL) cardinal points with the associated quadrature. Moreover, quadrangles allow the use of a tensor product to construct the Lobatto-Legendre cardinal interpolation basis (Bernardi & Maday 1992;Deville et al. 2002) in such a way that the mass matrix remains diagonal with a time-step restriction scaling as O(N −2 ), where N is the order of the polynomial basis, while the derivatives can be calculated using an efficient sum-factorization technique (Orszag 1980) at a cost of only O(E N α ) with α = 3 and α = 4 for 2-D and 3-D problems, respectively, and E the number of elements. Apart from that, this method inherits the nice property of the associated GLL quadrature allowing for the exact integration of polynomial functions of order 2N − 1, with only N + 1 points. Finally, straightforward C 0 continuity of the kinematic field can be insure while using conforming quadrangle meshes thanks to the local conforming maps and the cardinal Lobatto-Legendre interpolation basis.
The need for geometric flexibility is especially important in computational seismology when dealing with complex wave phenomena, such as the scattering by rough topographies of the Earth and sea bottom surfaces, or the seismic response of sedimentary basins with complex structures and fault geometries. Today, such flexibility is hardly achieved when using quadrangles for several reasons. First, when using non-conforming grids, continuity between mismatched elements must be enforced in a variational sense by resorting to the so-called mortar method (Bernardi et al. 1994;Ben Belgacem 1999), which adds both programming complexity and operation costs. Second, when resorting to unstructured meshes, that is meshes referenced using only one index, very few robust grid-generation software based on hexahedra are available today and all of them have difficulties to handle rapid multiscale mesh refinements and de-refinements. Third, the stability of explicit numerical schemes can be strongly affected when using quadrangles that are highly skewed, leading to non-uniform resolution and severe time-step restrictions. Higher geometrical flexibility is, therefore, a strong motivation for exploring the use of triangular elements and unstructured meshes, an extension that has been investigated by hp-finite element methods in these last decades using multivariate Lagrangian interpolation (Chen & Babuska 1995).
Extension of SEMs to unstructured triangular meshes hinges upon the application of more complex polynomial basis specifically tailored to accurately approximate smooth functions defined on the n d -simplex. Drawing on ideas from Orszag's work on spherical harmonics (Orszag 1974), and on orthogonal polynomial basis initially derived by Proriol (1957) and Koorwinder (1975), Dubiner (1991) first overcame these difficulties and obtained a tensor-product-like spectral method for triangles based on Jacobi polynomials. In the construction of such orthogonal polynomial basis, Dubiner makes use of a non-rotational collapsed coordinate system, induced by the mapping of each triangular region to a rectangle, and of the generalization of the direct tensor products by the so-called warped tensor product. Dubiner's basis has been extended to incompressible Navier-Stokes equations in 3-D by Sherwin & Karniadakis (1995a,b); , and to geophysical fluid dynamics problems by Wingate & Boyd (1996).
One possible approach relies on a modal hierarchical approximation that has been extensively investigated by Sherwin & Karniadakis (1995a,b, 1996; . This approach makes use of an almost orthogonal approximation basis, based on the modification of the original Dubiner's basis, which separate interior and boundary contributions, in order to insure C 0 continuity between the triangular subdomains. This so-called modified Dubiner's basis (Dubiner 1991), and its interior-orthogonal variant (Wingate & Boyd 1996), makes use of a triangular polynomial truncation in modes {x n y m , m + n ≤ N }, where N is the global polynomial order of the expansion, together with an efficient warped tensor product and a Gauss-type quadrature formula defined on a squared-grid of points. This allows the backward transform to be still computed in O(E N 4 ) in 3-D, and the forward transformation in O(E N q ), where N q is the order of the quadrature. All the computations are done in the modal space. However in practice, as noted in Hesthaven & Teng (2000), the quadrature formula is oversampled, it uses more grid points than required by the order of the polynomial basis. Even though the added complexity in the expression of the modified basis allows for C 0 continuity between the triangular elements, such a construction still depends on the orientation of the local coordinate system of the contributing elements. Moreover, the partial orthogonality of the basis destroys much of the sparsity of the SEM matrices and leads to a non-diagonal mass matrix (Warburton et al. 1999).
We investigate a second approach based on a nodal approximation Warburton et al. 2000;Hesthaven & Teng 2000;Bos et al. 2001). As in classical SEM, resorting to multivariate Lagrange polynomial interpolation allows straightforward C 0 continuity between triangular elements. It is worth noting here that a modal expansion can always be represented as a Lagrange polynomial expansion of equivalent order, defined on an appropriate set of points, a fact efficiently used in classical SEM on the n d -cube. The central part is however the identification of an appropriate set of nodal points on the n d -simplex that ensures good behaviour of the interpolating polynomials. While, in contrast to the n d -cube, the identification of an optimal set of nodal points remains an open problem in the general case, significant progress has been recently made (Chen & Babuska 1995, 1996Warburton et al. 2000;Hesthaven & Teng 2000;Heinrichs 2005;Blyth & Pozrikidis 2005), that allows for the construction of well-behaved high-order approximations of smooth functions defined on the n dsimplex. Unfortunately, these approximations have no underlying tensor-product construction. The computational complexity, for example, the spatial derivatives computation, scales asymptotically as O(E N 2n d ), which becomes prohibitive for large N. However, for small values of N as used in classical SEM, for example, N ≤ 8, the formulation can still be efficient.
Another important issue is the structure of the associated weight matrix. Classical SEM applications in computational seismology are generally explicit and make efficient use of the discrete orthogonality of the Lagrangian basis on the n d -cube with respect to the inherited GLL discrete inner product. In the case of the Triangular Spectral Element Method (TSEM), such an interesting property can only be preserved by resorting to a quadrature formula using the same grid of nodal points than the one used to define the Lagrangian interpolation basis at the cost of a low-order quadrature. This accuracy loss in the weight matrix has been shown (Komatitsch et al. 2001) to produce unacceptable dispersion errors in the simulation of elastic wave propagation.
Here, we explore the use of the unstructured triangle based SEM introduced in  and Warburton et al. (2000), relying on the use of the Dubiner's polynomials (Koorwinder 1975;Dubiner 1991) and multivariate Lagrangian interpolation based on Fekete points Hesthaven & Teng 2000), for 2-D seismological applications in time domain with special attention to dispersion errors. It is worth mentioning similar type of studies that have been recently conducted independently by Pasquetti & Rapetti (2004) in the context of Helmholtz equations, and Giraldo & Warburton (2005) in geophysical fluid dynamics. The organization of the paper is as follows. After introducing the elastodynamic problem and its variational formulation in Section 2, Section 3 presents the geometrical and functional discretization of the problem within the context of triangular spectral elements, and some of the actual issues of the TSEM formulation are reviewed. The TSEM is first explored in Section 4 by solving canonical examples of 2-D elastic wave propagation allowing for comparisons with analytical solutions and classical SEM solutions. Finally more realistic 2-D numerical examples, such as the scattering by an elastic cylinder in a homogeneous half-space, and the elastic wave propagation in a layered medium are provided in Section 5 to illustrate the enhanced geometrical flexibility of the TSEM. The paper ends up with a discussion of further research directions.

T H E E L A S T O D Y N A M I C P RO B L E M
Let us consider an elastic medium occupying an open domain Ω ⊂ R n d bounded by a smooth boundary ∂Ω withΩ = Ω ∪ ∂Ω, and a time interval of interest I = [0, T ] ⊂ R + . The displacement field is denoted by u :Ω × I → R n d , and the velocity field by v(x, t) =u(x, t) where the dot symbol indicates a differentiation with respect to time.
The Euler-Lagrange equations of elastodynamics may be written as with the initial conditions and, for sake of simplicity, let us consider only the following freesurface boundary condition where ρ(x) denotes the reference mass density. The symmetric Cauchy stress tensor is denoted by σ :Ω × I → S ⊂ R n d ×n d where S is the subspace of symmetric second-order tensors of dimension n d (n d + 1)/2. For linear elastic medium, where the linearized elastic stiffness tensor C(x) is symmetric, with minor and major symmetries, and positive definite restricted to S. The infinitesimal strain tensor is denoted ∈ S and is defined as, Although the method can be applied to general anisotropic media, in the present work only linear isotropic elasticity will be considered, with λ and μ the Lamé constants of the medium.
Finally, in the following f is either a collocated force, where F ∈ R n d is the vector of the force amplitude components; or a generalized body force, that is, it can be expressed as f = div[m], where m ∈ S is the symmetric second-order moment tensor, with classically where M ∈ S, δ(·) is the Dirac delta function, and G(·) the source time function. A variational formulation of the elastodynamic problem can be obtained by introducing the space of admissible kinematic fields (displacement or velocity) fields, as where H 1 (Ω) n d denotes the space of vector fields defined on Ω that are square integrable and have square integrable first-order partial derivatives in space, over the domain Ω. Considering the associated space of kinematic field variations at each time t, δS := {w(x):Ω → R n d | w ∈ H 1 (Ω) n d }, the variational formulation leads to search for u ∈ S t such that ∀w ∈ δS and ∀t ∈ I with the initial conditions and where the inner products refer to standard L 2 and energy inner products: (w, u) = Ω w·u dΩ, A(w, u) = Ω ∇w:C:∇u dΩ.

D I S C R E T I Z AT I O N O F T H E E L A S T O D Y N A M I C P RO B L E M
In the following, for sake of simplicity, we shall restrict the presentation to the 2-D case even though generalization to 3-D is straightforward.

Geometrical discretization
The original physical domainΩ is approximated by a decomposition into n e non-overlapping polygonal subdomainsΩ e as While classical SEM are based on a quadrangular domain decomposition, we shall assume here a domain decomposition based on simple straight-sided subdomainsΩ e allowing for a conforming unstructured triangulation I h ofΩ. It is worth noting here that such a triangulation is consistent with most of the modelling tools used in subsurface geophysics and allows to handle complex geological structures. Straight-sided triangular subdomainsΩ e ∈ I h (Ω) can be linearly transformed, see Fig. 1, into the reference simplex, locally compatible invertible map F e :Ω e → T 2 , is an affine map defined as where v 1 , v 2 , v 3 are the position vectors of the three vertices of Ω e . It is worth noting here that such a map yields to a constant Jacobian J e = | ∂(xy) ∂(rs) | over each subdomain, an essential feature for the remaining.
The variational formulation of the elastodynamic can be restated by simply replacingΩ with I h (Ω), together with the inner products: where u e and w e denotes the restriction of u and w, respectively, tō Ω e , that is u e : = u|Ω e :Ω e × I → R n d .

Piecewise polynomial approximation
Associated with the triangulation I h , let us consider a piecewise polynomial approximation of the kinematic fields with the finite dimensional subspaces where P N (T 2 ) denotes the space of polynomials of degree less or equal to N and defined on the reference simplex T 2 , that is, An appropriate starting point in developing the multidimensional expansion is to construct a system of orthogonal polynomials, with respect to the L 2 inner product, on the reference simplex T 2 . Interest in orthogonal polynomials originates from the fact that they naturally define a complete basis and allow for a straightforward computation of their derivatives. In classical SEM this is readily achieved by resorting to a tensor product construction of one-dimensional Jacobi polynomials p α,β n of degree n, depending of two real parameters α and β with α, β > −1 (Abramowitz & Stegun 1972). These are the only polynomial eigenvalue solutions of a singular Sturm-Liouville problem over , allowing for spectral accuracy as shown by Gottlieb & Orszag (1977). Such a tensor product construction is rooted in the bounded definition of the reference domain ᮀ = × in a Cartesian coordinate system. However, the reference simplex domain T 2 does not naturally inherit that property since the bounds of the Cartesian coordinates (r , s) are dependent on each other.
To develop a suitable tensorial type basis for unstructured triangular regions, a new coordinate system is introduced where local coordinates have independent bounds. This can be achieved while using the so-called collapsed coordinate system , a non-rotational coordinate system resulting from the singular mapping, with like in classical SEM, Such a mapping introduces a non-physical singularity at (−1, 1) much like cylindrical or spherical coordinate systems.

Spectral expansion
Drawing on ideas from Orszag's work on spherical harmonics (Orszag 1974), Dubiner (1991) first constructed an orthogonal polynomial basis in two dimensions that takes advantage of the orthogonality of Jacobi polynomials and which was further extended to 3-D by Sherwin & Karniadakis (1995a), . Dubiner's basis has been shown to be related to the polynomial solution of a separable Sturm-Liouville eigenvalue problem defined on T 2 achieving spectral accuracy (Owens 1998). Dubiner's orthogonal basis is given by where (η, ζ ) are the collapsed coordinates, and by definition p α,β k is the kth-order Jacobi polynomial defined on and orthogonal with respect to L 2 w inner product, with where (·) denotes the Gamma function (Abramowitz & Stegun 1972). The basis may be cast as the product of two polynomials in the (η, ζ ) coordinates where 1 n (η) = p 0,0 n (η) is the Legendre polynomial of degree n, and 2 nm (ζ ) = ( 1−ζ 2 ) n p 2n+1,0 m (ζ ) a polynomial introduced to maintain the orthogonality of the triangular basis. The Jacobi polynomial p α,0 k has zeros that are declustered away from ζ = −1 for increasing α, an important feature for maintaining the linear independence of the expansion and for well-conditioning numerical systems . Such factorization of 2-D functions, as the product of a principal function 1 n and a tensorial function 2 nm different for each principal function, is referred by Dubiner as a warped tensor product and results from the decomposition of the commutative algebra over T 2 into two independent linear subspaces. The warped product construction is an extension of the classical direct tensor product that retains the numerical efficiency of the sum-factorization process. The set of polynomials nm defined a complete orthogonal polynomial basis with respect to the uniform Legendre inner product thanks to the factor ( 1−ζ 2 ) n , with the normalization constant, The cardinality of this set is N = (N + 1)(N + 2)/2, and the truncated spectral expansion P N u at degree N of any function u ∈ L 2 (T 2 ) is given by where {û nm } (n,m)∈M denotes the modal coefficients of the spectral expansion of u in the Dubiner's basis. Taking into account the orthogonality relationship, yields to the relationsû nm = 1 γnm u , nm T 2 , and u − P N u , T 2 = 0, ∀ ∈ P N (T 2 ) which simply express that P N u is the orthogonal projection of u onto P N (T 2 ) with respect to the Legendre inner product.

Lagrangian interpolation
The orthogonality of the previous Dubiner's basis is attractive, but the impact of using such spectral expansion is that all modes are needed to evaluate P N u pointwise. Moreover, for the actual resolution of the elastodynamic problem such an expansion has to be tesselated while maintaining a C 0 continuity between subdomains of the triangulation I h and boundary conditions must be approximated. This is not easily handled using the modal spectral expansion, and requires splitting into interior and boundary modal contributions with the penalty that the modified expansion basis becomes only semiorthogonal and leads to higher computational complexity (Sherwin & Karniadakis 1995a;Wingate & Boyd 1996).
Another possible way is to resort, like in classical SEM, to a Lagrangian interpolation approximation (Bos et al. 2001;Warburton et al. 2000;Hesthaven & Teng 2000) based on a set of nodes defined inside and on the edges of the reference simplex. Function continuity is easily enforced by equating coincident nodal values at the nodes shared between subdomains.
The actual construction of a multivariate Lagrangian interpolation relies on the introduction of a complete polynomial basis. Let I N u ∈ P N (T 2 ) be the multivariate Lagrangian interpolation polynomial of u(r, s) ∈ L 2 (T 2 ), defined by a set of N nodes that we shall discuss later, then using the complete Dubiner's polynomial basis of eq. (22), we can write Using the interpolation condition, leads to the inverse discrete spectral transform ) where in the last term an ordering of the modal coefficientsũ j and the polynomial basis has been introduced for sake of simplicity, using a one-to-one relation j: = j(n, m). The generalized Vandermonde matrix V is defined by the elements V kj = j (r k , s k ). The associated discrete spectral transform which determines the modal coefficients from the prescribed function values u(r k , s k ), is therefore given as, Defining the cardinal Lagrangian function basis Combining eqs (30) and (34), and using eq. (33) leads to the genuine form of the multivariate Lagrangian interpolation functions as Taking into account the completeness property of the Dubiner's basis, any function of L 2 (T 2 ) can be decomposed into the formal series Inserting the spectral expansion of eq. (36) in the right-hand side of the discrete spectral transform (33) yields and, therefore where R N u is given by The aliasing error R N u ∈ P N (T 2 ) due to the interpolation is a low-order contribution that depends on the high-order terms of the spectral expansion of u. The interpolation error can be cast as For any function u ∈ L 2 (T 2 ), the error norm u − I N u T 2 is an upper bound of the truncation error u − P N u T 2 and inherits the spectral decay of Dubiner's expansion modesû k with respect to k.

Interpolation points and inner products
Actually the interpolation error may be important and a critical issue, besides the smoothness of the interpolated function, is the choice of a suitable distribution of nodal points {(r k , s k )} N k=1 in T 2 . For classical SEM, interpolation in the reference quadrangle ᮀ resorts to the tensor-product construction of a Lagrangian interpolation basis of P N ( ) defined by N + 1 GLL nodes in . This yields close to optimal interpolation properties, with points defined on the boundary of , a key feature for ensuring C 0 continuity between subdomains and to approximate essential boundary conditions. Extension to non-tensor-product domains, like the n d -simplex, still involves extensive search for optimal sets of interpolation points (Bos 1983(Bos , 1991Chen & Babuska 1995, 1996Hesthaven 1998;Heinrichs 2005;Blyth & Pozrikidis 2006), minimizing the constant L N appearing in the Lebesgue inequality |u − I N u| ≤ (1 + L N )|u − P N u| which provides a bound of the interpolation error. At this stage, a natural choice is the Fekete nodal set . They can be computed numerically by maximizing the determinant of the Vandermonde matrix V, independently of the choice of the ambient polynomial basis. For quadrangles, Fekete points (Bos 1991;Bos et al. 2001) are unique and they are given by the tensor-product of GLL points. For triangles, when the Vandermonde matrix is not singular, the Fekete points have been conjectured (Bos 1991) to have a maximum number of points on the boundary of the reference simplex. These boundary points were shown (Bos 1991; to be GLL points. This allows to enforce C 0 continuity and to use mixed conforming meshes combining triangles and quadrangles (Komatitsch et al. 2001). Fekete points have to be computed numerically (Bos 1983(Bos , 1991Chen & Babuska 1995; and are known up to degree N = 19. They are not actually optimal points, for example, they do not minimize the Lebesgue constant, but exhibit a good behaviour, for example, L N ≤ N and L N ∝ √ N . Moreover, the Lagrangian functions built on Fekete points are maximum at these points. The numerical quadrature associated with the set of interpolation points plays an essential role for collocated spectral methods. In classical SEM, the set of interpolation points and the associated numerical quadrature are built by a tensor-product construction of the GLL points defined on , which allow to exactly integrate polynomials of P 2N −1 over using N + 1 interpolation points. Together with the tensor-product construction, this allows to define an uniformly equivalent discrete inner product and a discrete orthogonality structure, an essential property that optimize the spectral transforms and lead to the popular collocated form of the SEM with a diagonal mass matrix. Today it is yet unclear, despite some interesting progress (Taylor et al. 2005;Heinrichs 2005;Blyth & Pozrikidis 2005), how to define Gauss-Lobatto type of quadratures for nontensor-product domains in order to extend collocated SEM with a diagonal mass matrix to triangles.
A Newton-Cotes quadrature formula can be associated with the N Fekete points {(r k , s k )} N k=1 . Let k be the Lagrangian polynomials of degree N associated with the Fekete points, then the associated quadrature weights are given by where the orthogonality of the Dubiner's basis has been used and 1 is the uniformly constant polynomial 1 (r , s) = 1, with γ 1 = γ 00 = 2 from eq. (28). By construction, the N Fekete points and the associated weights {ω k } N k=1 give a formula which exactly integrates a polynomial of P N (T 2 ), a minimal result.
In terms of the modal coefficients of the Dubiner's expansion, the L 2 inner product yields where Γ is the diagonal matrix that denotes the Dubiner's metric, In terms of the Lagrangian basis built on the Fekete points, and denoting u i = u(r i , s i ) the function evaluation at the Fekete points, the L 2 inner product yields where u − denotes the vector that gather all the Fekete point values of Finally, the matrix W denotes the inner product metric, for example, W =V −T ΓV −1 . Taking into account eqs (41) and (43), the quadrature weights can also be written as and the discrete inner product associated with the Newton-Cotes quadrature formula yields whereW is the lumped diagonal matrix associated with W . This collocated form used in Komatitsch et al. (2001) is however a minimal result since, using the Newton-Cotes quadrature, such a form is only exact if the integrand belongs to P N . This lack of accuracy of the Newton-Cotes quadrature associated with the Fekete points can have dramatic effect when using a diagonal mass formulation with the loss of the spectral convergence (Pasquetti & Rapetti 2004). It is worth noting here that improved multivariate quadrature points have been recently proposed by Taylor et al. (2005) to increase the strength of the numerical quadrature. However the new quadrature points do not coincide with the Fekete points, defined for interpolation and C 0 continuity, and therefore, do not provide for an extension of the collocated spectral element formulation to triangles but open interesting perspectives for more general curved triangulations and non-uniform material properties within a triangle.

Spectral element discretization
Within the framework of an unstructured triangulation I h (Ω) based on straight-sided triangles, the scalar product can be written as where u he := u h | Ωe denotes the restriction of the discrete approximation u h to the element Ω e , with u he • F e ∈ P N (T 2 ). Assuming piecewise constant material density ρ e within the triangulation, where u k − he denotes the vector that gather the values of u he k at the N Fekete points and M e = ρ e J e W is an where M e denotes the element mass matrix with M e = I 2 ⊗ M e , and I 2 is the 2 × 2 identity matrix. This induces a vector space norm that plays the role of the usual unweighted L 2 norm. It is worth noting here that in contrast with classical SEM, the element mass matrix is only block diagonal but it is exactly computed.
On the other hand, the bilinear form requires the evaluation of Cartesian derivatives. In the reference simplex, derivatives can be readily computed making use of the Dubiner's orthogonal basis as Making use of the affine geometrical mapping F e , the Cartesian derivatives can be readily computed where r ,l , and s ,l refers to the geometrical derivatives constants within eachΩ e . Assuming piecewise constant stiffnesses C e ijkl within the triangulation, the bilinear form is written as The external force term (w h , f h ), in the case of a collocated force (7), follows from eq. (46), while in the case of a generalized body force (8) it follows from eq. (53).
Finally the construction of S h t and δS h is completed by imposing C 0 continuity between the subdomains of I h . One of the major advantages of the multidimensional Lagrangian basis (35) built on Fekete points, is that function continuity is easily enforced by equating nodal values of the functions u k − he shared on the boundary of adjacent subdomains. The associated global numbering results from numbering only once nodes that are counted twice or more times when lying on a vertex or an edge of adjacent subdomains. Let us denoteN , the number of global degrees of freedom, that is, the number of distinct nodes of I h , and U − h ∈ R n dN , the vector form associated with the global numbering. Using the preceding results, the discrete variational form of the elastodynamics problem leads to a system of ordinary algebraic differential equations which can be written in standard form as follows, where M denotes the (n dN ) × (n dN ) global mass matrix, and F ext , F int the source and internal force vectors, respectively. The global mass matrix and vector forces are assembled from the element contributions as where A ne e=1 denotes the classical assembling operator, the action of which is the summation of the entries of coincident nodal values, and where

Time discretization
The semi-discrete momentum equation is simply enforced at time t n+1 using a classical Newmark stepping scheme (Hughes 1987;Zienkewicz & Taylor 1989), denotes the global displacement, velocity and acceleration fields at time t n+1 , respectively. The only member of the classical Newmark family that preserves exactly conservation of total angular momentum is the explicit central difference method obtained for β = 0 and γ = 1/2. This scheme is used hereafter for all the numerical simulations. Such an explicit scheme is second-order conditionally stable (Hughes 1987). The Courant condition to be satisfied for time steps is where t the time step, V p is the P wave velocity and x is the minimum distance between interpolation points. The stability value c h depends on the mesh geometry. Experience has shown that a value of c h = 0.5 is good enough to keep the simulations on unstructured meshes stable. In contrast with classical SEM, the resolution phase requires now solving the linear system of eq. (58) since the analytical mass matrix M is no more diagonal. Although this matrix is very sparse (≈1 per cent of non-zero entries), it may have full bandwidth depending on mesh connectivity. Since the goal of this study was only to explore the accuracy and the flexibility of unstructured TSEM, no attempt was made here to optimize the conditioning of the global mass matrix. For the iterative resolution of the linear system, a simple preconditioned conjugate gradient method with a Jacobi pre-conditioner was used. In all the numerical examples presented in this work, 5 to 10 iterations of the conjugate gradient algorithm were found satisfactory, for example, with a relative residual error of 10 −4 , even for long time evolutions. Clearly more evolved approaches remain yet to be explored in the context of elastodynamics, such as the overlapping Schwarz method proposed by Fisher (1997) that may allow less stringent stability criteria for the time stepping.

A comparison with classical SEM
A simple example of 2-D elastic wave propagation in homogeneous medium is first considered to investigate the long term stability and accuracy of the unstructured TSEM and to compare with classical SEM based upon structured quadrangular elements. The dimensions of the problem are shown in Fig. 2, with P wave velocity Vp = 3200 m s −1 , S wave velocity Vs = 1847.5 m s −1 S and a mass density of 2200 kg m −3 . The source is an explosive point source located exactly at the centre of the domain with a Ricker time wavelet of 16 Hz central frequency (f 0 ). The computational domain is subdivided in a structured grid of 40 by 40 quadrangular elements of 50 m side, which are subdivided using the diagonal of each element to create a mesh of triangular elements as shown in Fig. 2. A polynomial order of 5 is used for the interpolation, and therefore, both global grids, for example, the Gauss-Lobatto grid for quadrangular elements and the Fekete grid for triangular elements, have the same  number of 40401 interpolation points. In Fig. 2 both set of interpolation points are shown to coincide at the element edges. The number of points per minimum wavelength λ min = Vs/ f max with f max = 2.5 f 0 is between 4 and 5, which is the accuracy lower limit when using classical SEM (Seriani et al. 1992;Komatitsch & Vilotte 1998). Bi-periodic boundary conditions are imposed to propagate waves several times through the grid and a time step t = 0.5 ms is used, corresponding to a Courant number of 0.3. Simulations have been performed on the triangular mesh and on the quadrangular mesh, both up to 2 s total time (4000 time steps). Comparison of the horizontal displacement predicted by the TSEM and the classical SEM for the horizontal component at receiver 1 can be seen in Fig. 3.
Comparison with the analytical solution shows that both methods perform well even for waves travelling several times through the computational mesh, for example arrivals later than 0.5 s.

Unstructured triangulation
The main advantage of triangular elements is to resort to unstructured meshes for the discretization of complex geometries. In this section we assess the performance of the TSEM when using an unstructured mesh built with the help of the mesh generator Triangle of Shewchuk (1996). It is worth noting here that for homogeneous media, this mesh generator does not seem to be optimal, for example, it does not generate equilateral triangles within an homogeneous medium for bi-periodic boundary conditions. The numerical example of the previous section is repeated here but now using a unstructured mesh shown in Fig. 4. This mesh is composed of 3202 triangles and 40426 global grid points, when using polynomial order of degree 5 for the interpolation within each triangular element (quite similar to the 3200 elements and 40 401 global grid points used in the structured case). Bi-periodic boundary conditions are again imposed. Keeping the same source and time parameters as in the structured case, leads to 5 points per minimum wavelength and a Courant number of about 0.3. A snapshot of the displacement field at t = 0.5 s is shown Fig. 4 where no S wave is observed even in the case of such an unstructured mesh.
As in the previous example, the horizontal component of the displacement field calculated by the TSEM using the unstructured mesh is shown in Fig. 5 against the analytical solution. The result is quite accurate and it compares well with both structured mesh simulations shown in Fig. 3.
In a second test, both the consistent and the lumped mass TSEM are compared when using the unstructured triangular mesh. The results are shown in Fig. 6 where only the first P arrivals are compared to the analytical solution. The consistent TSEM formulation clearly exhibits much better accuracy than the lumped mass formulation.    clear ringing instability after the P wave arrival can be observed in the latter case, a fact that was already noticed by Komatitsch et al. (2001). Nevertheless it is worth mentioning that as expected, the computational time is greatly reduced by almost a factor of 5 when using the diagonal mass matrix formulation.

The Lamb's problem
The Lamb's problem (Lamb 1904) is a classical problem for assessing wave propagation algorithms. It involves elastic wave propagation in an homogeneous elastic half-space due to a point source. When the source is located just below the free surface, three main phases are excited: a direct P wave, a direct S wave and a nondispersive Rayleigh wave. The physical domain considered here is a homogeneous half-space with a P wave velocity of 3200 m s −1 , a S wave velocity of 1847.5 m s −1 and a density of 2000 kg m −3 . The numerical model has a width of 4000 m and a depth of 2000 m. The vertical point source is located just below the surface, at a computational node located at x s = 1500 m and z s = 1950 m, in order to generate a strong Rayleigh wave. The source time function is a Ricker wavelet of 10 Hz central frequency leading to a minimum wavelength of λ min = 74 m. In order to test the spectral convergence of the TSEM, we have used three different meshes composed of N el = 628, 2500 and 6882 triangles, with average triangle edges h of 180 m, 105 m and 65 m, respectively. The polynomial degree in each triangle is chosen from N = 3, 5, 7 and 9. Therefore, the spatial resolution goes from 2.4 points per minimum wavelength (for the coarsest mesh with the lowest polynomial degree) up to around 9 points per minimum wavelength (for the finest mesh with the highest polynomial degree). The snapshots of the displacement field at 0.5 s over each computational mesh are shown in Fig. 7. In order to minimize time discretization errors, we fix the Courant number at c h = 0.2 and calculate the time step for each simulation using eq. (61) and the minimum distance between nodal points in the corresponding mesh.
The analytical solution for this problem is calculated using the Cagniard-de Hoop method (de Hoop 1960) for two receivers located at the free surface in (2200, 2000 m) receiver 1, and (2700, 2000 m), receiver 2. It is worth mentioning that the accuracy obtained from the available code to calculate the analytical solution is around 10 −3 . In Fig. 8 the comparison between the analytical solution and the TSEM simulations for the mesh with N el = 2500 and N = 5, 7 and 9 are shown. It is clear the rapid convergence to the analytical solution while augmenting the polynomial interpolation degree. This illustrates the high accuracy of the TSEM for simulating Rayleigh wave propagation even in the case of an unstructured mesh, as shown in the zoomed image of Fig. 7.
In the following, the error is calculated as an average of the L ∞ error in the two components of receivers 1 and 2, that is, the maximum difference between analytical and numerical seismograms in the time window. Fig. 9 shows the semi-log plot of the error as a function of the polynomial degree N for the three different meshes. As expected, the convergence to the analytical solution is exponential in N for this example, confirming that the overall error is still dominated by the spatial discretization.
Let suppose we ask the method to simulate Rayleigh waves for 1 s with an error level around 2 per cent. From Fig. 9, two ( Both have similar sampling of the wavefield (around 5.6 points per minimum wavelength) and ts (around 0.3 ms). The difference lies on the number of degrees of freedomN in each simulation, a factor directly related to the computational cost of the algorithm. For (h, N ) = (105, 7) we obtain 61713, while for (h, N ) = (65, 5) we obtain 86531. Therefore, the pair with higher interpolation order in the coarser mesh should be chosen.

2 -D N U M E R I C A L S I M U L AT I O N S
More realistic 2-D simulations are now considered to illustrate the ability of the TSEM to take into account complex geometries using unstructured triangulations.

The elastic wedge problem
The problem considered is the reflection and transmission of a Rayleigh wave at a step of the free surface. This so-called 'elastic wedge problem' (Lapwood 1961;Hudson & Knopoff 1964) can be associated to a fault or a stiff cliff at the free surface. The physical domain is defined here as a heterogeneous quadrangle of 2000 m side separated by one of its diagonals into two different elastic media of different materials properties shown in Table 1  illustrates two well-known problems in the simulation of elastic wave propagation: accurate approximation of interface waves; and strong velocity contrasts between elastic media. The physical problem is, therefore, the impinging of a Rayleigh wave into a homogeneous elastic corner, corresponding seismograms can be seen in Fig. 11, for the homogeneous corner, and in Fig. 12, for the heterogeneous corner. When the waves hit the corners, complicated patterns of diffracted waves are generated. In the homogeneous (left) corner, it is essentially characterized by a P wave conversion, a reflected Rayleigh wave travelling backwards and a transmitted Rayleigh wave travelling downwards with an inverted phase as predicted by theory (Lapwood 1961;Hudson & Knopoff 1964), and numerically simulated using classical SEM in Komatitsch & Vilotte (1998). It is worth noting here that the Rayleigh wave is non-dispersive as expected and propagates with the right velocity. At the heterogeneous (right) corner, the diffraction pattern is more complicated due to strong reverberations generated at the wedge and to the dispersive characteristic of the Rayleigh wave in heterogeneous medium. The converted P wave, the transmitted Rayleigh wave, of much lower amplitude than in the homogeneous corner as expected when Material 2 is stiffer than Material 1, and the reflected dispersive Rayleigh wave, can be identified.

Elastic scattering by a cylinder in a homogeneous half-space
This example illustrates the scattering of a plane wave by an elastic cylinder buried in a homogeneous half-space. Such an example appears frequently in non-destructive testing and near surface seismic studies (cavities, gas inclusions). An analytical solution can be obtained from the expansion of the displacement field in series of Bessel functions, as detailed in the Appendix, and used to assess the accuracy of the TSEM built on an unstructured triangulation. From the numerical point of view, the main difficulty lies on meshing the curved internal interface, specially in the case of two media with high velocity contrast where both an accurate approximation of the body and the interface waves are mandatory. In order to illustrate the flexibility of an unstructured triangulation, three cases are considered here: first, a stiff cylindrical inclusion is included within a more compliant elastic space; second, a compliant cylindrical inclusion is embedded in a stiffer elastic space; and finally a more compliant inclusion is included in a stiffer elastic space in order to test the accuracy of the TSEM for extremely high velocity contrasts. In the first case the physical model is a circular inclusion of 500 m diameter with Vp = 3000 m s −1 , Vs = 1732 m s −1 and ρ = 2 kg m −3 , embedded in a more compliant elastic half-space of 2000 m by 2000 m with Vp = 2000 m s −1 , Vs = 1155 m s −1 and ρ = 1 kg m −3 . These parameters were borrowed from Priolo et al. (1994) who studied a two quarter space problem, that is two elastic half-spaces in contact along a vertical material discontinuity, with special emphasis in the simulation of the interface wave travelling along the vertical interface, a geometry that is well suited for classical SEM based on a quadrangular mesh. Here the scattering by an elastic cylinder of a plane wave travelling downwards is considered. In this problem, an interface wave is expected to travel along the circular boundary of the buried cylinder. The triangulation in this example is composed of 3402 triangles, with a total number of 42 846 global Fekete points when using a polynomial degree of order 5 for the interpolation. Using unstructured triangulation, with a maximum area criteria, uniform sampling of the wavefield can be easily achieved both within the elastic cylinder and the homogeneous space. Periodic boundary conditions in the left and right boundaries allow to simulate an incoming horizontal plane wave travelling downwards. The plane wave, with a P polarization and a Ricker wavelet modulation of 16 Hz central frequency, is prescribed by an initial displacement, velocity and acceleration fields. Along the curved interface, the sampling is about 4.5 interpolation points per minimum S wavelength. The wave field is propagated during 0.6 s using a time step of 0.2 ms, leading to a Courant number of the order of 0.2. Snapshots of the displacement field, at t = 0.04 s, 0.24 s and 0.44 s, are shown in Fig. 13. The incoming P plane wave travelling downwards, the PP reflection generated at the cylinder and the strong interface wave can be identified in the snapshot at t = 0.24 s. These events can also be seen in Fig. 14, where the displacement seismograms recorded at receivers located all along the interface are displayed. The P wave arrival can be identified with a kind of sinusoidal moveout, as well as the interface wave that exhibits a linear moveout as long as the receivers are equispaced along the circular interface. It is worth noting here the lack of numerical artefacts even though the circular interface is only approximated by straight-sided triangular elements. Accuracy of the TSEM for this problem, can be assessed when comparing with the analytical solution at receivers on the circular interface located at (1000 m, 750 m) for receiver 1 and (1237.8 m, 922.4 m) for receiver 11. These receivers are displayed in the rightmost snapshot of Fig. 13. The normalized analytical solution for each component of the displacement field and residuals computed using the numerical solution, amplified by a factor of 10 are displayed in Fig. 15. The L ∞ error is computed and shown on top of each seismogram. An error of less than 2 per cent for all the arrivals is observed for both receivers.
In the second case, the physical parameters are interchanged to obtain a more compliant cylindrical inclusion in a stiffer elastic space. The triangulation is now refined within the cylindrical inclusion, where the wave speeds are slower, in order to obtain a uniform sampling of the wavefield within the whole domain. The resulting triangulation is composed of 2370 triangles with a total number of 29 946 global interpolation points using again a polynomial degree of order 5 for the interpolation. The triangulation is shown in the leftmost panel of Fig. 16. The incident wave field is the same as in the previous example with periodic boundary conditions. Same time-stepping parameters as before are used and the total simulation time is 0.65 s. In Fig. 16   solution at the same receivers used for the first case, Fig. 18, ensures an error of less than 2 per cent between the analytical and numerical solutions.
In the third case, a much higher velocity contrast between the compliant cylindrical inclusion and the surrounding space is created by using Vp = 4000 m s −1 , Vs = 2310 m s −1 and ρ = 2 kg m −3 for the elastic space and Vp = 1400 m s −1 , Vs = 600 m s −1 and ρ = 1 kg m −3 for the inclusion. The incident wavefield is the same as in the previous examples with periodic boundary conditions. In order to keep a similar Courant number of 0.2, the time step must be reduced to 0.112 ms for this simulation due to the strong mesh refinement in the cylindrical inclusion. In Fig. 19 snapshots of the displacement field, at t = 0.0 s, 0.2 s and 0.3 s, are shown with the main phases identified in the seismograms of Fig. 20. It is interesting to note in Fig. 21 similar error levels as in the previous case (less than 2 per cent), showing the robustness of the method to deal with such strong velocity contrast within the elastic medium.

Elastic half-space with free-surface and interface topographies
The flexibility of an unstructured triangulation for wave propagation simulation can be illustrated when considering wave propagation in a layered 2-D elastic domain with arbitrary topography of the free  surface and internal material interfaces. In the example considered here, an upper layer is considered as filled by some kind of alluvial deposit made of soft materials. An extensive study of the seismic response of such an alluvial basin is beyond the scope of this paper and the interested reader can be referred to Sánchez-Sesma et al. (1993).
The main interest here is to illustrate the geometrical flexibility and the enhanced mesh adaptivity provided by the use of an unstructured triangulation. The TSEM built on an unstructured triangulation allows enhanced geometrical flexibility as well as the possibility of spatial mesh refinement in a rather natural and easy way. According to the wave velocity variations in the model, the element size can be locally reduced to keep uniform sampling of the wavefield. However, one should keep in mind that for small spectral elements, the distance between interpolation points may become prohibitively small with regard to the Courant number and conditional stability of an explicit time-stepping scheme. A consideration that should also be kept in mind while using highly deformed quadrangles. The spatial triangulation is again obtained here using a maximum area constraint, as provided by the meshing software Triangle, in order to obtain at least five interpolation points per minimum wavelength within the whole domain. Material parameters, as well as some other relevant parameters, are summarized in Table 2. An  explosive point source, with a Ricker time function of 10 Hz central frequency, is located in the upper layer just below the free surface so as to excite a strong Rayleigh wave along the curved topography. The dimension of the numerical domain has a width extension of 2500 m and an average depth of 1600 m. The discretization is made of 2418 triangles, and a polynomial degree of order 5 is used for the interpolation within each triangles. The numerical mesh shown in the upper-left panel of Fig. 22 is artificially extended to the left, right and bottom boundaries in order avoid spurious reflections. This can also be achieved in a more elegant way, using Perfectly Matched Layers as proposed by Komatitsch & Tromp (2003) and Festa & Vilotte (2005). The time step in this simulation is t = 0.25 ms, leading to a Courant number of 0.35 in the upper layer. Typical snapshots of the displacement field are displayed in Fig. 22 together with the unstructured mesh refined towards the upper layer and especially within the soft alluvial basin. The main observed phases here are: a strong Rayleigh wave propagating along the topography, a direct P wave, together with a PP phase generated at the free surface, a converted PS wave at the surface topography, and a strong headwave that propagates downwards into the medium. In the snapshot at t = 0.4 s, bottom-left panel, a Rayleigh to body wave conversion can be observed. Such a conversion is likely to happen when the Rayleigh wave reaches a strong curvature change of the surface topography (Jih et al. 1988;Komatitsch & Vilotte 1998). Within the alluvial valley, the refracted wave (Pr) can be seen travelling along the bottom of the alluvial basin with the P wave velocity of the upper layer, as well as a high amplification of the incoming Rayleigh wave. When considering the line of receivers shown in the upperleft panel of Fig. 22 and located at the surface of the alluvial basin, the previous phases can be identified directly from the displacement seismograms, as shown in Fig. 23.

C O N C L U S I O N S
In this study, a TSEM is investigated within the context of elastic wave propagation in complex geological media. Extension of spectral methods to problems defined on non-tensor-product domains, like the n d -simplex, is an important issue in computational seismology. This would allow expanding the applicability of the SEM to complex geological structures together with improving mesh adaptivity capability. Moreover, using n d -simplex spectral elements allows the use of existing software for automated unstructured mesh generation used by modern geological modellers. For the sake of completeness, we reviewed some of the actual issues on the formulation of such schemes. On n d -simplex subdomains, spectral expansion of the solution can be achieved when resorting to the Dubiner's orthogonal basis that results from a warped tensor product of Jacobi polynomials associated with a local collapsed coordinate system. Such a construction can be seen as a generalization of the classical tensor product construction in use with classical SEM built on tensor product domains. This orthogonal polynomial basis can be shown to be the solution of a Sturm-Liouville eigenvalue problem defined on the reference simplex, and therefore, to achieve spectral accuracy. Based on this orthogonal basis, a TSEM can be constructed by resorting, as in classical SEM, to a multivariate Lagrangian interpolation designed to enforce in a natural way the C 0 continuity requirement between subdomains of the triangulation. Today, even though no optimal set of points is known for interpolation on a non-tensor product-domain, such a multivariate Lagrangian interpolation can be built based on the set of Fekete points that are shown to endow good interpolation properties and to maximize the number of points on the boundary of the subdomains. In these boundaries, Fekete points turn out to coincide with the classical GLL points used in classical SEM. However, in contrast with the GLL set of points defined on tensor product domains, Fekete points do not define an optimal quadrature formula in the sense that the associated Newton-Cotes quadrature formula only integrates exactly polynomials of P N . This is clearly a minimal result that does not allow for the extension of the SEM collocated formulation to non-tensor product domains. As a consequence, the resulting consistent mass matrix is not diagonal, a feature that adds significant computational complexity for explicit time-stepping schemes typically used within the seismological community. Resorting to a lumped mass matrix in TSEM implies a significant accuracy loss that in some cases can even destroy the spectral convergence of the method. Moreover, there is no tensor product structure associated with Fekete points and the computation of spatial derivatives requires some careful analysis in order to avoid a prohibitive O(N 2n d ) asymptotic scaling. As it is, the TSEM is attractive enough to explore its possible use for numerical simulation of elastic wave propagation in seismology. A preliminary exploration is done in this paper through classical canonical examples and relevant 2-D elastic wave problems. Numerical results presented here do show the high accuracy and the geometrical flexibility of the TSEM. Comparison with classical SEM built on quadrangular grids shows that TSEM does share the same level of accuracy and stability. Surface and interface waves are shown to be accurately modelled even in the case of complex topography with the TSEM. Numerical results for specific problems such as 2-D elastic wave scattering by a cylinder embedded in a homogeneous half-space clearly illustrates the enhanced geometrical flexibility and mesh adaptivity of the TSEM and strongly suggests its use in computational seismology.
Although this paper does not address issues on efficient computational implementations of TSEM, nor precise comparisons in terms of computational costs between classical SEM and TSEM, it should be emphasized that without specific improvements, TSEM remains more expensive than classical SEM. Therefore, there is today clearly room for improvements of the TSEM.
Regarding the computational implementation of the TSEM presented here, the actual cost associated with derivative computations may be significantly reduced when resorting to efficient factorization techniques that exploit the symmetries inherent to the Fekete nodal set and a careful ordering of the nodes. This is clearly emphasized in the work of Hesthaven & Teng (2000), where such a careful implementation is shown to lead to a drastic reduction of the derivative computational costs that compare well with those obtained in classical SEM.
Another issue is related to the improvement of efficient strategies that can cope with the sparse global mass matrix resulting from the lack of an accurate quadrature formula associated with the Fekete nodal set. In the time domain, there is some room for improvements mainly with regard to the iterative resolution of the algebraic system that results from an explicit time discretization, therefore, governed by the mass matrix; or a semi-implicit time discretization, governed by a Helmholtz matrix. In this direction efficient iterative solvers using conjugate gradients preconditioned by the additive overlapping Schwarz method, introduced by Dryja & Widlund (1987) and extensively used by Fisher and co-authors (Fisher 1997;Lottes & Fisher 2005) in the context of computational fluid mechanics with SEM, remain to be explored for elastic wave propagation in the time domain, see also Pasquetti et al. (2005). Even though numerical results do suggest stability scaling in O(N −2 ) for the time step of the explicit time stepping in two dimensions, exploring in this context semi-implicit time stepping may be interesting. In the spectral domain, investigating the use of TSEM is of particular interest since a Helmholtz problem has to be solved for each frequency and eventually multiple sources. The geometrical flexibility and the spectral accuracy of the TSEM, in conjunction with an overlapping Schwarz iterative method, appear quite attractive with important implications for exploration geophysics. In both cases, an important related issue is to investigate the conditioning property of the TSEM for elastic wave simulation.

A C K N O W L E D G M E N T S
We are very grateful to Tim Warburton for providing us with his code to calculate Fekete points in the triangle. We would also like to thank Alexandre Fournier, Gaetano Festa and Martin Kaser for enlightening discussions during this work that help us improving the manuscript. The remarks and suggestions of Jacobo Bielak and an anonymous reviewer are greatly acknowledged as well. This work has been supported by the SPICE EEC training network and the ACI "Catastrophes et Risques Naturels" of the French Ministry of Research.

A P P E N D I X A : A N A LY T I C A L S O L U T I O N F O R T H E I N C I D E N C E O F P L A N E P WAV E S U P O N A N E L A S T I C C Y L I N D R I C A L I N C L U S I O N
This solution is well known after the classical work by Jing & Truell (1956). It was summarized later by Pao & Mow (1973) and can be obtained after the expansion of potentials in polar coordinates. Recently Liu et al. (2000) presented the treatment for viscoelastic materials and choose to give explicitly the coefficients of the where E (i)