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