Compared with standard Galerkin finite element methods, mixed methods for second-order elliptic problems give readily available flux approximation, but in general at the expense of having to deal with a more complicated discrete system. This is especially true when conforming elements are involved. Hence it is advantageous to consider a direct method when finding fluxes is just a small part of the overall modeling processes. The purpose of this article is to introduce a direct method combining the standard Galerkin Q1 conforming method with a cheap local flux recovery formula. The approximate flux resides in the lowest order Raviart-Thomas space and retains local conservation property at the cluster level. A cluster is made up of at most four quadrilaterals
