We present a novel multistage hybrid asymptotic–direct approach to the modeling of the nonlinear behavior of thin shells with piezoelectric patches or layers, which is formulated in a holistic form for the first time in this paper. The key points of the approach are as follows: (1) the asymptotic reduction in the three-dimensional linear theory of piezoelasticity for a thin plate; (2) a direct approach to geometrically nonlinear piezoelectric shells as material surfaces, which is justified and completed by demanding the mathematical equivalence of its linearized form with the asymptotic solution for a plate; and (3) the numerical analysis by means of an FE scheme based on the developed model of the reduced electromechanically coupled continuum. Our approach is illustrated by examples of static and steady-state analysis and verified with three-dimensional solutions computed with the commercially available FE code ABAQUS as well as by comparison with results reported in the literature.