## 1 Introduction  
Structural solid dynamics is usually formulated either in the time or in the frequency domains [1]. When considering a general linear viscoelastic behavior, momentum balance leads to the semi-discretized form  
\[\mathbf{M}\frac{d^{2}\mathbf{U}(t)}{dt^{2}}+\mathbf{C}\frac{d\mathbf{U}(t)}{dt }+\mathbf{K}\mathbf{U}(t)=\mathbf{F}(t), \tag{1}\]  
where \(\mathbf{M}\), \(\mathbf{C}\), and \(\mathbf{K}\) are, respectively, the mass, damping, and stiffness matrices; \(\mathbf{U}\) the vector that contains the nodal displacements; and \(\mathbf{F}\) the nodal excitations (forces). From now on, we assume the mechanical system described by N nodal degrees of freedom, which gives the size of the different matrices and vectors involved in Equation (1). The direct integration of the previous discrete form, which consists of N second order coupled ordinary differential equations, can be performed using either, explicit or implicit time integration schemes. Addressing fast transient dynamics can be usually accomplished using explicit integrations that require satisfying stability conditions affecting the largest integration time-step \(\Delta t_{max}\), closely related to the size of the elements involved in the mesh \(\mathcal{M}\) covering the domain \(\Omega\). Many robust integration schemas exist and are widely employed, from the very popular Newmark [2] to more advanced schemes able to preserve inherent mechanical properties [3]. In many cases, explicit integrations are combined with mass-lumping to recover a diagonal mass matrix, in which inverse becomes trivial, fact that improves significantly the integration efficiency. In the same spirit, modal analysis offers an appealing alternative route. The last is based on the fact that the properties (positivity and symmetry) of mass and stiffness matrices in Equation (1), imply the existence of a basis \(\{\boldsymbol{\phi}_{1},\ \boldsymbol{\phi}_{2},\ldots,\ \boldsymbol{\phi}_{ \mathbb{N}}\}\) (in which construction will be detailed later), such that the two following relations apply:  
\[\boldsymbol{\phi}_{i}^{\mathbb{T}}\mathbf{M}\boldsymbol{\phi}_{j}=\delta_{ij}, \tag{2}\]  
with \(\delta\) the Kroenecker delta, and  
\[\boldsymbol{\phi}_{i}^{\mathbb{T}}\mathbf{K}\boldsymbol{\phi}_{j}=\kappa_{i} \delta_{ij}. \tag{3}\]  
By defining matrix \(\mathbf{P}\), in which columns correspond to vectors \(\boldsymbol{\phi}_{i}\), the two previous relations can be rewritten in their more compact matrix counterparts  
\[\mathbf{P}^{T}\mathbf{M}\mathbf{P}=\mathbb{M}\equiv\mathbf{I}, \tag{4}\]  
and  
\[\mathbf{P}^{T}\mathbf{K}\mathbf{P}=\mathbb{K}, \tag{5}\]  
\(\mathbf{I}\) is the identity matrix, and \(\mathbb{K}\) is the diagonal matrix with entries \(\mathbb{K}_{ii}=\kappa_{i}\). Thus, if \(\mathbf{C}=\mathbf{0}\) (undamped dynamics) or \(\mathbf{C}=\alpha\mathbf{M}+\beta\mathbf{K}\) (proportional damping), using the coordinates \(\boldsymbol{\varphi}\), \(\mathbf{U}=\mathbf{P}\boldsymbol{\varphi}\), the problem (1) can be rewritten as  
\[\mathbb{M}\frac{d^{2}\boldsymbol{\varphi}(t)}{dt^{2}}+\mathbb{C}\frac{d \boldsymbol{\varphi}(t)}{dt}+\mathbb{K}\boldsymbol{\varphi}(t)=\mathbf{P}^{T }\mathbf{F}(t), \tag{6}\]  
with \(\mathbb{C}\) diagonal, i.e., \(\mathbb{C}_{ij}=\epsilon_{i}\delta_{ij}\) (\(\mathbb{C}=\mathbf{0}\) in the undamped case). Equation (6) constitutes a system of N uncoupled second order ordinary differential equations that can be integrated very efficiently for providing transient responses. Another alternative route consists in considering the frequency instead of the time as coordinate, which is applying the Fourier transform to problem (1), which constitutes a choice particularly appealing in the case of forced regimes and linear models. In the case of nonlinear models, such a route needs appropriate procedures detailed later.