In finite volume method one is using cell averages. In nonlinear conservation laws discontinuities can be created in the solution process. How to compute the flux when the flux contains a gradient and the Taylor approximation is not valid due to discontinuities?
Non-degenerate diffusion will prevent discontinuities in the true solution, but those may be poorly resolved. I will assume that your question is about how to evaluate diffusive fluxes of the form $f(u,\nabla u,p) = \kappa(u,\nabla u,p)\nabla u$ where $\kappa$ is a $d\times d$ SPD tensor and $p$ represents some spatially distributed parameters that are cell-centered and discontinuous. There are a few approaches.
Take the centroids $c_0$ and $c_1$ of the two cells and compute $\nabla u\cdot(c_1 - c_0) = u(c_1) - u(c_0)$. If the mesh is orthogonal (with face quadrature point on line segment connecting cell centroids and normal with aligned normal) and the coefficients are constant and isotropic, this is a consistent second order accurate flux (i.e., superconvergent at centroids and face centers). Unfortunately, it loses accuracy and has stability problems on general meshes and for variable coefficients. Two-point fluxes are often used because they are simple and compact, but it is a fragile methodology since the grids regularity and constant coefficient assumptions are so frequently violated.
The underlying principle is to reconstruct a gradient on each cell, typically with a weighted least squares procedure, and use that gradient evaluated at quadrature points (face centers) to define a flux. This procedure works simply for nonlinear fluxes. There are many variants that you can find in the literature, but you may want to look at these notes for background and the most commonly used variants. The downside of this approach is that the stencil is larger than strictly necessary and the resulting linear systems are generally not symmetric or M-matrices (thus allowing oscillations).
Reinterpretation as mixed finite element methods
Mixed finite element methods explicitly discretize the flux functions and work with inherently discontinuous (perhaps tensor-valued) coefficients and solution variables. The downside of mixed FEM is that there are more degrees of freedom and the resulting system of equations is a saddle point problem (which makes it more difficult to solve). After a special choice of quadrature, however, the "flux" degrees of freedom can be eliminated resulting in a (usually symmetric) positive definite system for the cell centered variables. An example of this is detailed in Wheeler and Yotov (2006) which uses an underlying BDM-1 space (linear discontinuous flux functions on edges, constant pressure functions on cells) and a cell quadrature at cell corners which makes the mass matrix block diagonal (with one block per corner). After eliminating the fluxes, the resulting system is SPD and has a minimal stencil (only connecting cells that share a corner). It is second order superconvergent for tensor coefficients and general grids, though the fluxes are only first order accurate for irregular grids.
Note that in addition to the relationship between finite volume and mixed finite element, there is a duality between mixed finite element methods and mimetic finite difference methods. In all cases I know of, it is possible to choose basis functions and quadratures to make the methods coincide, though one or the other may be more natural for certain types of analysis or for implementation in a given software package.