The “obvious” approach of directly evaluating \nabla u\cdot\mathbf{n} on the boundary turns out not to be the most accurate way to post-process fluxes from finite element solutions. You can look at my responses here and here for some discussion on extracting boundary fluxes that exactly satisfy the conservation law (including code examples and references).
1 Like