load "medit"
load "ioply"

// Parameters
int nn = 10; // Mesh quality
//example of operator for meshS 

// read the imput mesh, vertice/tetra/triangle/edge 
mesh3 Th3 = cube(nn,nn,nn);
Th3 = buildBdMesh(Th3); // build the surface mesh 
meshS ThS = Th3.Gamma;
//split the extrating surface
//plot(Th3);
meshS ThSsplit = trunc(ThS, 1, split=3);
mesh3 Th3split = trunc(Th3, 1, split=2);
savemesh(Th3split, "Th3split.mesh");

meshS Th3splitG = Th3split.Gamma;

medit("buildFromMesh3Splited",Th3splitG);
medit("buildFromMesh3Initial",ThSsplit);

mesh3 Th3mv = movemesh(Th3,[x+1,y/2,z+y]);
Th3mv=buildBdMesh(Th3mv);
meshS Th3mvG = Th3mv.Gamma;
meshS ThSmv = movemesh(ThS,[x+1,y/2,z+y]);

medit("Th3mvG",Th3mvG);
medit("ThSmv",ThSmv); 

meshS ThS2=plyloadS("plane.ply");
plot(ThS2);
//  version 4.11-1.. 13 sep 2022
meshS Ths = square3(10,10,[x,y,square(2*x-1)+square(2*y-1)]); 
plot(Ths,wait=1);
real[int] gzz;
mesh Th2 = movemesh(Ths,transfo=[x,y,z],getZ=gzz);//  get flat 2d mesh ...
fespace Vh2(Th2,P1);
Vh2 gz; gz[]=gzz;
plot(gz,wait=1);

