In this review we discuss self-consistent methods to calculate the global structure of strongly magnetised neutron stars within the general-relativistic framework. We outline why solutions in spherical symmetry cannot be applied to strongly magnetised compact stars, and elaborate on a consistent formalism to compute rotating magnetised neutron star models. We also discuss an application of the above full numerical solution for studying the influence of strong magnetic fields on the radius and crust thickness of magnetars. The above technique is also applied to construct a "universal" magnetic field profile inside the neutron star, that may be useful for studies in nuclear physics. The methodology developed here is particularly useful to interpret multi-messenger astrophysical data of strongly magnetised neutron stars.