In this work, we revisit the stress resultant elastoplastic geometrically exact shell finite element formulation that is based on the Ilyushin–Shapiro two-surface yield function with isotropic and kinematic hardening. The main focus is on implicit projection algorithms for computation of updated values of internal variables for stress resultant shell elastoplasticity. Four different algorithms are derived and compared. Three of them yield practically identical final results, yet they differ considerably in computational efficiency and implementation complexity, since they solve different sets of equations and they use different procedures that choose active yield surfaces. One algorithm does not provide acceptable accuracy. It turns out that the most simple and straightforward algorithm performs surprisingly well and efficiently. Several numerical examples are presented to illustrate the Ilyushin–Shapiro stress resultant shell formulation and the numerical performance of the presented integration algorithms.