## 1 Introduction  
The classic Kirchhoff-Love plate theory based on the assumption that a plane section perpendicular to the mid-plane of the plate before deformation remains plane and perpendicular to the deformed mid-plane after deformation is the simplest plate theory in engineering analysis. However, the Kirchhoff-Love plate theory is only applicable for thin plates due to the neglecting of the shear deformation effects. The most well-known and earliest plate theories that take into account the shear deformation effects were proposed by [PERSON] [1] and [PERSON] [2], in which the Mindlin plate theory was based on an assumption of a linear variation of in-plane displacements through the thickness of the plate, referred to as the First-order Shear Deformation Theory (FSDT). The plate elements derived from the FSDT only require C\({}^{0}\) continuity in approximation fields, have the advantages of physical clarity and simplicity of application [3], and hence were widely accepted and used to model thick-thin plates by scientists and engineers. Unfortunately, the FSDT elements suffer from the shear-locking problem when the thickness to length ration of the plate becomes very small, due to inadequate dependence among transverse deflection and rotations using an ordinary low-order finite element [4]. Quite a large number of techniques have been developed to overcome this problem, such as the assumed shear strain approach, the discrete Kirchhoff/Mindlin representation, the mixed/hybrid formulation, and the reduced/selected integration [5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15]. These formulations are free from shear locking and are applicable to a wide range of practical engineering problems, but in general, it is rather complex and time consuming to include the transverse shear effects for thick plates, which would also lead to complexity and difficulty in the programming. Moreover, the assumption of FSDT causes constant transverse shear strains and stresses across the thickness, which violates the conditions of zero transverse shear stresses on the top and bottom surfaces of plates. A shear correction factor is therefore required to properly compute the transverse shear stiffness. The finding of such a shear correction factor in FSDT is difficult since it depends on geometric parameters, material, loading and boundary conditions, etc. [16]. In recent years, High-order Shear Deformation Theories (HSDTs) have gained extensive recognition and made quite great progress [4; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43]. Based on polynomial or non-polynomial transverse shear functions, various HSDTs have been proposed to avoid the use of a shear correction factor, and to better predict the shear behavior of the plate, for instance, the third-order shear deformation theory [17; 18], the fifth-order shear deformation theory [19], the exponential shear deformation theory [20], the hyperbolic shear deformation theory [21], and the combined or mixed HSDTs [22; 23]. Please see [PERSON] and [PERSON] [16] and [PERSON] et al. [24] for a comprehensive review of HSDTs. In HSDTs, the bending angles of rotation and shear angles can be treated as independent variables, and the shear-locking problem encountered in FSDT can be well-solved [4]. In [30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43], two well-know HSDTs named as equivalent single layer (ESL) and layer-wise (LW) models are developed to evaluate the effective mechanical behavior of composite structures correctly. The accuracy and reliability of HSDTs have been illustrated by numerous examples in the literature [4; 17; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42; 43]. However, the general requirement of C\({}^{1}\) continuity in approximation fields in HSDTs brings difficulties for the numerical implementation of the standard Finite Element Method (FEM), which is similar to that of the classic Kirchhoff-Love plate theory.