  In this paper, we will use the interior functions of an hierarchical basis
for high order $BDM_p$ elements to enforce the divergence-free condition of a
magnetic field $B$ approximated by the H(div) $BDM_p$ basis. The resulting
constrained finite element method can be used to solve magnetic induction
equation in MHD equations. The proposed procedure is based on the fact that the
scalar $(p-1)$-th order polynomial space on each element can be decomposed as
an orthogonal sum of the subspace defined by the divergence of the interior
functions of the $p$-th order $BDM_p$ basis and the constant function.
Therefore, the interior functions can be used to remove element-wise all higher
order terms except the constant in the divergence error of the finite element
solution of $B$-field. The constant terms from each element can be then easily
corrected using a first order H(div) basis globally. Numerical results for a
3-D magnetic induction equation show the effectiveness of the proposed method
in enforcing divergence-free condition of the magnetic field.
