  Motivated by the observation of comets and asteroids rotating in
non-principal axis (NPA) states, we investigate the relaxation of a freely
precessing triaxial ellipsoidal rotator towards its lowest-energy spin state.
Relaxation of the precession arises from internal dissipative stresses
generated by self-gravitation and inertial forces from spin. We develop a
general theory to determine the viscoelastic stresses in the rotator, under any
linear rheology, for both long-axis (LAM) and short-axis (SAM) modes. By the
methods of continuum mechanics, we calculate the power dissipated by the stress
field and the viscoelastic material strain which enables us to determine the
timescale of the precession dampening. To illustrate how the theory is used, we
apply our framework to a triaxial 1I/2017 (`Oumuamua) and 4179 Toutatis under
the Maxwell regime. For the former, employing viscoelastic parameters typical
of very cold monolithic asteroids renders a dampening timescale longer by a
factor of $10^{10}$ and higher than the timescales found in the works relying
on the $\,Q$-factor approach, whilst the latter yields a significantly shorter
timescale as a consequence of including self-gravitation. We further reduce our
triaxial theory to bodies of an oblate geometry and derive a family of
relatively simple analytic approximations determining the NPA dampening times
for Maxwell rotators, as well as a criterion determining whether
self-gravitation is negligible in the relaxation process. Our approximations
exhibit a relative error no larger than $0.2\%$, when compared to numerical
integration, for close to non-dissipative bodies and $0.002\%$ for highly
energy dissipating rotators.
