// example40.edp from file tetgencube.edp
verbosity=2;
load "msh3"
load "tetgen"
load "medit"

// front (z=0) and back (z=1.5) faces
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];

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

// bottom (y=0) and top (y=2*pi) faces
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];

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

// Left (x=1) and right (x=2) faces
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];

mesh3 Th33h = movemesh23(Thsq3, transfo = [XX3max, YY3, ZZ3],
                label=ref33h, orientation=1);  
mesh3 Th33b = movemesh23(Thsq3, transfo = [XX3min, YY3, ZZ3],
                label=ref33b, orientation=-1); 

// glue surfaces together to make surface of rectangular parallelopiped
mesh3 Th33 = Th31h + Th31b + Th32h + Th32b + Th33h + Th33b;

medit("glumesh",Th33);  // plot using medit

plot(Th33);             // plot using FreeFem++

// savemesh(Th33,"Th33.mesh");

// build a mesh of a axis parallel box with TetGen
real[int] domaine = [1.5,pi,0.75,145,0.0025];
// Tetrahelize the interior of the cube with tetgen
mesh3 Thfinal = tetg(Th33,switch="pqaAAYYQ",nbofregions=1,regionlist=domaine);
medit("tetg",Thfinal);
//savemesh(Thfinal,"Thfinal.mesh");


// build a mesh of a cylindrical shell of interior radius 1. 
// and exterior radius 2 and heigh 1.5
func mv2x = x*cos(y);
func mv2y = x*sin(y);
func mv2z = z;
mesh3 Thmv2 = movemesh3(Thfinal, transfo=[mv2x, mv2y, mv2z]);
medit("cylindricalshell",Thmv2);
plot(Thmv2);
//savemesh(Thmv2,"cylindricalshell.mesh");