load "msh3"
load "medit"
// Parameters
searchMethod=1; // More safe seach algo
real a = 1, d = 0.5, h = 0.5;
int nnb = 7, nni = 10;
int nz = 3;
func zmin = 0;
func zmax = h;
// Mesh 2D
border b1(t=0.5, -0.5){x=a*t; y=-a/2; label=1;}
border b2(t=0.5, -0.5){x=a/2; y=a*t; label=2;}
border b3(t=0.5, -0.5){x=a*t; y=a/2; label=3;}
border b4(t=0.5, -0.5){x=-a/2; y=a*t; label=4;}
border i1(t=0, 2.*pi){x=d/2*cos(t); y=-d/2*sin(t); label=7;}
mesh Th = buildmesh(b1(-nnb) + b3(nnb) + b2(-nnb) + b4(nnb) + i1(nni));
{ // Cleaning the memory correctly
int[int] old2new(0:Th.nv-1);
fespace Vh2(Th, P1);
Vh2 sorder = x + y;
sort(sorder[], old2new);
int[int] new2old = old2new^-1; // Inverse permutation
Th = change(Th, renumv=new2old);
sorder[] = 0:Th.nv-1;
}
{
fespace Vh2(Th, P1);
Vh2 nu;
nu[] = 0:Th.nv-1;
plot(nu, cmm="nu=", wait=true);
}
// Mesh 3D
int[int] rup = [0, 5], rlow = [0, 6], rmid = [1, 1, 2, 2, 3, 3, 4, 4, 7, 7], rtet = [0, 41];
mesh3 Th3 = buildlayers(Th, nz, zbound=[zmin, zmax],
reftet=rtet, reffacemid=rmid, reffaceup=rup, reffacelow=rlow);
for(int i = 1; i <= 6; ++i)
cout << " int " << i << " : " << int2d(Th3,i)(1.) << " " << int2d(Th3,i)(1./area) << endl;
plot(Th3, wait=true);
medit("Th3", Th3);
fespace Vh(Th3, P2, periodic=[[1, x, z], [3, x, z], [2, y, z], [4, y, z], [5, x, y], [6, x, y]]);