In this work we develop a novel method for magnetohydrodynamics on convex domains robust with respect to the Reynolds and Hartmann numbers. The starting point is a weak formulation of magnetohydrodynamics where an H1-based form is considered for all diffusive terms. Both scalar- and vector-valued variables are discretized using polynomial functions inside the elements and on the faces in the spirit of \cite{one,two}. Concerning vector-valued fields, face variables are taken in full polynomial spaces, whereas element variables are taken in the Raviart-Thomas-Nédélec space. This choice leads, in particular, to the pointwise satisfaction of the mass conservation and Gauss's law. Crucial attention has been put into designing the trilinear form, as its skew-symmetry properties are crucial in the proof of the convergence theorem. A thorough stability and convergence analysis for the space semi-discrete problem is carried out. Specifically, we prove energy error estimates of order $h^{k+1}$ for the version of the scheme corresponding to a polynomial degree $k \ge 0$. Using techniques inspired by \cite{three}, such estimates are additionally robust with respect to the pressure and semirobust with respect to the Reynolds and Hartmann numbers, meaning that the constants in the right-hand side remain bounded in the whole range of value for such dimensionless numbers provided that the velocity and magnetic fields are smooth enough.