The last decade has witnessed the emergence of magneto-active polymers (MAPs) as one of the most advanced multi-functional soft composites. Depending on the magnetisation mechanisms and responsive behaviour, MAPs are mainly classified into two groups: i) hard magnetic MAPs in which a large residual magnetic flux density sustains even after the removal of the external magnetic field, and ii) soft magnetic MAPs where the magnetisation of the filler particles disappear upon the removal of the external magnetic field. Polymeric materials are widely treated as fully incompressible solids that require special numerical treatment to solve the associated boundary value problem. Furthermore, both soft and hard magnetic particles-filled soft polymers are inherently viscoelastic. Therefore, the aim of this paper is to devise a unified finite element method-based numerical framework for magneto-mechanically coupled systems that can work for compressible and fully incompressible materials and from hard to soft MAPs, including the effects of the time-dependent viscoelastic behaviour of the underlying matrix. First, variational formulations for the uncoupled problem for hard MAPs and the coupled problem for soft MAPs are derived. The weak forms are then discretised with higher-order Bézier elements while the evolution equation for internal variables in viscoelastic models is solved using the generalised-alpha time integration scheme, which is implicit and second-order accurate. Finally, a series of experimentally-driven boundary value problems consisting of the beam and robotic gripper models are solved in magneto-mechanically coupled settings, demonstrating the versatility of the proposed numerical framework. The effect of viscoelastic material parameters on the response characteristics of MAPs under coupled magnetomechanical loading is also studied.