Hard magnetorheological elastomers ( h -MREs) are essentially two phase composites comprising permanently magnetizable metallic inclusions suspended in a soft elastomeric matrix. This work provides a thermodynamically consistent, microstructurally-guided modeling framework for isotropic, incompressible h -MREs. Energy dissipates in such hard-magnetic composites primarily via ferromagnetic hysteresis in the underlying hard-magnetic particles. The proposed constitutive model is thus developed following the generalized standard materials framework, which necessitates suitable definitions of the energy density and the dissipation potential. Moreover, the proposed model is designed to recover several well-known homogenization results (and bounds) in the purely mechanical and purely magnetic limiting cases. The magneto–mechanical coupling response of the model, in turn, is calibrated with the aid of numerical homogenization estimates under symmetric cyclic loading. The performance of the model is then probed against several other numerical homogenization estimates considering various magneto–mechanical loading paths other than the calibration loading path. Very good agreement between the macroscopic model and the numerical homogenization estimates is observed, especially for stiff to moderately-soft matrix materials. An important outcome of the numerical simulations is the independence of the current magnetization to the stretch part of the deformation gradient. This is taken into account in the model by considering an only rotation-dependent remanent magnetic field as an internal variable. We further show that there is no need for an additional mechanical internal variable. Finally, the model is employed to solve macroscopic boundary value problems involving slender h -MRE structures and the results match excellently with experimental data from literature. Crucial differences are found between uniformly and non-uniformly pre-magnetized h -MREs in terms of their pre-magnetization and the associated self-fields.