// 
// level 0
//
//   On block n side walls are 1n-2n-3n-4n
//
//  Bottom lowest block label 100
//  Top lowest block label 200
//
x0=-.5;x1=0.5;y0=0;y1=1;
mesh S0=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
int[int] rdowng=[0,100],rupg=[0,200],rtet=[0,0];
int[int] rmidg=[1,10,2,20,3,30,4,40];

mesh3 B0=buildlayers(S0,nn, zbound=[0,1],region=rtet, labelup=rupg, labeldown=rdowng,labelmid=rmidg);

//
// TRIANGULACION
//
load "msh3"
load "medit"  
//
int nn;
cout<<" MAKE THERMAL FIN TRIANGULATION.... GIVE n  "<<endl;
cin>>nn;

ofstream fileThsize("Meshsize.txt");
fileThsize<<nn;
border OX(t=1,3){x=t;y=0;z=0;}
border OY(t=-1,1){x=2;y=t;z=0;}
//border OZ(t=0,3){x=2;y=0;z=t;}
//
// Define the basis of the block coordinates
//
real x0,x1,y0,y1;
// 
// level 0
//
//   On block n side walls are 1n-2n-3n-4n
//
//  Bottom lowest block label 100
//  Top lowest block label 200
//
x0=-.5;x1=0.5;y0=0;y1=1;
mesh S0=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
int[int] rdowng=[0,100],rupg=[0,200],rtet=[0,0];
int[int] rmidg=[1,10,2,20,3,30,4,40];

mesh3 B0=buildlayers(S0,nn, zbound=[0,1],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);


//Level 1
//   Block from level 0 top label is 200
// Bottom of structure is 199-200-201 (200 must macht)
// Top of structure is 299-300-301
//
// Block 1
//
x0=-1.5;x1=-0.5;y0=0;y1=1;
mesh S1=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,199];rupg=[0,299];rtet=[0,1];
rmidg=[1,11,2,21,3,31,4,41];

mesh3 B1=buildlayers(S1,nn, zbound=[1,2],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);
// Block 2         
x0=-0.5;x1=0.5;y0=0;y1=1;
mesh S2=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,200];rupg=[0,300];rtet=[0,2];
rmidg=[1,12,2,22,3,32,4,42];

mesh3 B2=buildlayers(S2,nn, zbound=[1,2],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);
// Block 3         
x0=0.5;x1=1.5;y0=0;y1=1;

mesh S3=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,201];rupg=[0,301];rtet=[0,3];
rmidg=[1,13,2,23,3,33,4,43];

mesh3 B3=buildlayers(S3,nn, zbound=[1,2],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);  

         //
// Level 2: one block
// 
// Block from level 1 top of structure is 299-300-301
// Bottom of block on level 2 is now 300 (300 must macht)
// Top of block is 400
//       
x0=-0.5;x1=0.5;y0=0;y1=1;
mesh S4=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,300];rupg=[0,400];rtet=[0,4];

rmidg=[1,14,2,24,3,34,4,44];               
mesh3 B4=buildlayers(S4,nn, zbound=[2,3],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);

//Level 3
// top from Level 2 is 400
//    bottom of structure is 399-400-401
//   First level of structure is 499-500-501
// 
// Block 5
//
x0=-1.5;x1=-0.5;y0=0;y1=1;
mesh S5=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,399];rupg=[0,499];rtet=[0,5];
rmidg=[1,15,2,25,3,35,4,45];

mesh3 B5=buildlayers(S5,nn, zbound=[3,4],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);
// Block 6         
x0=-0.5;x1=0.5;y0=0;y1=1;
mesh S6=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,400];rupg=[0,500];rtet=[0,6];
rmidg=[1,16,2,26,3,36,4,46];

mesh3 B6=buildlayers(S6,nn, zbound=[3,4],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);
// Block 7         
x0=0.5;x1=1.5;y0=0;y1=1;

mesh S7=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,401];rupg=[0,501];rtet=[0,7];
rmidg=[1,17,2,27,3,37,4,47];
mesh3 B7=buildlayers(S7,nn, zbound=[3,4],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);
//
// Level 4: one block
// 
// Block from level 3 top of structure is 500
//    bottom of structure is 500
//   Top of structure is 600
// 
//       
x0=-0.5;x1=0.5;y0=0;y1=1;
mesh S8=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,500];rupg=[0,600];rtet=[0,8];

rmidg=[1,18,2,28,3,38,4,48];               
mesh3 B8=buildlayers(S8,nn, zbound=[4,5],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);

//Level 5
// top from Level 2 is 600
//    bottom of structure is 599-600-601
//   Top of structure is 699-700-701
// 
// Block 9
//
x0=-1.5;x1=-0.5;y0=0;y1=1;
mesh S9=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,599];rupg=[0,699];rtet=[0,9];
rmidg=[1,19,2,29,3,39,4,49];

mesh3 B9=buildlayers(S9,nn, zbound=[5,6],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);
// Block 10         
x0=-0.5;x1=0.5;y0=0;y1=1;
mesh S10=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,600];rupg=[0,700];rtet=[0,10];
rmidg=[1,110,2,210,3,310,4,410];

mesh3 B10=buildlayers(S10,nn, zbound=[5,6],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);
// Block 11         
x0=0.5;x1=1.5;y0=0;y1=1;

mesh S11=square(nn,nn,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

rdowng=[0,601];rupg=[0,701];rtet=[0,11];
rmidg=[1,111,2,211,3,311,4,411];
mesh3 B11=buildlayers(S11,nn, zbound=[5,6],region=rtet,
         labelup=rupg, labeldown=rdowng,labelmid=rmidg);        

mesh3 th0=B0+B1+B2+B3+B4+B5+B6+B7+B8+B9+B10+B11;



//savemesh(th0,"holath0.mesh");
//plot(th0,OX(1),OY(1),wait=1,cmm="FIRST ");
mesh3 Th3=movemesh3(th0,transfo=[y,-x,z]);
mesh3 Th3c=movemesh3(Th3,transfo=[x,5+y,z]);
savemesh(Th3,"Th3.mesh");
savemesh(Th3c,"Th3copy.mesh");
//cout<<"GUARDAMOS LA MALLA "<<endl;
//savemesh(th,"holath.mesh");

plot(Th3,OX(1),OY(1),wait=1,
cmm="COMPUTATIONAL DOMAIN AND COPY. POINTS PER BLOCK ..."+nn);
//medit("ThermalFin  ",Th3);

//plot(Th3,Th3c,OX(1),OY(1),wait=1,
//cmm="COMPUTATIONAL DOMAIN AND COPY. POINTS PER BLOCK ..."+nn);
//medit("ThermalFin  ",Th3);

//////////////////////////////////////////////  

//mesh3 th1=B0+B1+B2+B3;
//mesh3 Th3=movemesh3(th1,transfo=[y,-x,z]);
//cout<<"triangulacion terminada! "<<endl;