Summary: The lay-out of a numerical solution procedure for the strain gradient flow (rate-independent) theory by N. A. Fleck and J. R. Willis [J. Mech. Phys. Solids 57, No. 7, 1045–1057 (2009; Zbl 1173.74316)] has been an open issue, and its finite element implementation is yet to be completed. Only recently, a sound solution procedure that allows for elastic-plastic loading/unloading and the interaction of multiple plastic zones has been put forward, and demonstrated within a 1D model set-up. The aim of the present work is to extend this procedure to form a basis for 2D and 3D calculations, and thereby to broaden its use within engineering applications. Focus is on the numerical implementation that adopts image analysis techniques to identify individual plastic zones and to treat none-regular mesh configurations. The developed finite element model is demonstrated through; i) tensile loading of a finite homogeneous material slab that offers a simple interpretation of results, ii) tensile loading of a double edge notch tension specimen involving symmetry considerations, and iii) shearing of a periodically voided structure that displays the ability of the procedure to treat periodic boundary conditions. A comparison between the implemented flow theory model and the corresponding rate-dependent visco-plastic version of the model shows coinciding results in the limit of zero rate-sensitivity.