In this work, we describe a finite element (FE) implementation of an elasto-plastic self-consistent (EPSC) polycrystal plasticity model termed FE-EPSC, which is intended for simulations of metal forming. To this end, we present an analytical Jacobian, which is necessary for the implicit coupling and ensuring a fast convergence. Every FE integration point is a material point that can be represented either by a single crystal or a polycrystalline material. The constituent crystal can deform by a combination of anisotropic elasticity, crystallographic slip, and deformation twinning. The model is validated and applied to a suite of tests, including monotonic compression, cyclic forward loading, unloading and reverse loading and non-monotonic four-point bending, and materials, such as different alloy compositions, crystal structures, and initial microstructures. The same FE-EPSC framework is applied for all these cases with the main differences pertaining to intrinsic properties, such as the available slip and twinning deformation modes, and the material parameters for activating and hardening of these modes. Full characterization for these parameters for high-purity α-Ti is presented here for the first time. Through these examples we show that, in addition to being predictive with great accuracy, the key advantage of this model lies in its versatility. It accounts for the development of backstress aided dislocation glide, thermally activated storage of dislocations, elastic anisotropy, crystallographic slip and deformation twinning via multiple modes, and de-twinning as well as multi-level homogenizations.