1load "distance"
2mesh Th=buildmesh("g.gmesh",nbvx=10000);
3fespace Vh(Th,P1);
4real[int] bb(4);
5boundingbox(Th,bb);
6// b[0] xmin, xmax, ymin ymax
7real xc1 = bb[0]*0.6 + bb[1]*0.4;
8real xc2 = bb[0]*0.2 + bb[1]*0.8;
9real yc = bb[2]*0.3 + bb[3]*0.7;
10real r2 = ((bb[3]-bb[2])*.1)^2;
11cout << xc1 << " "<< xc2 << " " << yc << "r2 = " << r2 << endl;
12cout << bb << endl;
13func f0 = ((x-xc1)^2 + (y-yc)^2- r2)*((x-xc2)^2 + (y-yc)^2 - r2);
14Vh u = f0, v;
15
16plot(u, dim=3, value=1,wait=1);
17
18verbosity=1;
19
20
21distance(Th,f0,v[]);
22
23real[int] v0 = [-10,-1,-0.01,0,0.01,1,10];
24plot(u,wait=1,viso=v0);
25plot(v,wait=1,cmm="distance");
26