  In order to solve the magnetohydrodynamics (MHD) equations with a
$\mathbf{\mathcal{H}}(\mathbf{div})$-conforming element, a novel approach is
proposed to ensure the exact divergence-free condition on the magnetic field.
The idea is to add on each element an extra interior bubble function from
higher order hierarchical $\mathbf{\mathcal{H}}(\mathbf{div})$-conforming
basis. Four such hierarchical bases for the $\mathbf{\mathcal{H}}%
(\mathbf{div})$-conforming quadrilateral, triangular, hexahedral and
tetrahedral elements are either proposed (in the case of tetrahedral) or
reviewed. Numerical results have been presented to show the linear independence
of the basis functions for the two simplicial elements. Good matrix
conditioning has been confirmed numerically up to the fourth order for the
triangular element and up to the third order for the tetrahedral element.
