Branched discussion from Using another solver to modify OpenFOAM-Nutils CHT tutorial
That’s a very classical question I had 2 years ago as well. That’s the good way of computing fluxes. Directly differentiating the shape functions is the bad way.
You find a good explanation on page 3 of
Toselli, Andrea, and Olof Widlund. Domain decomposition methods-algorithms and theory. Vol. 34. Springer Science & Business Media, 2006."
More details also here: http://dx.doi.org/https://doi.org/10.1115/1.4005187