In the present work, a non-linear finite element model for the analysis of composite and sandwich plates is proposed. The underlying variational formulation is based on the kinematics of the refined zigzag theory (RZT) as well as on a modified version of the two-field Hellinger-Reissner (HR) functional, where the displacements and the transverse shear stresses are involved as independent field variables. In fact, a consistent expression for each shear stress function is obtained by integrating the local equilibrium equations written in terms of derivated strain components. Regardless of the number of material layers, the proposed four-node plate element, which is labeled HR-RZT, exhibits only seven nodal degrees of freedom. The additionally introduced variables are effectively eliminated on the element level. One key advantage of the HR-RZT formulation is the fact that interlayer-continuity of the transverse shear stresses is automatically obtained without resorting to any post-processing procedures. The same applies to the zero stress conditions at the outer surfaces of the laminate. Further, due to the enhanced kinematics, the layerwise distortion of the laminate's cross-section can also be accurately captured. The performance of the novel element formulation is studied by several numerical applications, where the solutions of the two-dimensional HR-RZT elements are verified by comparison with fully three-dimensional finite element models.