In recent years there has been an increasing interest in the theoretical and experimental study of field responsive, functional composite materials. Magneto-active polymers (MAPs) are a special class of field responsive solids that comprise of a polymeric matrix with dispersed micro-sized magnetizable particles. Based on the magnetic properties of the underlying ferromagnetic filler particles, MAPs can be classified into two categories: (i) soft and (ii) hard MAPs. Soft MAPs comprising magnetically soft particles, e.g. carbonyl iron, exhibit negligible hysteresis loss and demagnetize completely after the removal of the external magnetic field which in consequence leads to reversible deformation mechanisms. In contrast, NdFeB particle-filled hard MAPs exhibit distinct nonlinear, dissipative material behavior, such as the characteristic magnetic and ”butterfly” field-induced strain hysteresis. In the present work, we present a comprehensive variational-based modeling framework for hard MAPs including the response of soft MAPs as a limiting case. We outline ingredients of the constitutive theory based on the concept of generalized standard materials, that necessitates suitable definitions of (i) the total energy density function and (ii) the dissipation potential [1]. Key idea of the constitutive approach is an additive split of the material part of the total energy density function into three contributions associated with (i) an elastic ground stress, (ii) the magnetization and (iii) a magnetically induced mechanical stress, respectively. We propose suitable constitutive functions in an energy-based setting that allow to accurately capture the highly nonlinear material behavior of hard MAPs with stochastic microstructures. Subsequently, the two constitutive functions are embedded in a time-discrete minimization principle that is supplemented by a conforming finite element method. Furthermore, we discuss polyconvexity of the constitutive model which implies quasiconvexity and rank-one convexity and thus ensures material stability. The performance of the developed variational-based modeling framework is demonstrated by solving some application-oriented boundary value problems. The main emphasis of the numerical studies lies on the investigation of the magnetostrictive effect of hard MAPs at the macroscale level as well as on the in depth analysis of pre-magnetized beam structures undergoing large deformations.