tetgencube.edp

// file tetgencube.edp 
load "msh3" 
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]; 

mesh3 Th31h = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1max],label=ref31h,orientation=1); 
mesh3 Th31b = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1min],label=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]; 

mesh3 Th32h = movemesh23(Thsq2,transfo=[XX2,YY2max,ZZ2],label=ref32h,orientation=-1); 
mesh3 Th32b = movemesh23(Thsq2,transfo=[XX2,YY2min,ZZ2],label=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]; 

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

//////////////////////////////// 
mesh3 Th33 = Th31h+Th31b+Th32h+Th32b+Th33h+Th33b; // "gluing" surface meshs to obtain the surface of cube 
//medit("glumesh",Th33); 
savemesh(Th33,"Th33.mesh"); 

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


// build a mesh of a half 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]); 
//savemesh(Thmv2,"halfcylindricalshell.mesh"); 
//verbosity=2; 
mesh3 Thmv2surf = movemesh3(Th33, transfo=[mv2x,mv2y,mv2z], facemerge=0); 
medit("maillagesurf",Thmv2surf,wait=1); 
//savemesh(Thmv2surf,"maillagesurfacecylindre.mesh"); 
//medit("maillageplein",Thmv2);