// file tetgencube.edp
// this is  a BUGGUS  file work in progress F. Hecht .... Oct 2011
load "tetgen"
load "medit"

real x0,x1,y0,y1;
x0=1.; x1=2.; y0=0.; y1=2*pi;
mesh Thsq1 = square(5,35,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

func ZZ1min = 0;
func ZZ1max = 1.5;
func XX1 = x;
func YY1 = y;

int[int] ref31h = [0,12];
int[int] ref31b = [0,11];

meshS Th31h = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1max],region=ref31h,orientation=1);
meshS Th31b = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1min],region=ref31b,orientation=-1);   

//medit("haut",Th31h);
//medit("bas",Th31b);

/////////////////////////////////
x0=1.; x1=2.; y0=0.; y1=1.5;
mesh Thsq2 = square(5,8,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

func ZZ2 = y;
func XX2 = x;
func YY2min = 0.;
func YY2max = 2*pi;

int[int] ref32h = [0,13];
int[int] ref32b = [0,14];

meshS Th32h = movemesh23(Thsq2,transfo=[XX2,YY2max,ZZ2],region=ref32h,orientation=-1);  
meshS Th32b = movemesh23(Thsq2,transfo=[XX2,YY2min,ZZ2],region=ref32b,orientation=1);

/////////////////////////////////
x0=0.; x1=2*pi; y0=0.; y1=1.5;
mesh Thsq3 = square(35,8,[x0+(x1-x0)*x,y0+(y1-y0)*y]);
func XX3min = 1.;
func XX3max = 2.;

func YY3 = x;
func ZZ3 = y;

int[int] ref33h = [0,15];
int[int] ref33b = [0,16];
cout << "tyty\n";
verbosity = 4; 
meshS Th33h = movemesh23(Thsq3,transfo=[XX3max,YY3,ZZ3],region=ref33h,orientation=1);  
cout << "tyty\n";
meshS Th33b = movemesh23(Thsq3,transfo=[XX3min,YY3,ZZ3],region=ref33b,orientation=-1); 

////////////////////////////////
meshS Th33 = Th31h+Th31b+Th32h+Th32b+Th33h+Th33b; // "gluing" surface meshs to obtain the surface of cube
cout << "premier check" << endl;
checkbemesh(Th33,manifolds=[ [[11,1],[12,1],[13,1],[14,1],[15,1],[16,1]] ]);


func mv2x = x;
func mv2y = y;
func mv2z = -1.5+z;
cout << "movemesh ::new mesh" <<endl;
meshS Th33mv = movemesh(Th33, [mv2x,mv2y,mv2z]);
int[int] r1=[12,11,11,21,13,23,14,24,15,25,16,26];
cout << "change ::new mesh" <<endl;
//medit("Th33mv",Th33mv);
if(0) {
Th33mv=change(Th33mv,label=r1);  // bug change if non tet in mesh .... HERE ...
cout << "add ::new mesh" <<endl;
meshS Thall = Th33+Th33mv;  // The vertex of faces of labels 11 of Th33 and Th33mv must be same. The faces of label 11 of Thall correspond to the faces of the first in the sum : Th33.

medit("Thall",Thall);
// check manifold of the two gluing mesh
cout << "check gluing 1" << endl; 
checkbemesh(Thall,manifolds=[ [[11,1],[12,1],[13,1],[14,1],[15,1],[16,1]] ]);
cout << "check gluing 2" << endl; 
checkbemesh(Thall,manifolds=[ [[11,-1],[21,1],[23,1],[24,1],[25,1],[26,1]] ]);   
cout << "check exterior" << endl;
checkbemesh(Thall,manifolds=[ [[12,1],[13,1],[14,1],[15,1],[16,1],[21,1],[23,1],[24,1],[25,1],[26,1]] ]); // check the exterior manifold 
}
