//  created with md2edp examples/examples/dist-projection.md
load "msh3"
border C(t=0,2*pi){ x= cos(t); y=sin(t); label=1; region =1;}
meshL Lh =buildmeshL(C(20));
mesh Th= square(20,20,[(x-0.5)*3,(y-0.5)*3]);
fespace Vh(Th,P1);
Vh d=dist(Lh);
Vh Ux = projection(Lh).x-x;
Vh Uy = projection(Lh).y-y;
plot(d, [Ux,Uy],Lh, wait=1);
int nu;
R3 ph;
mesh Th2= square(10,10,[(x-0.5)*3,(y-1.06)*3]);
Ux = projection(Th2,nu=nu,Phat=ph).x-x;
Uy = projection(Th2,nu=nu,Phat=ph).y-y;
x=0;y=0;z=0;
 cout << projection(Th2,nu=nu,Phat=ph) <<", " << nu <<","<< ph << endl; 
 cout << Th2[nu][0].x << endl;
 cout << Th2[nu][0].y  << endl;
 cout << Th2[nu][1].x << endl;
 cout << Th2[nu][1].y  << endl;
 cout << Th2[nu][2].x << endl;
 cout << Th2[nu][2].y  << endl;
plot([Ux,Uy],Th2, wait=1);
load "msh3";
//meshL Lh2 = extract(Th2);
meshS ThS= square3(10,10,[(x-0.5)*3,(y-1.06)*3,(x+y)/2]);
mesh3 Th3=cube(10,10,10,[(x-0.5)*3,(y-1.06)*3,z*2-1]);
fespace Uh(Th3,P1);
Uh d3=signeddist(ThS);
plot(d3,wait=1);
