load "distance"
mesh Th=buildmesh("g.gmesh",nbvx=10000);
fespace Vh(Th,P1);
real[int] bb(4); 
boundingbox(Th,bb);
// b[0] xmin, xmax, ymin ymax 
real xc1 = bb[0]*0.6 + bb[1]*0.4;
real xc2 = bb[0]*0.2 + bb[1]*0.8;
real yc = bb[2]*0.3 + bb[3]*0.7;
real r2 = ((bb[3]-bb[2])*.1)^2;
cout << xc1 << " "<< xc2 << " " << yc << "r2 = " << r2 << endl; 
cout << bb << endl;
func f0 = ((x-xc1)^2 + (y-yc)^2- r2)*((x-xc2)^2 + (y-yc)^2 - r2);
Vh u = f0, v;

plot(u, dim=3, value=1,wait=1); 

verbosity=1; 


distance(Th,f0,v[]);

real[int] v0 = [-10,-1,-0.01,0,0.01,1,10];
plot(u,wait=1,viso=v0);
plot(v,wait=1,cmm="distance");
