
// creates a cube
// defines all lines, surfaces and volume as transfinite alternate, check result to see if there are surface elements not conciding with volume elements

// extrusion does not extend the physical definitions, so we have to do it by hand in a semi-automated way

Include "gmsh-3d-structured-transfinite-alternate.geo.opt.back" ;

Point(1) = {0,0,0};
Point(2) = {1,0,0};
Point(3) = {1,1,0};
Point(4) = {0,1,0};

Point(5) = {0,0,0.5};
Point(6) = {1,0,0.5};
Point(7) = {1,1,0.5};
Point(8) = {0,1,0.5};

Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {4,3};
Line(4) = {1,4};

Line(5) = {5,6};
Line(6) = {6,7};
Line(7) = {8,7};
Line(8) = {5,8};

Line(9) = {1,5};
Line(10) = {2,6};
Line(11) = {3,7};
Line(12) = {4,8};

Geometry.PointNumbers = 0;
Geometry.LineNumbers = 1;
Geometry.Surfaces = 1;

Line Loop(1) = {1,2,-3,-4}; Plane Surface(1) = {1}; // normal of surface pointing inwards
Line Loop(2) = {-5,8,7,-6}; Plane Surface(2) = {2};
Line Loop(3) = {-1,9,5,-10}; Plane Surface(3) = {3};
Line Loop(4) = {3,11,-7,-12}; Plane Surface(4) = {4};
Line Loop(5) = {4,12,-8,-9}; Plane Surface(5) = {5};
Line Loop(6) = {-2,10,6,-11}; Plane Surface(6) = {6};

Geometry.LineNumbers = 0;
Geometry.SurfaceNumbers = 1;

Surface Loop(1) = {1:6};
Volume(1) = {1};

Geometry.SurfaceNumbers = 0;

h_mesh = 0.2 ;

trans_len = 1.0 ;
trans_n = Ceil((trans_len/2.0)/h_mesh)*2+1; // we do this to have an odd number and then transfinite AlternateLeft produces nicer symmetric results on corners
Transfinite Line {1:8} = trans_n ; // all lines

trans_len = 0.5 ;
trans_n = Ceil((trans_len/2.0)/h_mesh)*2+1; // we do this to have an odd number and then transfinite AlternateLeft produces nicer symmetric results on corners
Transfinite Line {9:12} = trans_n ; // all lines

method = 0 ; // works! but interface surface has quads, anyway meshing unstructured works with that
//method = 1; // does not work: some volume faces do not coincide with surface faces
//method = 2; // does not work: some volume faces do not coincide with surface faces
//method = 3; // does not work: internally it is not structured, generates tets and pyramids: quads are mantained which is good and no face elements have all points on boundary
method = 31; // works: internally it is structured, generates tets, hexa and pyramids: quads are mantained which is good and no face elements have all points on boundary
//method = 4; // does not work: some volume faces do not coincide with surface faces
//method = 5; // does not work: some volume faces do not coincide with surface faces
//method = 51; // does not work: some volume faces do not coincide with surface faces
//method = 6; // does not work: some volume faces do not coincide with surface faces
//method = 7; // does not work: some volume faces do not coincide with surface faces
//method = 70; // does not work: some volume faces do not coincide with surface faces
//method = 71; // does not work: some volume faces do not coincide with surface faces
//method = 8; // works! generates tets only, but quads are destroyed to triangles and some elements have three points on boundary which is bad for boundary conditions

If(method==0)
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {1:6} = 10; // not recombining is giving problems of coincidence of volume faces and surface faces
   Transfinite Volume {1};
EndIf

If(method==1)
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Transfinite Volume {1};
EndIf
   
If(method==2)
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf

If(method==31)
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
   TransfQuadTri {1} ;
EndIf

If(method==3)
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   TransfQuadTri {1} ;
EndIf

If(method==4)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 0; Mesh.Recombine3DConformity = 0;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf

If(method==5)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 1; Mesh.Recombine3DConformity = 0;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf

If(method==51)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 1; Mesh.Recombine3DConformity = 1;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf

If(method==6)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 1; Mesh.Recombine3DConformity = 3;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf

If(method==70)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 2; Mesh.Recombine3DConformity = 0;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf

If(method==71)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 2; Mesh.Recombine3DConformity = 1;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf

If(method==7)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 2; Mesh.Recombine3DConformity = 3;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   Transfinite Volume {1};
EndIf



If(method==8)
   Mesh.Recombine3DAll = 1; Mesh.Recombine3DLevel = 1; Mesh.Recombine3DConformity = 0;
   Transfinite Surface {1:6} AlternateLeft; // AlternateLeft + odd number produces bi-symmetric mesh in 2D
   Recombine Surface {3:6} = 10;
   TransfQuadTri {1} ;
EndIf

Mesh 3;

Geometry.Points = 0;
Geometry.Lines = 0;
Geometry.Surfaces = 0;

Printf("**** Finished Meshing");
Printf("NbTetrahedra : %g",Mesh.NbTetrahedra);
Printf("NbNodes      : %g",Mesh.NbNodes);
Printf("NbHexahedra  : %g",Mesh.NbHexahedra);
Printf("NbPrisms     : %g",Mesh.NbPrisms);
Printf("NbPyramids   : %g",Mesh.NbPyramids);
Printf("NbTrihedra   : %g",Mesh.NbTrihedra);

Abort;

//Physical Surface("bottom",1) = {3} ;
//Physical Surface("top",2) = {4} ;
//Physical Surface("zmin",3) = {1} ;
//Physical Surface("zmax",4) = {2} ;
//Physical Surface("xmin",5) = {5} ;
//Physical Surface("xmax",6) = {6} ;

//thevalue = GetNumber(Srt(SensitivityParameter));

// .msh
//$Elements
//9
//1 15 2 0 1 1       one node point
//4 1 2 0 101 1 4    two node line

