## 4 Generalized Differential Quadrature Method  
In the present paper, the governing equations shown in matrix form in (38) are solved by means of the GDQ method. According to this numerical technique, it is possible to approximate the \(n\)-th derivative of a sufficiently smooth function at a generic point within the reference domain as a weighted linear sum of the values that the function assumes in some chosen grid points. Starting from this definition, it is clear that the strong form of the governing equations is solved. Let us consider a two-dimensional domain, as the one needed to describe any doubly-curved surface, in which a generic function \(f(\alpha_{1},\alpha_{2})\) is defined. Three types of derivative can be computed at a generic point \(\big{(}\alpha_{1i},\alpha_{2j}\big{)}\) for \(i=1,2,\ldots,I_{N}\) and \(j=1,2,\ldots,I_{M}\): the \(n\)-th derivative with respect to the first coordinate \(\alpha_{1}\)  
\[\frac{\partial^{n}f(\alpha_{1},\alpha_{2})}{\partial\alpha_{1}^{n}}\bigg{|}_ {\alpha_{1}=\alpha_{1i},\alpha_{2}=\alpha_{2j}}\cong\sum_{k=1}^{I_{N}}\xi_{ \alpha_{1}(ik)}^{(n)}f\big{(}\alpha_{1k},\alpha_{2j}\big{)} \tag{42}\]  
the \(m\)-th derivative with respect to the second coordinate \(\alpha_{2}\)  
\[\frac{\partial^{m}f(\alpha_{1},\alpha_{2})}{\partial\alpha_{2}^{m}}\bigg{|}_ {\alpha_{1}=\alpha_{1i},\alpha_{2}=\alpha_{2j}}\cong\sum_{l=1}^{I_{M}}\xi_{ \alpha_{2}(jl)}^{(m)}f(\alpha_{1i},\alpha_{2l}) \tag{43}\]  
and the mixed derivative of order \((n+m)\)  
\[\frac{\partial^{n+m}f(\alpha_{1},\alpha_{2})}{\partial\alpha_{1}^{n}\partial \alpha_{2}^{m}}\bigg{|}_{\alpha_{1}=\alpha_{1i},\alpha_{2}=\alpha_{2j}}\cong \sum_{k=1}^{I_{N}}\xi_{\alpha_{1}(ik)}^{(n)}\left(\sum_{l=1}^{I_{M}}\xi_{ \alpha_{2}(jl)}^{(m)}f(\alpha_{1k},\alpha_{2l})\right) \tag{44}\]It is pointed out that \(I_{N}\), \(I_{M}\) denote the overall number of discrete points within the reference domain. The symbols \(\zeta_{a_{1}(ik)}^{(n)}\) and \(\zeta_{a_{2}(jl)}^{(m)}\) are introduce to specify the weighting coefficients for the derivatives at issue. It is clear that the key points of the present numerical approach are the choice of a proper discrete grid point distribution and the evaluation of the weighting coefficients. If the Lagrange polynomials are chosen to approximate the function, there is no restriction on the choice of the grid points. This aspect represents one of the main advantages of the GDQ method. In addition, the weighting coefficients can be easily evaluated if the recursive formulae given by [PERSON] are employed [90]. Further details concerning this technique and its major applications can be found in the review paper by [PERSON] et al. [91]. As far as the choice of the grid points is concerned, several distributions can be used, as highlighted in the work [91]. Some examples of grid distributions typically employed in many structural problems are listed in Table 1. It should be recalled that a discrete grid must be applied along the two principal coordinates \(\alpha_{1},\alpha_{2}\) as shown in the following relations  
\[\alpha_{1i}=\frac{\alpha_{1}^{1}-\alpha_{1}^{0}}{r_{I_{N}}-r_{1}} (r_{i}-r_{1})+\alpha_{1}^{0} \tag{45}\] \[\alpha_{2j}=\frac{\alpha_{2}^{1}-\alpha_{2}^{0}}{r_{I_{M}}-r_{1}} (r_{j}-r_{1})+\alpha_{2}^{0} \tag{46}\]  
The values of \(r_{i}\) and \(r_{j}\), for \(i=1,2,\ldots,I_{N}\) and \(j=1,2,\ldots,I_{M}\), can be found in Table 1.