  We present a virtual element method (VEM) for the numerical approximation of
the electromagnetics subsystem of the resistive magnetohydrodynamics (MHD)
model in two spatial dimensions. The major advantages of the virtual element
method include great flexibility of polygonal meshes and automatic
divergence-free constraint on the magnetic flux field. In this work, we
rigorously prove the well-posedness of the method and the solenoidal nature of
the discrete magnetic flux field. We also derive stability energy estimates.
The design of the method includes three choices for the construction of the
nodal mass matrix and criteria to more alternative. We present a set of
numerical experiments that independently validate theoretical results. The
numerical experiments include the convergence rate study, energy estimates and
verification of the divergence-free condition on the magnetic flux field. All
these numerical experiments have been performed on triangular, perturbed
quadrilateral and Voronoi meshes. Finally, We demonstrate the development of
the VEM method on a numerical model for the Hartmann flow.
