## 5 Generalized Integral Quadrature Method  
A numerical approach has to be introduced to evaluate the integrals needed to compute the elastic constants \(A_{nm\;(pq)}^{(\tau\eta)}\), \(A_{nm\;(pq)}^{(\tau\eta)}\), \(A_{nm\;(pq)}^{(\tau\overline{\eta})}\), \(A_{nm\;(pq)}^{(\tau\overline{\eta})}\) defined in (30). In the present paper, the Generalized Integral Quadrature (GIQ) method is employed for this purpose. In fact, a numerical solution to the integrals in (30) must be found since it is impossible to solve them analytically due to their definition. According to this approach, the integral of a smooth function \(f(\zeta)\) defined in a close domain \([\zeta_{l},\zeta_{e}]\) \begin{table}
\begin{tabular}{c c} \hline
**Chebyshev I Kind (Cheb I)** & **Chebyshev II kind (Cheb II)** \\ \hline \(r_{k}=\cos\left(\frac{2(N-k)+1}{2N}\pi\right)\) & \(r_{k}=\cos\left(\frac{N-k+1}{N+1}\pi\right)\) \\ for \(k=1,2,\ldots,N,\quad r\in[-1,1]\) & for \(k=1,2,\ldots,N,\quad r\in[-1,1]\) \\ Chebyshev III kind (Cheb III) & Chebyshev IV kind (Cheb IV) \\ \(r_{k}=\cos\left(\frac{2(N-k)+1}{2N+1}\pi\right)\) & \(r_{k}=\cos\left(\frac{2(N-k+1)}{2N+1}\pi\right)\) \\ for \(k=1,2,\ldots,N,\quad r\in[-1,1]\) & for \(k=1,2,\ldots,N,\quad r\in[-1,1]\) \\ Chebyshev-Gauss (Cheb-Gau) & Chebyshev-Gauss-Lobatto (Cheb-Gau-Lob) \\ \(r_{1}=-1\), \(r_{N}=1\), \(r_{k}=\cos\left(\frac{2(N-k)-1}{2(N-2)}\pi\right)\) & \(r_{k}=\cos\left(\frac{N-k}{N-1}\pi\right)\) \\ for \(k=2,3,\ldots,N-1,\quad r\in[-1,1]\) & for \(k=1,2,\ldots,N,\quad r\in[-1,1]\) \\ Legendre-Gauss (Leg-Gau) & Legendre-Gauss-Lobatto (Leg-Gau-Lob) \\ \(r_{k}=\text{roots of }\left(1-r^{2}\right)\cdot L_{N-1}(r)\) & \(r_{k}=\text{roots of }\left(1-r^{2}\right)\cdot\frac{d}{dr}(L_{N}(r))\) \\ for \(k=1,2,\ldots,N,\quad r\in[-1,1]\) & for \(k=1,2,\ldots,N,\quad r\in[-1,1]\) \\ \hline \end{tabular}
\end{table} 
Table 1: Eight different kinds of grid distributions. The indices \(k\) and \(N\) specify respectively \(i\) and \(I_{N}\), when \(\alpha_{1}\) direction is taken into account. On the other hand, they correspond to \(j\) and \(I_{M}\) when \(\alpha_{2}\) direction is considered. Finally, if the grid distribution has to be applied along the coordinate \(\zeta\), one gets \(k=m\) and \(N=I_{T}\). It should be noted that \(L_{N}\) denotes the Legendre polynomials. can be evaluated as a weighted linear sum of the values that the function \(f(\zeta)\) assumes in each point of the reference domain. In other words, this definition can be written as follows  
\[\int\limits_{\zeta_{l}}^{\zeta_{e}}f(\zeta)d\zeta=\sum\limits_{k=1}^{I_{T}}{(w_{ek }-w_{lk})f(\zeta_{k})} \tag{47}\]  
where \(I_{T}\) represents the number of discrete points defined within the reference domain, whereas \(w_{ek},w_{lk}\) are the weighting coefficients for the integration. These last quantities can be computed directly from the definitions of the corresponding weighting coefficients of the first order derivatives, as illustrated in depth in the review paper by [PERSON] et al. [91]. It is clear that the GIQ method shares the same fundamentals with the GDQ method. Thus, the same recursive formulae given by [PERSON] can be used once again to evaluate the weighting coefficients at issue, once an arbitrary grid distribution is chosen to discretize the reference domain. It should be noticed that differently from the well-known Gaussian integration schemes, this method is able to compute the weighting coefficients without any restriction on the position of the points within the reference domain. Therefore, all the different kinds of grid distributions listed in Table 1 can be used with this aim. For the sake of clarity, it should be specified that the one-dimensional domain in which the integrals in (30) must be evaluated is defined along the coordinate \(\zeta\) of the local reference system. In other words, the integration procedure is performed through the shell thickness.