// file buildlayermesh.edp 
load "msh3" 
load "tetgen" 
// Test 1 

int C1=99, C2=98; // could be anything 

border C01(t=0,pi){ x=t;  y=0;      label=1;} 
border C02(t=0,2*pi){ x=pi; y=t;  label=1;} 
border C03(t=0,pi){ x=pi-t;  y=2*pi;    label=1;} 
border C04(t=0,2*pi){ x=0;    y=2*pi-t; label=1;} 

border C11(t=0,0.7){ x=0.5+t;  y=2.5;      label=C1;} 
border C12(t=0,2){ x=1.2;    y=2.5+t;  label=C1;} 
border C13(t=0,0.7){ x=1.2-t;  y=4.5;     label=C1;} 
border C14(t=0,2){ x=0.5;    y=4.5-t; label=C1;} 

border C21(t=0,0.7){ x= 2.3+t;     y=2.5;  label=C2;} 
border C22(t=0,2){        x=3;   y=2.5+t;  label=C2;} 
border C23(t=0,0.7){   x=3-t;     y=4.5;  label=C2;} 
border C24(t=0,2){       x=2.3;   y=4.5-t; label=C2;} 

mesh Th=buildmesh(    C01(10)+C02(10)+ C03(10)+C04(10) 
                    + C11(5)+C12(5)+C13(5)+C14(5) 
                    + C21(-5)+C22(-5)+C23(-5)+C24(-5)); 

mesh Ths=buildmesh(    C01(10)+C02(10)+ C03(10)+C04(10) 
                    + C11(5)+C12(5)+C13(5)+C14(5) ); 

// construction of a box with one hole and two regions 
func zmin=0.; 
func zmax=1.; 
int MaxLayer=10; 

int[int] r1=[0,41]; 
int[int] r2=[98,98,99,99,1,56]; 
int[int] r3=[4,12];    //  The triangles of upper surface generated by the triangle in the 2D region of mesh Th of label 4 as label 12. 
int[int] r4=[4,45];    //  The triangles of lower surface generated by the triangle in the 2D region of mesh Th of label 4 as label 45. 

mesh3 Th3=buildlayers(Th, MaxLayer, zbound=[zmin,zmax], region=r1, labelmid=r2, reffaceup = r3, reffacelow = r4 ); 

// Construction of a sphere with tetgen 
func XX1 = cos(y)*sin(x); 
func YY1 = sin(y)*sin(x); 
func ZZ1 = cos(x); 

real [int] domain = [0.,0.,0.,0,0.001]; 
string test="paACQ"; 
cout << endl; 
cout << "test=" << test << endl; 
mesh3 Th3sph=tetgtransfo(Ths,transfo=[XX1,YY1,ZZ1],switch=test,nbofregions=1,regionlist=domain);