A new least-squares finite element method is developed for the curl-div magnetostatic problem in Lipschitz and multiply connected domains filled with anisotropic nonhomogeneous materials. In order to deal with possibly low regularity of the solution, local $L^2$ projectors are introduced to standard least-squares formulation and applied to both curl and div operators. Coercivity is established by adding suitable mesh-dependent bilinear terms. As a result, any continuous finite elements (lower-order elements are enriched with suitable bubbles) can be employed. A desirable error bound is obtained: $||{\bf u}-{\bf u}_h||_0\le C\,||{\bf u}-\tilde{{\bf u}}||_0$, where ${\bf u}_h$ and $\tilde{{\bf u}}$ are the finite element approximation and the finite element interpolant of the exact solution ${\bf u}$, respectively. Numerical tests confirm the theoretical results.