  We introduce energy-preserving integrators for nonholonomic mechanical
systems. We will see that the nonholonomic dynamics is completely determined by
a triple $({\mathcal D}^*, \Pi, \mathcal{H})$, where ${\mathcal D}^*$ is the
dual of the vector bundle determined by the nonholonomic constraints, $\Pi$ is
an almost-Poisson bracket (the nonholonomic bracket) and $\mathcal{H}:{\mathcal
D}^{*}\rightarrow \mathbb{R}$ is a Hamiltonian function. For this triple, we
can apply energy-preserving integrators, in particular, we show that discrete
gradients can be used in the numerical integration of nonholonomic dynamics. By
construction, we achieve preservation of the constraints and of the energy of
the nonholonomic system. Moreover, to facilitate their applicability to complex
systems which cannot be easily transformed into the aforementioned
almost-Poisson form, we rewrite our integrators using just the initial
information of the nonholonomic system. The derived procedures are tested on
several examples: A chaotic quartic nonholonomic mechanical system, the
Chaplygin sleigh system, the Suslov problem and a continuous gearbox driven by
an asymmetric pendulum. Their performace is compared with other standard
methods in nonholonomic dynamics, and their merits verified in practice.
