I have a Galerkin solution for a heat equation
$$ u_t = \Delta u + f $$
with Dirichlet conditions $$ u=0, \qquad x \in \partial\Omega $$
The time discretization is done using a BDF scheme. How can I accurately and efficiently compute the heat flux on the boundary $$ \sigma = \frac{\partial u}{\partial n}, \qquad x \in \partial\Omega $$ from the Galerkin solution $u_h$.