SUMMARYWe present a finite element method for the analysis of ductile crystals whose energy depends on the density of geometrically necessary dislocations (GNDs). We specifically focus on models in which the energy of the GNDs is assumed to be proportional to the total variation of the slip strains. In particular, the GND energy is homogeneous of degree one in the slip strains. Such models indeed arise from rigorous multiscale analysis as the macroscopic limit of discrete dislocation models or from phenomenological considerations such as a line‐tension approximation for the dislocation self‐energy. The incorporation of internal variable gradients into the free energy of the system renders the constitutive model non‐local. We show that an equivalent free‐energy functional, which does not depend on internal variable gradients, can be obtained by exploiting the variational definition of the total variation. The reformulation of the free energy comes at the expense of auxiliary variational problems, which can be efficiently solved using finite element approximations. The addition of surface terms in the formulation of the free energy results in additional boundary conditions for the internal variables. The proposed framework is verified by way of numerical convergence tests, and simulations of three‐dimensional problems are presented to showcase its applicability. A performance analysis shows that the proposed framework solves strain‐gradient plasticity problems in computing times of the order of local plasticity simulations, making it a promising tool for non‐local crystal plasticity three‐dimensional large‐scale simulations. Copyright © 2012 John Wiley & Sons, Ltd.