The paper outlines a numerical procedure for solving physically and geometrically nonlinear problems of statics for thin shells based on three mesh-based methods: finite-difference, variational difference, and finite-element methods. The methodological, algorithmic, and analytical aspects of implementing the Kirchhoff-Love hypotheses are analyzed. The algorithmic approach employs Lagrangian multipliers. The advantages and disadvantages of these methods are evaluated Keywords: finite-difference method, variational difference method, finite-element method, Kirchhoff-Love hypotheses, Lagrangian multipliers, nonlinear elastic composite, elastoplastic stateIntroduction. Strength analysis of structural members such as shells of complex geometry is mainly conducted on a computer using numerical mesh-based methods, including the finite-difference method (FDM), the variational difference method (VDM), and the finite-element method (FEM). To derive the basic equations of the theory of thin shells, it is sufficient to use the classical model based on the Kirchhoff-Love hypotheses. The Kirchhoff-Love kinematics is implemented differently in these three methods.The system of governing equations in the FDM is known to be derived from the equilibrium equations for an element of the shell, while the Kirchhoff-Love constraints are incorporated analytically, by substituting the corresponding relations into the formulas describing the distribution of displacements throughout the thickness of the shell. This way causes no complications. Contrastingly, the system of governing equations in the VDM and FEM is derived from the extremum condition for some functional. An attempt to implement Kirchhoff-Love constraints analytically gives rise to higher-than-first-order derivatives of displacements in the functionals. Note that the complications in the VDM are only restricted to more awkward equations and the necessity of deriving them because the VDM imposes no special requirements upon the smoothness of the integrands and functions. In the FEM, however, certain requirements are placed upon the smoothness of the coordinate functions, and the construction of an element is complicated if the integrands include higher-order derivatives.These circumstances led to a search for alternative ways of incorporating the Kirchhoff-Love constraints such as those that would not produce derivatives of higher than the first order in the integrands. Noteworthy are two groups of approaches that employ Lagrangian multipliers and penal functions. Applying them requires certain correctness [17,52] associated with the variation order and additional conditions. However, the use of Lagrangian multipliers to implement the Kirchhoff-Love hypotheses increases the number of varied functions and requires much computing. This is obviously why the Kirchhoff-Love constraints were first (1967) used [46] in the FEM for a plate by specifying them at discrete points of an element (discrete Kirchhoff element). Later, this method was extended to thin shells [45,70] and is wi...