We introduce an accurate cell-centered method for modeling Darcy flow on general quadrilateral, hexahedral, and simplicial grids. We refer to these discretizations as the multipoint flux mixed finite element (MFMFE) method. The MFMFE method is locally conservative with continuous fluxes and can be viewed within a variational framework as a mixed finite element method with special approximating spaces and quadrature rules. We study two versions of the method: with a symmetric quadrature rule on smooth grids and a non-symmetric quadrature rule on rough grids. The framework allows for handling hexahedral grids with non-planar faces defined via trilinear mappings from the reference cube. Moreover, the MFMFE method allows for local elimination of the velocity, which leads to a cell-centered pressure system. Theoretical and numerical results demonstrate first order convergence on rough grids. Second order superconvergence is observed on smooth grids. We also discuss a new splitting scheme for modeling multiphase flows that can treat higher order transport discretizations for saturations. We apply the MFMFE method to obtain physically consistent approximations to the velocity and a reference pressure on quadrilateral or hexahedral grids, and a discontinuous Galerkin method for saturations. For higher order saturations, we propose an efficient postprocessing technique that gives accurate velocities in the interior of the gridblocks. Computational results are provided for flow in highly heterogenous reservoirs, including different capillary pressures arising from different rock types.