A hyperbolic type integro-differential equation, modeling fractional order viscoelasticity, is considered. Well-posedness of the model problem is studied within two frameworks: the semigroup theory of linear operators and the Galerkin approximation approach. Then spatial semidiscretization of the problem is formulated by means of finite element method. Stability estimates are proved and optimal order a priori error estimates are obtained. Full discretization of the problem is also studied using continuous and discontinuous Galerkin methods. Stability and a priori error analysis are also studied
