We discuss the numerical simulation of visco-elastodynamics based on finite element schemes for an inherent (singular hyperbolic) relaxation approximation. Those methods use piecewise polynomials of arbitrary degree, so their consistency error is of high order. Stability properties of the relaxation finite element schemes are derived from estimating suitable modifications of the standard energy functional, the so-called modulated energy, that is suggested by the presence of relaxation terms. Similar techniques are also applied to prove error estimates and deduce the convergence to the (smooth) solutions of the equations of visco-elastodynamics. The computational performance of these schemes improves when combined with an adaptive meshing strategy, that yields an extra stabilization mechanism against the effects of the nonlinear response over shock regions of the solution.