1load "distance"
2load "msh3"
3load "medit"
4mesh3 Th=cube(10,10,10);
5fespace Vh(Th,P1);
6Vh u = z;//-0.501;
7func fxy=max(abs(x-y),abs(-1.+x+y));
8func fxz=max(abs(x-z),abs(-1.+x+z));
9func fyz=max(abs(y-z),abs(-1.+y+z));
10func fxyz = max(max(fxy,fxz),fyz);
11func ff=max(max(abs(x-0.5),abs(y-0.5)),abs(z-0.5));
12u = 0.25-ff;
13real[int] v0 = [0.001];
14plot(u, dim=3, value=1,viso=v0,wait=1);
15Vh v=0,zz=z;
16verbosity=3;
17//cout << v[] << endl;
18
19distance(Th,u,v[],distmax=1);
20
21
22plot(u,wait=1,viso=v0);
23plot(v,wait=1,cmm="distance",nbiso=10);
24cout << " vmin/ vmax  " << v[].min << " "<< v[].max << endl;
25verbosity=0;
26/*
27for( int k=0; k< Th.nt; ++k)
28for( int j=0; j< 4; ++j)
29{
30	int i= Vh(k,j);
31if(abs(v[][i]-zz[][i]) > 0.001)
32cout << " i "<< i << " k = "<< k << " " << j << endl;
33}
34//medit("Th-v",Th,v);
35*/