Summary 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 nonsymmetric quadrature rule on rough grids. The framework allows for handling hexahedral grids with nonplanar faces defined by 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 post-processing technique that gives accurate velocities in the interior of the gridblocks. Computational results are provided for flow in highly heterogeneous reservoirs, including different capillary pressures arising from different rock types.