## Basic Principle
### Galerkin Finite Element Method  
Many methods can be used to solve Equation (5), and Ref. [31] uses the time FEM based on the Galerkin weak form to solve Equation (5). Setting \(\mathbf{v}\in H_{v}^{1}\) to be the test function and \(\mathbf{u}^{h}\in H_{E}^{1}\) to be the trial function, the weighted residual weak form of Equation (5) in \([t_{n},t_{n+1}]\) is established as follows:  
\[\int_{t_{n}}^{t_{n+1}}\mathbf{v}^{T}(\mathbf{M}\mathbf{\ddot{u}}^{h}+\mathbf{ \overline{C}}\mathbf{\dot{u}}^{h}+\mathbf{\overline{K}}\mathbf{u}^{h})\mathrm{ d}t=\int_{t_{n}}^{t_{n+1}}\mathbf{v}^{T}\mathbf{f}\mathrm{d}t \tag{6}\]  
where \(H_{E}^{1}\) and \(H_{v}^{1}\) are the trial function space and the test function space, respectively. Regarding the time interval \([t_{ntime interval \([t_{n},t_{n+1}]\) is divided into \(m\) (\(m\) = 1, 2, 3,... ) intervals so that the \(m\)-degree interpolation polynomial is used.  
\[\mathbf{u}^{h}=\mathbf{N}\mathbf{d}=\sum\limits_{i=1}^{m+1}N_{i}(t)\mathbf{d}_ {i},\ \mathbf{v}^{h}=\mathbf{N}\mathbf{d}_{v}=\sum\limits_{i=1}^{m+1}N_{i}(t) \mathbf{v}_{i} \tag{7}\]  
where \(\mathbf{N}=[N_{1}(t)\mathbf{I},N_{2}(t)\mathbf{I},\ldots,N_{m+1}(t)\mathbf{I}]\) is the shape function matrix and \(m\) is the degree of time elements. \(N_{i}(t)\) (\(i\) = 1, 2,..., \(m+1\)) is the shape function of Lagrange polynomial form and \(\mathbf{I}\) is the unit matrix. \(\mathbf{d}=[\mathbf{d}_{1}^{T},\mathbf{d}_{2}^{T},\ldots,\mathbf{d}_{m+1}^{T} ]^{T}\) is the matrix of nodal displacement vectors in time interval \([t_{n},t_{n+1}]\) and \(\mathbf{v}=[\mathbf{v}_{1}^{T},\mathbf{v}_{2}^{T},\ldots,\mathbf{v}_{m+1}^{T} ]^{T}\) is the matrix of nodal test displacement vectors.  
Substituting Equation (7) into Equation (6) and implementing partial integration, the following equation can be derived by coordinate mapping, in which time \(t\in[t_{n},t_{n+1}]\) is converted to \(\zeta\in[-1,1]\).  
\[(\frac{-4}{\Delta t^{2}}\mathbf{{}^{m}}\mathbf{{}^{K}}+\frac{2}{\Delta t} \mathbf{{}^{c}}\mathbf{{}^{K}}+\mathbf{{}^{k}}\mathbf{{}^{K}})\mathbf{d}+ \frac{2\mathbf{M}}{\Delta t}(\mathbf{{}^{N}}(1)^{T}\dot{\mathbf{u}}_{n+1}- \mathbf{{}^{N}}(-1)^{T}\dot{\mathbf{u}}_{n})=\mathbf{P} \tag{8}\]  
where  
\[\mathbf{{}^{m}}\mathbf{{}^{K}}=\int_{-1}^{1}\dot{\mathbf{{}^{N}}}\mathbf{{}^{ T}}\mathbf{{}^{M}}\mathbf{{}^{M}}\mathbf{{}^{\zeta}},\ \mathbf{{}^{c}}\mathbf{{}^{K}}=\int_{-1}^{1}\mathbf{{}^{T}}\mathbf{{}^{\overline {C}}}(\zeta)\dot{\mathbf{{}^{N}}}\mathbf{{}^{d}}\zeta \tag{9}\]  
\[\mathbf{{}^{k}}\mathbf{{}^{K}}=\int_{-1}^{1}\mathbf{{}^{N}}\mathbf{{}^{T}} \mathbf{{}^{\overline{K}}}(\zeta)\mathbf{{}^{N}}\mathbf{{}^{d}}\zeta,\mathbf{{ }^{P}}=\int_{-1}^{1}\mathbf{{}^{T}}\mathbf{{}^{\underline{R}}}\mathbf{{}^{ \underline{\zeta}}} \tag{10}\]  
From Equation (8), for any positive integer \(m\), we can derive and obtain a time integration scheme. In fact, considering both the calculation efficiency and accuracy, letting \(m\) = 1 or 2 or 3 is sufficient [31]. At the same time, it can be proved that when reduced Gaussian integration is used, the derived scheme is unconditionally stable for linear problems [33]. When \(m\) = 1 and using the Gaussian integration, the time integration scheme from Equation (8) can be presented as  
\[\left[\begin{array}{cc}\mathbf{{}^{K}_{12}^{(k)}}&0\\ \mathbf{{}^{K}_{22}^{(k)}}&\frac{2\mathbf{M}}{\Delta t}\end{array}\right] \left\{\begin{array}{c}\mathbf{{}^{u}_{n+1}}\\ \mathbf{{}^{\underline{\zeta}}_{n+1}}\end{array}\right\}=\left\{\begin{array}[ ]{cc}\mathbf{{}^{p}_{1}}^{(k)}-\mathbf{{}^{K}_{11}^{(k)}}\mathbf{{}^{