Hi everyone!

I tried to implement the indicator function of a rectangle Q as follows

// Q = [xq1,xq2] x [yq1,yq2], xq1<xq2, yq1<yq2

func real chiQ(real xq1,real xq2,real yq1,real yq2){

if(x>=xq1 & x<=xq2 & y>=yq1 & y<=yq2) return 1;

else return 0;

}

but I think I am wrong with something.

Indeed, by computing the integral of the indicator function of Q = [0.25,0.75] x [0.25,0.75], I should obtain the area of Q, that is 0.25, but I got results, such as 0.26, 0.22, varying based on the mesh on which I’m computing the integral.

Here the way I compute this integral:

cout<<int2d(Th) (chiQ(0.25,0.75,0.25,0.75))<<endl;

Do you have any idea about why?