Thank you for your replies. I finally found the solution which was how I defined the boundary. I should have used meshtags as in How to compute flux at particular boundary facets?.
I do not need so far to compute the flux at every point of the boundary, only its integral, so I do not need the normal vector approximation. But thanks for pointing that out, it will certainly be useful someday.
1 Like