Computing normal heat flux density

Dear all,

I’m currently working on a thermic problem (Laplace equation) and i want to compute the normal heat flux density along the borders.
I tried to compute it with a variationnal form, like this:
qx = -khdx(u); // Fourier’s Law
qy = -kh
varf vphi(f,g) = on(0, 1, 2, 3, 4, 5, 6, 7, 8 , 9, 10, 11, 12, f = qxN.x + qyN.y);
phi[] = vphi(0., Vh, tgv = 1);
where u is the temperature field. It seems fine on simple problem but on more complex ones (complex geometry with mixed boundary condition) i can’t find zero when i integrate it even with a fine mesh (global outflux = 0 with no sinks or sources in the domain).

The problem seems to occur where a Dirichlet condition meets a Robin condition ( Newton’s Law of cooling: -k.grad(u).n = h ( u - us) ), i have a singularity, my temperature field is fine but my heat flux density is infinite.
But i thought that this problem was local and wouldn’t influence too much the global outflux.

I tried to refine the mesh with the residual error indicator presented in the doc. It is maximum around the singularity as expected but even after 10 iterations i cannot get it to an acceptable value.

Have you any idea how to fix this ?

I apologize in advance, i’m not really experienced and i’m probably just facing standard issue with the finite element method.

By the way, how does freefem calculate the partial derivatives on borders ? I tried, with no results to calculate the normal heat flux density by hand without using the functions dx() and dy()…

Best regards,
Gayrard Hugo