Efficiency improvement of the frequencydomain BEM for rapid transient elastodynamic analysis
Abstract
The frequencydomain fast boundary element method (BEM) combined with the exponential window technique leads to an efficient yet simple method for elastodynamic analysis. In this paper, the efficiency of this method is further enhanced by three strategies. Firstly, we propose to use exponential window with large damping parameter to improve the conditioning of the BEM matrices. Secondly, the frequency domain windowing technique is introduced to alleviate the severe Gibbs oscillations in timedomain responses caused by large damping parameters. Thirdly, a solution extrapolation scheme is applied to obtain better initial guesses for solving the sequential linear systems in the frequency domain. Numerical results of three typical examples with the problem size up to 0.7 million unknowns clearly show that the first and third strategies can significantly reduce the computational time. The second strategy can effectively eliminate the Gibbs oscillations and result in accurate timedomain responses.
Keywords:
Elastodynamics Fast boundary element method Wave propagation Exponential window method Frequency domain windowing∎
1 Introduction
Transient analysis is encountered in almost every field in science and engineering. For most realworld problems, numerical simulation is the only viable approach for obtaining accurate solutions. The boundary element method (BEM) has been wellrecognized as a powerful numerical method to treat wave propagation problems. Its reduced dimensionality, high accuracy and its ability in capturing rapid transitions and steep gradients of the fields have made this method particularly suitable for problems with complex geometry, for example, dynamic analysis of porous solids and fracture analysis. In view of the increasingly complex systems that we have encountered nowadays, the development of efficient BEM for largescale wave propagation simulations has become indispensable.
According to the different approximation solution strategies in time space, the BEM for treating elastodynamic wave propagation problems generally follows two approaches, namely, timedomain approaches and frequencydomain approaches (see e.g. the reviews by Beskos Bde97 () and Costabel CM04 ()). Timedomain approaches can be further classified into timestepping methods and the spacetime integral equation method. In these methods the physical problems are directly solved in the real time domain, thus one can observe the phenomenon as it evolves. However, such methods require an adequate choice of the time step size. An improper time step could lead to instability or numerical damping. For recent development of timedomain methods, see e.g. SA97 (); BS12 (); ADF12 ().
In the frequencydomain approach, the frequencydomain boundary integral equation is established and solved at a series of sampling frequencies. The results are then transformed back into time domain by inverse Laplace transform which can be numerically evaluated by the discretized Fourier transform (DFT) AM87 (); PTB05 (); PGS10 (). This approach is attractive for transient analysis because it is: (1) stable in obtaining timedomain solutions, (2) simple in implementation and (3) suitable for parallel solution. In addition, one notes that the frequencydomain problems by itself already has a broad range of applications such as harmonic wave propagation in solids Bde97 (), air damping on micro resonators DY04 (), heat transfer SST12 (), to name a few.
There are two issues need to be dealt with when computing the transient response using the inverse Fourier transform. The first issue is regarding to the discretization and truncation errors introduced by the DFT MR08 (). The discretization in DFT causes wraparoud (or aliasing) error. For damped systems this error can be removed by adding trailing zeroes to damp out the free vibration before the end of the period. However, this scheme is inefficient for low damped systems, and more importantly, it does not work for nondissipative systems such as the elastic solids considered in this work. To circumvent this problem, the exponential window method (EWM) KR92 (); Humar02 () was employed in the BEM transient elastodynamic analysis PGSG11 (); XYCZ12 (). This method adds an artificial damping to the system to effectively damp out the free vibrations. By properly choosing the damping parameter, denoted by , in the EWM, the response period can be arbitrarily set independent of the actual timedomain response. In addition, it has been found that a large can improve the conditioning of the BEM coefficient matrix. However one cannot arbitrarily increase because it also amplifies the truncation error of DFT and therefore causes unacceptable Gibbs oscillations in the later time response. In this paper, a frequencydomain windowing technique is introduced to eliminate the Gibbs oscillations, and consequently accurate time response can be obtained in all simulation periods. The frequencydomain windowing technique in combine with the EWM leads to a simple but efficient modified Fourier transform approach for BEM elastodynamic transient analysis.
The second issue concerns the computational efficiency. In the frequencydomain approach the transient problem is transformed into independent timeharmonic problems. The computational cost grows almost linearly with . For many practical problems sampling frequencies are often required to obtain accurate timedomain solutions. This implies a huge computational cost because solving one frequencydomain eqation in a threedimension domain using classic BEM is already computationally very costly. To resolve this issue, a series of work has been recently published on the development of efficient integral equation approaches for frequencydomain elastodynamics. Examples include, but are not limited to, fast multipole SBD08 (); TC09 (); CBS08 (), matrix BA10 () and precorrectedFFT (pFFT) accelerated BEM YZY10 (), as well as efficient BEM with adaptive cross approximation MS08 (). In this paper, the pFFT technique is employed to accelerate the frequencydomain BEM analysis XYCZ12 () due to its ease in implementation and high computational efficiency.
A key factor that affects the computational cost of the frequency domain approach is the number of iterations required to solve each linear system. In order to achieve further speedup, a least square extrapolation method is proposed to obtain better initial guesses for the iterative solution procedure.
The basic setting of the frequencydomain approach for BEM transient analysis is briefly recalled in Section 2. The modified Fourier transform method for transient elastodynamic anaylsis is proposed in Section 3. In Section 4 the precorrectedFFT and the extrapolation method for initial guess are presented. The performance of the methods put forward in this paper is tested by three examples in Section 5, and a summary is drawn in Section 6.
2 Frequency domain BEM for transient analysis
2.1 Frequency domain boundary element method
The first step in the frequencydomain approach for transient elastodynamic analysis is to transform the time domain into the frequency domain via Fourier transform. Let denote the region of space occupied by a threedimensional elastic solid with isotropic constitutive properties defined by Lamé constants and , Poisson’s ratio and mass density . In a mixed boundary value problem, the boundary consists of a displacement boundary and a traction boundary such that . The displacement governing equation for linear elasticity in the frequency domain is formulated as
(1)  
where, is the circular frequency; is the displacement in the frequency domain; and denote the Nabla and the Laplace operators, respectively; is the traction operator defining the stressstrain relation based on Hooke’s law; and are prescribed boundary displacement and traction on and , respectively.
The boundary integral equation corresponding to equation (1) is formulated as
(2)  
where, (P.V.) indicates a Cauchy principal value (CPV) of the singular integral; the freeterm is equal to for a smooth boundary at ; and denote the th components of the displacement and traction generated at by a unit point force applied at along the direction , given by
(3) 
(4) 
where,
(5)  
is the unit normal to at directed outwards of ;
and are the respective wavenumbers of S and P elastic waves, so that
The boundary integral equation (2) is numerically solved using a surface discretization with triangular boundary elements. Piecewiseconstant approximation is employed in this work. A standard collocation approach yields a discretized system of the BIE (2) of the form
(6) 
where matrices and with dimensions of by , , consist of boundary integrals of kernel functions and on each element respectively; the vectors and contain the nodal displacement and traction at each element. Note that the contribution of the freeterm has been absorbed into the matrix . By enforcing the boundary conditions, the linear system (6) becomes
(7) 
where, the vector collects the unknown nodal displacement and traction components and can be obtained by solving the above linear system. All the matrix and vectors in (7) are functions of frequency . The computational cost for an iterative solver is . This can be greatly reduced by using the pFFT method in section 4.1.
2.2 Conventional frequencydomain approach
Equation (1) or (2) is solved at a series of discrete frequencies. The resulting displacements and tractions are then transformed back to the time domain by using the inverse DFT to obtain the timedomain response.
Let and be the circular frequency and time resolutions, respectively; and be the time period of the transient response. Given the number of sampling points , one has the basic relations and . Note that the Nyquist frequency should be chosen such that responses corresponding to frequencies higher than are insignificant and can thus be discarded.
The procedure for obtaining the transient responses using the frequencydomain approach can be summarized as:

Determine a frequency resolution and the number of sampling points . needs to be small enough to minimize the loss of useful frequency information.

Sample boundary condition , which could be the given boundary displacement or the traction, and perform DFT to obtain the frequencydomain boundary conditions at frequencies ,
(8) 
Conduct analysis for the first frequencies, that is, , to obtain the frequencydomain response for the first frequencies. The responses corresponding to the last can be determined utilizing the conjugate symmetric property of the DFT as,
where, and denotes the complex conjugate of .

Perform IDFT for to obtain timedomain response , i.e., the desired boundary displacement and/or traction,
(9)
Discretization and truncation errors arise when the forward and inverse Fourier transforms are numerically evaluated using DFT and IDFT in (8) and (9). The discretization error is also known as aliasing error. To eliminate it, the timedomain signal is required to be zero for . The truncation error is due to the truncation of the frequency spectrum. It causes high frequency oscillations near the discontinuities of the signals, known as Gibbs oscillations.
In next sections, techniques for reducing the two errors are introduced. Consequently, a general and efficient modified Fourier transform (MFT) method for elastodynamic transient analysis is obtained. The MFT method itself is not new. It was first developed for analyzing power system transients by Mullineux’s group about four decades ago; see MR08 () and the references therein. In this paper, it is applied in elastodynamic transient analysis.
3 Modified Fourier transform method
In the MFT method, two data processing techniques, one with respect to the timedomain data and another to the frequencydomain data, are introduced to reduce the discretization and truncation errors associated with the DFT and IDFT.
3.1 Discretization error & timedomain damping
As aforementioned, the periodic nature of the DFT requires that the timedomain response tends to zero at the end of a period. This however often cannot be satisfied in the dynamic analysis using the frequency domain approach. A common way to circumvent this problem for a dissipated system is to add trailing zeroes to prolong the duration so that the free vibration can be damped out before the end of the period. Such an approach is clearly not applicable for a nondissipated system. Another more general method suitable for all systems is the exponential window method (EWM) KR92 (); Humar02 (). In this method, a complex frequency shift, amounting to an artificial damping, is introduced to efficiently damp out the free vibrations within the selected period. The actual response is obtained after a proper scaling of the damped response.
Recently, the EWM was used in 2D and 3D BEM elastodynamic transient analysis PGSG11 (); XYCZ12 (). Accurate time domain response can be obtained by a proper choice of the damping parameter . To be specific, let be the displacement in the time domain and its Fourier transform is in (1). In applying the EWM, the displacement is scaled by an exponential window function to obtain a damped displacement , i.e.,
(10) 
Since is a bounded function, the damped response tends to be zero for large . The Fourier transform of is , meaning that can be obtained by simply replacing the frequency in by , where . To obtain the transient response, one needs to solve the problem at a series of complex frequencies . Once the damped response is obtained, the actual response can be obtained as
(11) 
The choice of the damping parameter determines the performance of the EWM. On one hand, it cannot be too large since the factor in (11) acts as an amplifier which magnifies the cumulative errors due to the truncation of the frequency space which will be discussed in section 3.2 and BEM error. On the other hand, given the factor is used to damp the timedomain response, it could be supposed that a high value for is required to sufficiently eliminate the aliasing error. A compromise between the two ends leading to the following ruleofthumb for the selection of KR92 ()
(12) 
This rule was used in the work described in PGSG11 (); XYCZ12 (). It was found in XYCZ12 () that in BEM analysis a choice of the coefficient between could often lead to good timedomain results.
In addition to reducing the aliasing error, another important effect of lies in that a positive can improve the conditioning of the system matrices in (6) and (7). This is a nature result of the Gershgorin circle theorem, because the exponential function reduces the absolute values of the nondiagonal entries of the matrices, thus making the eigenvalues more clustered. Hence, the damping factor actually acts as a perconditioner (called artifical damping perconditioner, ADP, in this paper) whose performance can be enhanced by increasing . As such, a larger is always desirable.
Matrix preconditioning is a crucial problem in BEM analysis. This is particularly true for high frequency and multidomain problems, in which the coefficient matrices tend to be much more illconditioned. Various preconditioners exist (e.g., INYN12 (); CSB12 ()) and can be used. However, it is well known that the performance of many existing preconditioners is generally problemdependent and there always is a tradeoff between their performance and computational cost. On the contrary, the ADP is simple yet efficient, because its implementation requires no extra computation and memory. Whereas, it can largely increase the convergence speed of the iterative solvers; see Section 5.
Meanwhile, the ADP can be used in combination with other preconditioners. For instance, in this work a block diagonal preconditioner is used together with the timedomain damping technique described in this section to solve (7).
3.2 Truncation error & frequencydomain windowing
It is observed that choosing in BEM analysis can cause unacceptable Gibbs oscillations in timedomain responses at later times XYCZ12 (). Here, a frequencydomain data windowing technique is introduced to alleviate the the Gibbs oscillations. Consequently, a larger value of can be set. This implies a considerable improvement in the conditioning of the BEM matrices when compared with ; see evidences in Section 5.
The implementation of this technique is rather simple MR08 (). Before the damped frequencydomain response is transformed back into the time domain using IDFT, it is first multiplied by a suitable data window ; that is, the function is transformed in lieu of itself. The function of the frequencydomain windowing is to reduce the Gibbs oscillations by averaging out the oscillations within one period.
3.3 A summary
To summarize, the timedomain damping in conjunction with the frequencydomain windowing techniques leads to a MFT method for elastodynamic transient analysis. The basic implementation procedure of the MFT method is similar to that described in section 2.2. The distinction lies in that in the MFT method the sampling frequencies are instead of , and the DFT formula (8) and IDFT formula (9) are replaced by
(13) 
and
(14) 
respectively.
4 Rapid solution of linear systems
The main computational cost for the frequencydomain approach is the solution of linear systems (7) for a sequence of sampling frequencies , i.e., the parameterized linear systems
(15) 
where, , and are defined analogously. In this section, two methods are described to considerably reduce the computational cost, namely, a pFFT method for accelerating the BEM matrixvector multiplication, and a frequency extrapolation method for obtaining good initial guesses and thus reducing the total number of iterations in solving the sequential systems.
4.1 PFFT accelerated elastodynamic BEM
The computational cost of a matrixvector multiplication in traditional BEM is , which could be prohibitively expensive for largescale elastodynamic simulations. Within the last two decades, significant progress has been made in the development of accelerated BEM techniques. With the accelerated techniques, the computational complexity, including both the CPU time and memory usage, can be considerably reduced, therefore enabling BEMbased largescale computation.
The first accelerated technique used in elastodynamic BEM is the FMM F98 (), which is further developed by several groups; see, e.g., TC09 (); CBS08 (); TNK03 (). In BA10 () the matrix method is applied in frequencydomain elastodynamic BEM. Recently the pFFT method was used to accelerate elastodynamic BEM analysis by the present authors XYCZ12 (). The pFFT method is very easy to implement and to be parallelized, so a high computational efficiency can be achieved with ease. For more information about the pFFT technique and the details of its implementation, the readers are referred to XYCZ12 ().
Theoretically, for problem domains with small surfacetovolume ratios such as porous solids, optimal order for CPU time and for memory can be achieved by pFFT, which is comparable to the FMM and the matrix method. For other problems, the CPU time of pFFT is between and . In these cases, the FMM and the matrix method would be advantageous. However numerically, due to the small overhead, the actually computational cost of pFFT could be less than that of FMM and the matrix method for problems with up to order ; see XYCZ12 () for comparisons.
4.2 Solution extrapolation method
Iteratively solving each linear system in the sequence (15) needs an initial guess, denoted by . Typically, is independently chosen as zero or one vectors. Here a least square extrapolation method is proposed to determine a better initial guess for the solution of based on the calculated solutions at previous frequencies in the sequence; this approach is termed as the solution extrapolation method (SEM).
The motivation is inspired by the fact that solution depend on the frequency . If is a smooth function of and that the frequency step is sufficiently small, it should be possible to obtain by extrapolation from the previous solutions; that is,
(16) 
where, , . Let . Then the coefficient vector can be obtained by minimizing the square error
(17) 
By computing the QRdecomposition , one gets the least square solution of (17)
Thus, the initial guess is can be computed as
(18) 
The main computational cost for evaluating is the computation of matrixvector products .
Optimal value for the number of preceding solutions in (16) clearly depends on the smoothness of the solutions and thus is problem dependent. For the cases studied in this work, it is found that is optimal. Generally, a large does not improve the accuracy of the initial guess, because the solution is often not smooth over a large frquency range and the frequency step may not be so small as required.
5 Numerical examples
In this section, three examples are selected to validate and demonstrate the performance of the MFT method in Section 3 and the SEM for initial guess in Section 4.2. The first example is a classical benchmark for transient elastodynamic analysis: a prismatic rod subject to a sudden applied uniform pressure at one end. The analytic solution is available and is employed to evaluate the accuracy and the efficiency of the method. The second case is an aluminum plate subject to an impact loading. This example is chosen to test the method for problems with a timedependent arbitrary loading. Finally to test the ability of the method for solving largescale problems, a unit cube containing 64 spherical cavities subject to a sudden applied uniform pressure is simulated.
In the following simulations, constant elements are used. The pFFT method in XYCZ12 () is used to accelerate the BEM analysis. All the linear systems are preconditioned using the block diagonal preconditioner and solved by the generalized minimal residual method (GMRES) SS86 () with convergence tolerance . Calculations are performed on a workstation with a Xeon 3.0GHz CPU and 30GB RAM.
5.1 Example 1: prismatic rod
As the first test case, a prismatic rod with dimensions of as sketched in Figure 1 is simulated. The rod is fixed at its left end and subject to a Heaviside traction = at its right end, where N/m and is the Heaviside function. By setting the Poisson’s ratio of the material to zero, the 3D problem is reduced to a 1D problem, of which the analytical solution is available, see XYCZ12 (). The Young’s modulus and the density of the material are set to be Pa and kg/m, respectively.
5.1.1 Effects of the ADP and the SEM
The boundary of the prismatic rod is partitioned into 12800 triangular elements. The response period and the number of samplings are set to be s and , respectively. To test the effects of the damping coefficient and the SEM on accelerating the convergence of the linear system solver, simulations are carried out under the three cases specified by the first three columns of Table 1. For cases without using the SEM, the initial guesses are set to be zero vectors.
Cases  SEM  Total its  Total time (h)  

Case 1  2  ✗  29237  16.0 
Case 2  4  ✗  21758  14.4 
Case 3  4  ✓,  17562  12.1 
The number of iterations for solving the sequential linear systems (15) corresponding to the first 129 sampling frequencies are plotted in Figure 2. The total number of iterations and the total CPU time for solving all the 129 systems in each of the three cases are listed in the last two columns of Table 1. One can see that by using a larger damping with , a 25% reduction in the total number of iterations can be achieved when compared with the case of . By further using the SEM, the reduction is up to about 40%. However, the reduction of the total CPU time is not as considerable as the number of iterations. This is because the matrixvector multiplication in this problem is relatively fast and the preprocessing time (i.e., time for computing the nearfield interations) occupies a large portion.
5.1.2 Effects of the frequency windowing technique
As described in Section 3.2, the frequency windowing technique can be employed to eliminate the Gibbs oscillations in the timedomain response. Here this is verified using the results of case 3, i.e. , in Table 1. Our numerical experience shows that both the Blackman and Hanning windows obtains almost the same results. Thus only the results of Hanning window are presented here.
Figure 3 and 4 demonstrate both the analytical and simulated displacements at the free end and the tractions at the fixed end as functions of time; the nondimensional time (, with and being the wave speed and length of the rod, respectively) is employed in the plots. It is observed that, before using the frequency windowing technique, both the displacement and traction responses show drastic oscillations in later times. The magnitude of the traction oscillations is about 4 times larger than the largest traction value; see Figure 4 (a). But after the frequency windowing technique is used, the red curves obviously show that the oscillations are effectively removed and the resulting responses agree very well with the analytical solutions in all the simulation period; see Figure 4 (b).
5.1.3 Long time behavior of the method
It is well known that timedomain BEM analysis often suffers from either strong numerical damping or instabilities when the calculation time period is long. Here, the longtime behavior of the present frequency domain method is studied. The simulation conditions are the same as the case 3 in Section 5.1.1; but the response period and the number of samplings are set to be s and , respectively. Thus, the corresponding time step is s and 257 frequency domain BEM analyses are needed.
The computed traction history is shown in Figure 5. Clearly, the proposed method shows a very nice behavior in the entire simulation period. Computed traction responses of similar problems can be found in Figure 12 of BMS12 () and Figure 9 of ADF12 (). It seems that the total numbers of time steps required in their work are much larger than that in the present simulation. For example, the number of time steps for obtaining Figure 9 of ADF12 () is about , whereas in the present method only 257 frequency steps corresponding to 257 time steps are needed.
5.2 Example 2: a plate with a hole subject to an impact loading
The performance of the present method is further demonstrated by a more realistic problem; i.e., an aluminum plate with a hole subject to a vertical impact force as illustrated in Figure 6. The material properties of the aluminum plate are: GPa, kg/m and . For more details of the simulation conditions and the load history, the readers are referred to Section 4.2 of XYCZ12 (). In the present computation, a mesh with 14072 triangular elements is used, s, . The damping coefficient is set to be 4, the SEM with is used. The Hanning window is employed for frequency windowing.
The responses of displacement and strain in the loaddirection at points S (see Figure 6) are computed using the present method and compared with the FEM and/or experiment results in Section 4.2 of XYCZ12 (); see Figure 7. By using the frequency windowing technique the Gibbs oscillations in both displacement and strain responses in the later times are effectively removed. The windowed results coincide very well with the FEM results. The agreement of the computed strain with the experimental measurement is also quite good overall.
5.3 Example 3: a cube with 64 spherical cavities
To demonstrate the performance of the presented method for solving largescale elastodynamic wave propagation problems, a model consisting of a solid cube with 64 identical spherical cavities was built and simulated; see Figure 8. The nonoverlapped spherical cavities are randomly distributed inside the cube and their radius is 0.072 m. The cube is fixed at the bottom surface () and is subject to a uniform step traction at the upper surface (). The load is the same as that in Section 5.1. All other surfaces are traction free. The material is aluminum and the properties are the same as those defined in Section 5.2. This example was also used in XYCZ12 (), but here a finer mesh with elements is used; hence, the total degrees of freedoms is nearly 0.7 million. The response period and the sampling size are set to be s and , respectively. Damping coefficient and the SEM with are used. The Hanning window is employed in frequency windowing.
The simulation takes about 101 hours, and consumes 21 GB memory. The CPU time still seems to be quite large. Further reduction of computational expense can be achieved by employing a more effective preconditioning method and by parallel computing. The time responses of the displacement and traction in direction at four boundary points are plotted in Figure 9. Compared with the 1D solutions of the solid rod presented in Section 5.1, the oscillation patterns of the displacement and traction seem reasonable.
6 Summary
Fast solvers are indispensable for largescale wave propagation analysis. Recently, a fast frequency domain pFFT BEM incorporated with the exponential window technique was developed for solving transient elastodynamic problems XYCZ12 (). It has been shown that this method is stable and general in the sense that: (1) the method can be used to obtain accurate time responses for an arbitrarily chosen time period; (2) the method is applicable to solving problems with damping, either small or large, as long as the damping term is linear and the resulting frequencydomain equation has an equivalent integral formulation
In this work, the computational efficiency of the method in XYCZ12 () is further improved via three approaches:

Introducing large damping via the exponent window function. It can be proved that the damping actually plays the role of preconditioner, denoted as the artificial damping preconditioner (ADP) in this paper. Applying a large damping can in general improve the conditioning of the BEM matrices, thus reducing the total number of iterations required to achieve convergence. Numerical results show that the reduction is around 25% for compared with .

Using frequency windowing technique in Fourier inversion. Exponential function with a large damping coefficient would cause unacceptable oscillations in the later time responses. The oscillations stem from the Gibbs phenomenon which is inherent in expanding nonperiodic functions using Fourier series. The frequency domain windowing technique is an effective way to abate the Gibbs phenomenon and is employed in the Fourier inversion of the frequency domain BEM results. It is found that the wellknown Blackman and Hanning windows perform equally well in removing the Gibbs oscillations in the BEM solutions. The performance of the Hanning window is demonstrated by numerical results in Section 5. The timedomain damping in conjunction with the frequencydomain windowing techniques results in a modified Fourier transform method for elastodynamic wave propagation simulation.

Using solution extrapolation method (SEM) for obtaining good initial guess. In this method, the initial guess for the current linear system is obtained from the solutions of the linear systems corresponding to several previous sampling frequencies via extrapolation. The aim is to reduce the number of iterations for iterative solvers. Numerical results indicate that the reduction can be up to 15% and is generally problem dependent.
The effectiveness of the above strategies is demonstrated by three examples. The largest model has nearly 0.7 million DOFs, and is simulated on a workstation with a Xeon 3.0GHz CPU and 30GB RAM.
Acknowledgements.
JX and LW were supported by the NSFC under Grant 11102154 and 11074201, and the New Teacher Fund for Doctor Station from the Chinese Ministry of Education under Grant 20106102120009. WY was supported by Hong Kong Research Grants Council under Competitive Earmarked Research Grant 621411.References
 (1) Beskos DE. Boundary element methods in dynamic analysis: Part II (19861996). Applied Mechanics Reviews 1997; 50(3):149–197.
 (2) Costabel M. Timedependent problems with the boundary integral method, in: Ency clopedia of Computational Mechanics, Stein E, de Borst R., Hughes TJR (eds), Wiley, 2004.
 (3) Schanz M, Antes H. Application of ‘operational quadrature methods’ in time domain boundary element methods. Meccanica 1997; 32(3):179–186.
 (4) Banjai L, Schanz M. Wave Propagation Problems Treated with Convolution Quadrature and BEM. Lecture Notes in Applied and Computational Mechanics 2012; 63:145–184.
 (5) Aimi A, Diligenti M, Frangi A, Guardasoni C. A stable 3D energetic Galerkin BEM approach for wave propagation interior problems. Engineering Analysis with Boundary Elements 2012; 36(3):1756–1765.
 (6) Ahmad S, Manolis GD. Dynamic analysis of 3D structures by a transformed boundary element method. Comput Mech 1987; 2:185–196.
 (7) Polyzos D, Tsepoura KG, Beskos DE. Transient dynamic analysis of 3D gradient elastic solids by BEM. Computers and Structures 2005; 83:782–792.
 (8) Phan AV, Gray LJ, Salvadori A. Transient analysis of the dynamic stress intensity factors using SGBEM for frequencydomain elastodynamics. Computer Methods in Applied Mechanics and Engineering 2010; 199(4548):3039–3050.
 (9) Ding J, Ye W. A fast integral approach for drag force calculation due to oscillatory slip stokes flows. International Journal for Numerical Methods in Engineering 2004; 60(9):1535–1567.
 (10) Simões I, Simões N, Tadeu A, Reis M, Vasconcellos CAB, Mansur WJ. Experimental validation of a frequency domain BEM model to study 2D and 3D heat transfer by conduction. Engineering Analysis with Boundary Elements 2012; 36(11):1686–1698.
 (11) Moreno P, Ramirez A. Implementation of the Numerical Laplace Transform: A Review. IEEE Transactions on Power Delivery 2008; 23(4):2599–2609.
 (12) Kausel E, Roësset JM. Frequency domain analysis of undamped systems. Journal of Engineering Mechanics 1992; 118(4):721–734.
 (13) Humar JL. Dynamics of structures. 2nd edition. Balkema, 2002.
 (14) Phan AV, Guduru V, Salvadori A, Gray LJ. Frequency domain analysis by the exponential window method and SGBEM for elastodynamics. Computational Mechanics 2011; 48(5):615–630.
 (15) Xiao J, Ye W, Cai Y, Zhang J. Precorrected FFT accelerated BEM for largescale transient elastodynamic analysis using frequencydomain approach. International Journal for Numerical Methods in Engineering 2012; 90(1):116–134.
 (16) Sanz JA, Bonnet M, Domínguez J. Fast multipole method applied to 3D frequency domain elastodynamics. Engineering Analysis with Boundary Elements 2008; 32:787–795.
 (17) Tong MS, Chew WC. Multilevel fast multipole algorithm for elastic wave scattering by large threedimensional objects. Journal of Computational Physics 2009; 228(1):921–932.
 (18) Chaillat S, Bonnet M, Semblat JF. A multilevel fast multipole BEM for 3D elastodynamics in the frequency domain. Comput Methods Appl Mech Eng 2008; 197:4233–4249.
 (19) Benedetti I, Aliabadi MH. A fast hierarchical dual boundary element method for threedimensional elastodynamic crack problems. International Journal for Numerical Methods in Engineering 2010; 84(9):1038–1067.
 (20) Yan ZY, Zhang J, Ye W. Rapid solution of 3D oscillatory elastodynamics using the pFFT accelerated BEM. Engineering Analysis with Boundary Elements 2010; 34(11):956–962.
 (21) Messner M, Schanz M. Adaptive cross approximation in an elastodynamic boundary element formulation. PAMM 2008; 8:10309–10310.
 (22) Isakari H, Niino K, Yoshikawa H, Nishimura N. Calderon’s preconditioning for periodic fast multipole method for elastodynamics in 3D. International Journal for Numerical Methods in Engineering 2012; 90(4):484–505.
 (23) Chaillat S, Semblat J, Bonnet M. A preconditioned 3D multiregion fast multipole solver for seismic wave propagation in complex geometries. Commun Comput Phys 2012; 11(2):594–609.
 (24) Chen LY, Chen JT, Hong HK, Chen CH. Application of Cesaro mean and the Lcurve for the deconvolution problem. Soil Dynamics and Earthquake Engineering 1995; 14(5):361–373.
 (25) Chen JT, Chen CH. Analytical study and numerical experiments for Laplace equation with over specified boundary conditions. Appl. Math. Modelling 1998; 22:703–725.
 (26) Fujiwara H. The fast multipole method for the integral equations of seismic scattering problems. Geophysical Journal International 1998; 133(3):773–782.
 (27) Takahashi T, Nishimura N, Kobayashi S. A fast BIEM for threedimensional elastodynamics in time domain. Engineering Analysis with Boundary Elements 2003; 27:491–506.
 (28) Saad Y, Schultz MH. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci and Stat Comput 1986; 7(3):856–869.
 (29) Banjai L, Messner M, Schanz M. RungeKutta convolution quadrature for the Boundary Element Method. Comput. Methods Appl. Mech. Engrg. 2012; 245246(3):90–101.