Magnetic hysteresis affects the performance of electromagnetic devices, e.g., motors, generators, and transformers. However, due to complex, non-linear, and multi-valued nature of this phenomenon, its accurate coupling to Finite Element Analysis (FEA) of these devices has been always a challenging task. A novel approach has been presented in this paper for linking FEA to the Preisach model, which is known as the most accurate hysteresis model. An individual Preisach module has been considered for each field component of each element of the hysteresis material mesh. Hysteresis characteristics between each two successive time steps have been linearly approximated. An iterative algorithm has been proposed for obtaining field distributions along with parameters of these lines, simultaneously. The proposed method has been applied to a general magnetic circuit to predict its behavior over a given time span. Space distributions of flux density at some time steps, time variations of flux density and field intensity for one element, induced voltage, and hysteresis characteristics for some elements have been obtained. In contrast with most previous works, approach of this paper could reflect the details of hysteresis phenomenon, including minor loops, into the FEA. Also, it is applicable to problems with non-uniform and rotating field distributions.