I have a very simple question. I am doing flux calculations and want to calculate line integral on the “outlet” boundary (the right boundary) of a square mesh. I know my solution i.e u_h is correct because it converges but it seems like I am doing something wrong with flux calculation since it is not giving me the results I expect.

The analytical flux I expect is grad(u_analy)*n and I choose u_analytical to be (x^2+y^2,0) in two dimensions. Since I am in right boundary of a square mesh my n is (1,0). So I expect flux analytical to be (2x , 0). The square has dimensions of [-1,1]^2 so the analytical result is 0. But mu numerical result is far from that… I am doing the boundary integration below…

With the same logic, on the circular boundary defined by level set function x^2+y^2-0.5^2=0 (using the same analytical solution) i calculate n to be <2x/sqrt(4x^2+4y^2) , 2y/sqrt(4x^2+4y^2)> and I calculate grad(u_analy)=<2x , 2y>
that makes analytical flux around the circle which is grad(u_analy).n to be
<4x^2/sqrt(4x^2+4y^2) + 4y^2/sqrt(4x^2+4y^2)> = 1. Am I also wrong here ?

Let the domain [-1,1]^2 and u = (x^2+y^2,0). Then grad(u)= (2x, 2y; 0, 0). The right boundary is given by the set { (1,y) : -1<=y<=1} and the according normal vector is n=(1,0). Therefore, grad(u)*n = (2x, 0) on the right boundary. Further x=1 on the right boundary => grad(u)*n = (2, 0) on the right boundary. Integrating this flux over the right boundary (it has length 2) yields int(2, y=-1…1) = 4.
This is what I also observe numerically.

Let the domain given by the level-set function x^2+y^2-0.5^2=0. The normal vector is then given by n= (2x,2y) (simplifying your expression by using the level-set equation). Thus, grad(u)n = (4x^2 + 4y^2, 0) = (1,0). Integrating over the boundary of the circle yields int(1) = 2r*pi= pi.
The numerical result is also near pi.

Attached you find the according python-code for both examples.