In this work, we describe in detail a model for small strain elasto-visco-plasticity with convex, non-smooth, yield functions and associative nonlinear kinetic laws, restricted to linear hardening. Using concepts of non-smooth convex geometry, numerical methods are developed to integrate the evolution equations of the model. These algorithms are analyzed and shown to inherit a discrete dissipation inequality, irrespective of the smoothness of the yield function. The results are applied to models based on Tresca's and Drucker-Prager's yield criteria, which include all possible types of non-smoothness. Numerical examples are shown to illustrate the performance of the methods. (C) 2017 Elsevier B.V. All rights reserved.