A mixed variational formulation for the solution of three dimensional magnetostatic problems, based on both the magnetic field H and the magnetic induction field B, has been developed. The approach is based upon the minimisation of the residual in the fulfilment of the constitutive relationship, imposing Maxwell's equations as constraints. A finite element discretisation has been realised by using face and edge elements to interpolate B and Fl., respectively, in order to guarantee the physical continuity properties. The mathematical formulation, its code development and computational results are presented and discussed in detail.