Most of mechanical systems and complex structures involve shell components (for instance, in additive manufacturing, in body-in-white manufacturing or in sheet metal forming). Therefore, the 2D simulation based on shell theory, appears as an appealing choice in structural analysis since it reduces the computational complexity. Nevertheless, this 2D framework may fail to capture rich physics compromising the usual assumptions made in shell theories. In these scenarios, a 3D discretization seems compulsory to ensure accurate results. Unfortunately, the resulting meshes usually involve too many degrees of freedom and, as a result, solving a model requires a significant amount of CPU time. This work deals with a novel hybrid finite element formulation for thin to thick structures based on enriching the first-order Reissner-Mindlin theory. Such formulation preserves the in-plane kinematics of the classical first-order theory, while enriching the out-of-plane one via higher-order polynomials. This is done without increasing the number of degrees of freedom of the standard shell element. The out-of-plane polynomial coefficients are predicted using a machine learning (ML) model. Indeed, the new finite element does not neglect the effect of the transverse stress through the normal direction. As usual, the machine learning framework consists of two phases. A preliminary offline phase where, for a specific configuration, the out-of-plane rich kinematics is learned from fully 3D computations. In the online phase, such rich out-of-plane kinematics is integrated into a 2D shell computation. The proposed methodology enhances the quality of classic shell-based simulations when the hypotheses of the 2D theory are compromised, while keeping the same computational cost. First results focus on elastic plates where classical shell elements give inaccurate mechanical predictions. The study opens new perspectives in computational inelasticity and fracture of thin and thick structures.