//load "medit"

//load "gmsh"
//load "iovtk"

int n= 10;
int nvb = (n+1)^3 - (n-1)^3;// Nb boundary vertices
int ntb = n*n*12; // Nb of Boundary triangle 
mesh3 Th=cube(n,n,n);
Th = buildBdMesh(Th); // build the surface mesh


int[int] ref2=[1,0,2,0];
mesh3 Thm=change(Th,label=ref2);


meshS Th1 = Th.Gamma;

assert(Th1.nv==602); 
assert(Th1.nt==1200);
assert(Th1.nbe==120);

plot(Th1);
savemesh(Th1, "Th1.mesh");

meshS Th2 = readmeshS("Th1.mesh");

assert(Th1.nv==Th2.nv); 
assert(Th1.nt==Th2.nt);
assert(Th1.nbe==Th2.nbe);
