  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.
