mesh Th=square(100,100,[x*2-1,y*2-1]);
load "isoline"
fespace Vh(Th,P1);
fespace Ph(Th,P0); 
func rr= sqrt(x^2+y^2);
int k =5;
Vh u= sin(k*pi*x)*sin(k*pi*y);//rr*(1-rr);
plot(u, wait=1);
Ph r; 

verbosity=3;
int[int] lm=findalllocalmin(Th,u[],r[],convex=1,eps=0);
plot(r,u,fill=1,wait=1);
// lm.n == number of min
// r  bassin number 
// lm[int(r[k])] is the minimal vertex of the bassin associated 
//cout << r[]<<endl; 

cout << lm << endl; 
