Compute flux from weak form

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:

Thank you! This makes a lot clear.