A three dimensional adaptive finite element refinement algorithm is developed for non-linear magnetostatic field problems. In the method, the edge elements are used for finite element formulation, and the local error in each element is estimated from the fact that the tangential components of magnetic field intensity and the normal components of magnetic flux density should be continuous at the interface of the two adjacent elements. Based on the estimated error, the elements which have big error are divided into several elements using bisection method. The effectiveness of the developed algorithm is proved through numerical examples.