The 2-D approximation functions based on a general exact 3-D plate solution are used to derive locking-free, rectangular, 4-node Mindlin (i.e., first-order plate theory), Levinson (i.e., a third-order plate theory), and Full Interior plate finite elements. The general plate solution is defined by a biharmonic mid-surface function, which is chosen for the thick plate elements to be the same polynomial as used in the formulation of the well-known nonconforming thin Kirchhoff plate element. The displacement approximation that stems from the biharmonic polynomial satisfies the static equilibrium equations of the 2-D plate theories at hand, the 3-D Navier equations of elasticity, and the Kirchhoff constraints. Weak form Galerkin method is used for the development of the finite element model, and the matrices for linear bending, buckling and dynamic analyses are obtained through analytical integration. In linear buckling problems, the 2-D Full Interior and Levinson plates perform particularly well when compared to 3-D elasticity solutions. Natural frequencies obtained suggest that the optimal value of the shear correction factor of the Mindlin plate theory depends primarily on the boundary conditions imposed on the transverse deflection of the 3-D plate used to calibrate the shear correction factor.