An alternative approach is proposed and applied to solve boundary value problems within the strain gradient elasticity theory. A mixed variation formulation of the finite element method (FEM) based on the concept of the Galerkin method is used. To construct finite-dimensional subspaces separate approximations of displacements, deformations, stresses, and their gradients are implemented by choosing the different sets of piecewise polynomial basis functions, interrelated by the stability condition of the mixed FEM approximation. This significantly simplifies the pre-requirement for approximating functions to belong to class C1 and allows one to use the simplest triangular finite elements with a linear approximation of displacements under uniform or near-uniform triangulation conditions. Global unknowns in a discrete problem are nodal displacements, while the strains and stresses and their gradients are treated as local unknowns. The conditions of existence, uniqueness, and continuous dependence of the solution on the problem’s initial data are formulated for discrete equations of mixed FEM. These are solved by a modified iteration procedure, where the global stiffness matrix for classical elasticity problems is treated as a preconditioning matrix with fictitious elastic moduli. This avoids the need to form a global stiffness matrix for the problem of strain gradient elasticity since it is enough to calculate only the residual vector in the current approximation. A set of modeling plane crack problems is solved. The obtained solutions agree with the results available in the relevant literature. Good convergence is achieved by refining the mesh for all scale parameters. All three problems under study exhibit specific qualitative features characterizing strain gradient solutions namely crack stiffness increase with length scale parameter and cusp-like closure effect. [ABSTRACT FROM AUTHOR]