A surface integration method for modeling of permanent magnets in three-dimensional space using edge finite elements is presented, where the value of the coercive magnetic force Hc is integrated only over the entire surface area of the permanent magnet. This allows direct computation of the equivalent source current values, which are afterwards assigned to the surface edges of the permanent magnet. The main advantages of the proposed method are: accurate results with less computational effort and improved convergence rate of the iterative solver. Verification of this surface integration method is carried by comparing the results with those obtained by the Biot–Savart law and the traditional current sheet method. Finally, a comparison between the numerical results obtained through the surface integration method and the measured results of a complex electromagnetic device with permanent magnets is given.