load "tetgen"

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

// Mesh
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)
	);

plot(Th, cmm="Th");
plot(Ths, cmm="Ths");

// 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, labelup=r3, labeldown=r4);
savemesh(Th3, "box2region1hole.mesh");
plot(Th3, cmm="Th3");

// 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 << "test=" << test << endl;
mesh3 Th3sph = tetgtransfo(Ths, transfo=[XX1, YY1, ZZ1], switch=test,
	nbofregions=1, regionlist=domain);
savemesh(Th3sph,"sphere2region.mesh");
plot(Th3sph, cmm="Th3sph");

