  In this paper, we develop a class of mixed finite element methods for the
ferrofluid flow model proposed by Shliomis [Soviet Physics JETP, 1972]. We show
that the energy stability of the weak solutions to the model is preserved
exactly for both the semi- and fully discrete finite element solutions.
Furthermore, we prove the existence and uniqueness of the discrete solutions
and derive optimal error estimates for both the the semi- and fully discrete
schemes. Numerical experiments confirm the theoretical results.
