1load "msh3"
2load "medit"
3searchMethod=1; // more safe seach algo .. (FH for PICHON ??)
4verbosity=1;
5real a=1, d=0.5, h=0.5;
6border b1(t=0.5,-0.5) {x=a*t; y=-a/2; label=1;};
7border b2(t=0.5,-0.5) {x=a/2; y=a*t; label=2;};
8border b3(t=0.5,-0.5) {x=a*t; y=a/2; label=3;};
9border b4(t=0.5,-0.5) {x=-a/2; y=a*t; label=4;};
10border i1(t=0,2*pi) {x=d/2*cos(t); y=-d/2*sin(t); label=7;};
11int nnb=7, nni=10;
12mesh Th=buildmesh(b1(-nnb)+b3(nnb)+b2(-nnb)+b4(nnb)+i1(nni));//, fixedborder=true);
13//Th=adaptmesh(Th,0.1,IsMetric=1,periodic=[[1,x],[3,x],[2,y],[4,y]]);
14int nz=3;
15{ // for cleanning  memory..
16int[int] old2new(0:Th.nv-1);
17fespace Vh2(Th,P1);
18Vh2 sorder=x+y;
19sort(sorder[],old2new);
20int[int]  new2old=old2new^-1;   // inverse the permuation
21//for(int i=0;i< Th.nv;++i) // so by hand.
22//  new2old[old2new[i]]=i;
23Th= change(Th,renumv=new2old);
24sorder[]=0:Th.nv-1;
25}
26{
27  fespace Vh2(Th,P1);
28  Vh2 nu;
29  nu[]=0:Th.nv-1;
30  plot(nu,cmm="nu=",wait=1);
31}
32int[int] rup=[0,5], rlow=[0,6], rmid=[1,1,2,2,3,3,4,4,7,7], rtet=[0,41];
33func zmin=0;
34func zmax=h;
35mesh3 Th3=buildlayers(Th, nz, zbound=[zmin,zmax],
36reftet=rtet,reffacemid=rmid, reffaceup=rup, reffacelow=rlow);
37for(int i=1;i<=6;++i)
38  cout << " int " << i << " :  " << int2d(Th3,i)(1.) << " " << int2d(Th3,i)(1./area) << endl;
39savemesh(Th3,"Th3.mesh");
40plot(Th3,wait=1);
41medit("Th3",Th3);
42
43fespace Vh(Th3,P2, periodic=[[1,x,z],[3,x,z],[2,y,z],[4,y,z],[5,x,y],[6,x,y]]);