FEM simulations are a standard step in the design of accelerator magnets. It is custom for accelerator applications to characterize the field quality in terms of field expansion coefficients. Expansion coefficients are usually calculated by means of a Fourier transform of the local FEM solution evaluated at points on a circle (2D) or cylinder (3D case). The accuracy of the coefficients calculated this way depends strongly on the FEM mesh configuration and simple refinement of the mesh does not always improve accuracy. The accuracy of the expansion coefficients calculation can be improved by using the data on the magnetization of elements in the magnet yoke, obtained in the solution, instead of using directly the local solution. Since currents and the yoke magnetization are the only sources of the field, with these data the field expansion coefficients can be calculated at any remote point. We derive closed forms for calculating expansion coefficients and implemented these results in the ANSYS® add-on. Results for a case study are presented, which demonstrate that expansion coefficients can be calculated with good accuracy even for a rather coarse mesh.