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