// example41.edp from cube.idp and Example 5.20, p. 124
load "msh3"
load "medit"

verbosity = 3;
int nx=10, ny=10, nz=10; // number of steps in each direction

real x0=0., x1=1.;
real y0=0., y1=1.;
real z0=0., z1=1.;

// build one square
// area = 0
mesh Thx = square(nx, ny, [x0+(x1-x0)*x, y0+(y1-y0)*y] );
  
// renumber regions 
// 2D: bottom = 1, right = 2, top = 3, left = 4
// Following gives pairs for top/bottom: [2D number, 3D number]
// 4 pairs for side faces
int[int] rup=[0, 6],  rdown=[0, 5], 
    rside=[1, 3,  2, 2,  3, 4,  4, 1];

mesh3 Th=buildlayers(Thx, nz, zbound=[z0,z1], 
	       labelmid = rside, labelup = rup, labeldown = rdown);

medit("Cube",Th);
plot(Th, wait=true);