Finding the maximum of a function on the mesh

Hello to all of you FreeFem users,

After solving my system of PDE, I would like to compute the maximum value of my solution and the corresponding coordinates in my domain.

I found on another thread that I can do “u[].max” to find the maximum on the mesh (where u is my solution) but then how can I find the corresponding (x,y) coordinates ?

I know that the result of this procedure won’t be the true maximum if I use P2 elements because it gives the maximum value on vertices of the mesh but since I’m using a fine mesh it’s not really important, I don’t need a huge precision.

Thank you in advance.

mesh Th = square(20, 20, [x+0.2,y+0.2]);
fespace Vh(Th, P2);
Vh u = cos(6*pi*x)*tanh(y) + sin(2*pi*y)*tanh(x);
int i = u[].imax;
cout << i << " " << u[](i) << " " << u[].max << endl;
Vh cx = x, cy = y;
cout << cx[](i) << " " << cy[](i) << endl;