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