## 2 Basic Theory of HSDTPEs  
### Incompatible Polynomial Approximation over Each Triangular Element  
Consider a linear elastic plate with a length \(a\), width \(b\), and thickness \(h\) undergoing infinitesimal deformation, as illustrated in Figure 1. The mid-plane of the plate is divided into arbitrary triangular elements, as shown in Figure 2. The displacement function of the most well-known HSDTs [17] for each triangular element \(e_{i}\) is generally defined by:  
\[\left\{ \begin{array}{rl}u(x,y,z)&=u_{0}(x,y)-z\theta_{x}\ +g(z)\Big{(}\frac{ \partial w_{0}}{\partial x}-\theta_{x}\Big{)}\\ v(x,y,z)&=v_{0}(x,y)-z\theta_{y}\ +g(z)\Big{(}\frac{\partial w_{0}}{ \partial y}-\theta_{y}\Big{)}\\ w(x,y,z)&=w_{0}(x,y)\end{array} \right. \tag{1}\]  
where \(u_{0},\ v_{0}\) and \(w_{0}\) are the in-plane and transverse displacements at the mid-plane, respectively; \(\mathbf{u}=[u,v,w]^{\text{T}}\) denotes the displacements of a point \(\mathbf{x}\) on the plate; \(\theta_{x}\) and \(\theta_{y}\) are the rotations of the normal to the cross section; \(z\) is the coordinate in the transverse direction; and \(g(z)\) describes the distribution of shear effect in the thickness direction. A review of transverse shear functions \(g(z)\) can be found in [PERSON] et al. [56]. For isotropic plates with infinitesimal strains, the in-plane displacements \(u_{0}(x,y)\) and \(v_{0}(x,y)\) can be neglected because the thickness \(h\) is much smaller than the characteristic length \(a\) and \(b\), and the transverse displacement is much smaller than the thickness \(h\), which leads to \(u_{0}(x,y)\approx 0\) and \(v_{0}(x,y)\approx 0\) at the mid-plane. The transverse normal displacement \(w\) can also be assumed as \(w=w(x,y,z)\), which is not a constant along the \(z\) axis, and can be defined by the ESL or LW models [30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43] to capture the effective mechanical behavior along the thickness of composite structures well. Figure 1: A linear elastic plate. To demonstrate the performance of the present theory for various transverse shear functions, the third-order shear function \(g(z)=-\frac{4\textbf{e}^{3}}{3\hbar^{2}}\)[17] and the fifth-order shear function \(g(z)=-\frac{z}{8}-\frac{2\textbf{e}^{3}}{\hbar^{2}}+\frac{2\textbf{e}^{5}}{ \hbar^{4}}\)[19] are used to develop the Third-order Shear Deformation Triangular Plate Element (TfSDTPE) and Fifth-order Shear Deformation Triangular Plate Element (FfSDTPE), respectively. For the special case \(g(z)=0\), Equation (1) is actually the expression of Mindlin plate theory (or FSDT). The corresponding plate element is referred to as FiSDTPE (First-order Shear Deformation Triangular Plate Element) for comparison in the paper. We assume that:  
\[\left\{ \begin{array}{l}u_{0}=\textbf{P}^{2}\textbf{a}^{u_{0}}\\ v_{0}=\textbf{P}^{2}\textbf{a}^{v_{0}}\\ \theta_{x}=\textbf{P}^{2}\textbf{a}^{\theta_{x}}\\ \theta_{y}=\textbf{P}^{2}\textbf{a}^{\theta_{y}}\\ w_{0}=\textbf{P}^{3}\textbf{a}^{w}\end{array} \right. \tag{2}\]  
where \(\textbf{a}^{u_{0}}=\left[ \begin{array}{cccc}a_{1}&a_{2}&\cdots&a_{6}\\ \end{array} \right]^{\text{T}}\), \(\textbf{a}^{v_{0}}=\left[ \begin{array}{cccc}a_{7}&a_{8}&\cdots&a_{12}\\ \end{array} \right]^{\text{T}}\), \(\textbf{a}^{\theta_{x}}=\left[ \begin{array}{cccc}a_{13}&a_{14}&\cdots&a_{18 }\\ \end{array} \right]^{\text{T}}\), \(\textbf{a}^{\theta_{y}}=\left[ \begin{array}{cccc}a_{19}&a_{20}&\cdots&a_{24 }\\ \end{array} \right]^{\text{T}}\), \(\textbf{a}^{w}=\left[ \begin{array}{cccc}a_{25}&a_{26}&\cdots&a_{34}\\ \end{array} \right]^{\text{T}}\) are the vector of generalized approximation DOFs (degrees of freedom) of the triangular element \(e_{i}\), \(\textbf{P}^{2}\) is the second-order polynomial basis function, and \(\textbf{P}^{3}\) is the third-order polynomial basis function in which:  
\[\textbf{P}^{2}(\textbf{x})=\left[ \begin{array}{cccccc}1&x_{0}&y_{0}&x_{0}^ {2}&x_{0}y_{0}&y_{0}^{2}\\ \end{array} \right] \tag{3}\]  
\[\textbf{P}^{3}(\textbf{x})=\left[