load "distance"
mesh Th;
Th=square(100,100);
if(1)
{
	mesh Th2;
	mesh3 Th3;
	int ltube=1,lup=3,lout=4,lin =2;
	int l3d = 5;
	real D3d=1; 
	{
	int Nbtube= 4; 
    
	real Dx=1,Dy=1,La=3*Dx,Lh=5*Dy;

	int Ndx=19,Ndy=19, Nh=(Lh-Dy)*Ndy/Dy,Na=(La-Dx)*Ndx/Dx ;     

	mesh Th0;
	{
	    int[int] ll0=[1,0,0,1];
	    int[int] llv=[0,1,2,1];
	    int[int] llvv=[0,1,4,1];
	    int[int] llh=[1,3,1,0];
	 Th0 = square(Ndx,Ndy,[x*Dx,y*Dy],label=ll0); 
	mesh Thv = square(Ndx,Nh,[x*Dx,Dy+y*(Lh-Dy)],label=llv); 
	mesh Thvv = square(Ndx,Ndy,[x*Dx,Lh+y*(Dy)],label=llvv); 
	mesh Thh = square(Na,Ndy,[Dx+x*(La-Dx),y*Dy],label=llh); 

	Th=Th0+Thv+Thh;
	Th0=Th;
	Th=Th+Thvv;
	plot(Th,wait=1);
	}

	for(int i=1; i< Nbtube; ++i)
	   Th=Th + movemesh(Th0,[x+La*i,y]);
	   int[int] llabel=[1,ltube,2,lup,3,lout,4,lin] ;
	 Th=change(Th,label=llabel,rmInternalEdges=1);  
	}
	plot(Th,wait=1);
}
fespace Vh(Th,P1);
Vh u = y-0.500; 
u = 0.400-max(abs(x-y),abs(-1.+x+y));
plot(u, dim=3, value=1); 
Vh v=u;
//verbosity=100; 
//cout << v[] << endl; 

distance(Th,u,v[],distmax=100);

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