Least-squares finite element methods for first-order formulations of the Poisson equation are not subject to the inf-sup condition and lead to stable solutions even when all variables are approximated by equal-order continuous finite element spaces. For such elements, one can also prove optimal convergence in the ``energy'' norm for all variables and optimal L^2 convergence for the scalar variable. However, showing optimal L^2 convergence for the flux has proven to be impossible without adding the redundant curl equation to the first- order system of partial differential equations. In fact, numerical evidence strongly suggests that nodal continuous flux approximations do not posses optimal L^2 accuracy. In this talk, we show that optimal L^2 error rates for the flux can be achieved without the curl constraint, provided that one uses the div-conforming family of Brezzi-Douglas-Marini or Brezzi-Douglas-Duran-Fortin elements. Then, we proceed to reveal an interesting connection between a least- squares finite element method involving div-conforming flux approximations and mixed finite element methods based on the classical Dirichlet and Kelvin principles. We show that such least- squares finite element methods can be obtained by approximating, through an L^2 projection, the Hodge operator that connects the Kelvin and Dirichlet principles. Our principal conclusion is that when implemented in this way, a least-squares finite element method combines the best computational properties of mixed-Galerkin finite element methods based on each of the classical principles.