  A new weak Galerkin (WG) method is introduced and analyzed for the second
order elliptic equation formulated as a system of two first order linear
equations. This method, called WG-MFEM, is designed by using discontinuous
piecewise polynomials on finite element partitions with arbitrary shape of
polygons/polyhedra. The WG-MFEM is capable of providing very accurate numerical
approximations for both the primary and flux variables. Allowing the use of
discontinuous approximating functions on arbitrary shape of polygons/polyhedra
makes the method highly flexible in practical computation. Optimal order error
estimates in both discrete $H^1$ and $L^2$ norms are established for the
corresponding weak Galerkin mixed finite element solutions.
