To obtain a finite element algorithm in the three-field formulation, a conditional functional was used, based on the equality of the actual work of the internal force factors (internal forces and moments) at the strains and curvatures of the middle surface. As an addition, the functional assumed the condition that the work of the residual of internal forces at the strains and curvatures of the middle surface had been equal to zero. The difference between the adopted internal forces and the internal forces represented through the strains and curvatures of the middle surface according to Hooke's law was used as the residual of internal forces. A quadrilateral fragment of the middle surface of a thin shell was used as the finite element. Kinematic quantities (displacements and their first order derivatives), strain quantities (strains and curvatures of the middle surface) and force quantities (internal forces and moments) were taken as the nodal unknowns. Approximating expressions with Hermite polynomials of third degree were used to approximate kinematic quantities. The sought strain and force quantities were approximated through the corresponding nodal unknowns by bilinear shape functions. A finite element stiffness matrix (with the dimensions 36x36) with respect to the kinematic nodal unknowns was formed through the functional minimization. Specific examples showed the efficiency of the three-field finite element algorithm in determining the displacements, strains and internal forces.