(eng): The quasistatic, Prandtl-Reuss perfect plasticity at small strains is combined with a gradient, reversible (i.e., admitting healing) damage which influences both the elastic moduli and the yield stress. Existence of weak solutions of the resulting system of variational inequalities is proved by a suitable fractional-step discretization in time with guaranteed numerical stability and convergence. After finite-element approximation, this scheme is computationally implemented and illustrative two-dimensional simulations are performed. The model allows, e.g., for application in geophysical modeling of reoccurring rupture of lithospheric faults. Resulting incremental problems are solved in MATLAB by quasi-Newton method to resolve the elastoplasticity component of the solution, while the damage component is obtained by solving a quadratic programming problem.