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 resulted system of variational inequalities is proved by a suitable fractional-step discretisation in time with guaranteed numericalstability and convergence. After finite-element approximation, this scheme is computationally implemented and illustrative 2-dimensional simulations are performed. The model allows e.g. for application in geophysical modelling of re-occurring rupture of lithospheric faults. Resulted incremental problems are solved in MATLAB by quasi-Newton method to resolve elastoplasticity component of the solution while damage component is obtained by solution of a quadratic programming problem.1 2. The model, its weak formulation, and existence result. Hereafter, we suppose that the damageable elasto-plastic body occupies a bounded smooth domain Ω ⊂ R d , d = 2 or 3. We denote by n the outward unit normal to ∂Ω. We further suppose that the boundary of Ω splits aswith Γ D and Γ N open subsets in the relative topology of ∂Ω, disjoint one from each other, each of them with a smooth ((d−1)-dimensional) boundary, and covering ∂Ω up to (d−1)-dimensional zero measure. Considering T > 0 a fixed time horizon, we set Q := (0, T )×Ω, Σ := (0, T )×Γ, Σ D := (0, T )×Γ D , Σ N := (0, T )×Γ N .