## 3 Implementation  
### Discretization Using Isoparametric and High-Order Spectral Elements in the SBFEM  
The displacement expressions, discretized over the element domain, can be written in terms of the shape functions \(Ni(\eta)\) as shown in Equation (5), in which the displacements of each point along the \(z\)-axis is a function of the vertical position of the plate thickness, and the nodal unknown displacements, \((U_{yi},U_{zi})\), are expressed in cartesian coordinates,  
\[u(z,\eta)=N(\eta)u_{n}(z). \tag{5}\]  
Because the shape functions depend only on \(\eta\), the derivatives can be separated in terms of the two spatial coordinates for the strain (Equation (6)) and for the stress-displacement relationship (Equation (7),  
\[\varepsilon(z,\eta)=B_{1}u_{,z}(z)+B_{2}u_{n}(z), \tag{6}\]  
\[\sigma=D(B_{1}u_{n,z}(z)+B_{2}u_{n}(z), \tag{7}\]  
where \(B_{1}\) and \(B_{2}\) include the shape functions:  
\[B_{1}=b_{1}N\,\ B_{2}=\frac{1}{y,n}\mathrm{b}_{2}N_{\eta}, \tag{8}\]  
\[N(\eta)= \begin{bmatrix}N_{1}(\eta)&0&N_{2}(\eta)&0&N_{3}(\eta)&0\\ 0&N_{1}(\eta)&0&N_{2}(\eta)&0&N_{3}(\eta)\end{bmatrix} , \tag{9}\]  
with \(N_{1}\), \(N_{2}\) and \(N_{3}\) expressed in the local coordinates \(\eta\). Isoparametric elements with three nodes equally spaced are located at local coordinates \(\eta 0=-1\), \(\eta 1=0\) and \(\eta 2=+1\), so that an element of order \(p\) contains \(p+1\) nodes. The curves in Figure 2 illustrate the Lagrange shape functions of a one-dimensional isoparametric element corresponding to the polynomial of order \(p=2\) and the quadratic interpolation functions used as indicated in Equation (10),  
\[N_{1}=\frac{1}{2}(\eta^{2}-\eta),N_{2}=(1-\eta^{2}),N_{3}=\frac{1}{2}(\eta^{2 }+\eta). \tag{10}\]  
Although Lagrange interpolation polynomials may, in principle, be defined for any set of distinct nodes, their properties are critically dependent on node positions. On the other hand, other different nodal distributions can be proposed, so this study employs an example of the Gauss-Lobatto-Legendre nodal distribution, as shown in Figure 2. In this case, the two nodes are located at the extremes \(\eta 1=-1\) and \(\eta_{p\ +1}=1\)[31]. The \(p-1\) nodes, which correspond to internal points according to Gauss-Lobatto-Legendre, GLL, are obtained when they meet the condition of Equation (11),  
\[\frac{d}{d\eta}P_{p}(\eta_{i})=0,i=2\ldots p, \tag{11}\]  
where \(P_{p}\) denotes the Legendre polynomial of order \(p\). The Gauss-Lobatto-Legendre quadrature points are calculated for an element of the given order as follows: For a one-dimensional element of order \(p\) with local coordinates \(\eta\), ranging from \(-1\) to \(1\), the internal GLL nodes are the roots of the Lobatto polynomial (\(L\)) of order \(p-1\), which is defined as the first derivative of the Legendre polynomial (\(P\)) of order \(p\), so that:  
\[L_{o_{p-1}}(\eta)=\frac{d}{d\eta}P_{p}(\eta_{i}) \tag{12}\]The derivative of the Legendre polynomial \(P_{p}\) of order \(p-1\) yields the Lobatto polynomial, and its roots correspond to the location \(\eta\) of the nodes--Equation (12). These Legendre polynomials for the given order can be obtained as follows, in Equations (13)-(15),  
\[P_{0}(\eta)=1, \tag{13}\]  
\[P_{1}(\eta)=\eta, \tag{14}\]  
\[P_{k}(\eta)=\frac{2k-1}{k}\eta L_{k-1}(\eta)-\frac{k-1}{k}L_{k-2}(\eta),k=2,3,4\ldots \tag{15}\]  
Once the location of the nodes is known, the corresponding shapes functions \(N_{i}(\eta)\) of p order can be constructed using the Lagrange interpolation polynomials [40], obtained from Equation (16),  
\[N_{i}(\eta)=\prod_{ \begin{array}{c}j=1\\ j\
eq i\end{array} }^{p+1}\frac{\eta-\eta_{j}}{\eta_{i}-\eta_{j}},i=1,2,\ldots,p +1. \tag{16}\]  
According to the proposed scheme, the first five Legendre polynomials and the corresponding Lobatto polynomials have been obtained as examples.