1load "msh3" load "tetgen" load "mshmet" load "medit"
2//build initial mesh
3int nn  = 6;
4int[int] l1111=[1,1,1,1],l01=[0,1],l11=[1,1];//   label numbering
5mesh3 Th3=buildlayers(square(nn,nn,region=0,label=l1111),
6      nn,  zbound=[0,1],  labelmid=l11,   labelup = l01,  labeldown = l01);
7Th3 = trunc(Th3,(x<0.5) | (y < 0.5) | (z < 0.5) ,label=1);// remove the $]0.5,1[^3 cube$
8//end of build initial mesh
9fespace Vh(Th3,P1);
10Vh u,v,usol;
11
12macro Grad(u) [dx(u),dy(u),dz(u)] // EOM
13
14problem Poisson(u,v,solver=CG) = int3d(Th3)( Grad(u)'*Grad(v) )  // ') for emacs
15  -int3d(Th3)( 1*v ) + on(1,u=0);
16
17real errm=1e-2;// level of error
18
19for(int ii=0; ii<5; ii++)
20{
21  Poisson;
22  cout <<" u min, max = " <<  u[].min << " "<< u[].max << endl;
23  Vh h ;
24  h[]=mshmet(Th3,u,normalization=1,aniso=0,nbregul=1,hmin=1e-3,hmax=0.3,err=errm);//loptions=MSHloptions,doptions=MSHdoptions);
25  cout <<" h min, max = " <<  h[].min << " "<< h[].max << " " << h[].n << " " << Th3.nv << endl;
26  // FFCS: add 3D view parameters
27  plot(u,wait=1,fill=0,boundary=0,CutPlane=0,ShowMeshes=1,LabelColors=0);
28  errm*= 0.8;// change the level of error
29  cout << " Th3" << Th3.nv < " " << Th3.nt << endl;
30  Th3=tetgreconstruction(Th3,switch="raAQ",sizeofvolume=h*h*h/6.);
31  medit("U-adap-iso-"+ii,Th3,u,wait=1);
32}
33
34
35
36