We develop a class of mixed finite element schemes for stationary magnetohydrodynamics (MHD) models, using the magnetic field $\mathbfit {B} $ and the current density $\mathbfit {j} $ as discretization variables. We show that Gauss’s law for the magnetic field, namely $\nabla\cdot\mathbfit {B}= 0$, and the energy law for the entire system are exactly preserved in the finite element schemes. Based on some new basic estimates for finite elements, we show that the new finite element scheme is well-posed. Furthermore, we show the existence of solutions to the nonlinear problems and the convergence of the Picard iterations and the finite element methods under some conditions. References