// Variables
MeshSize =  0.5; // avg meshsize when not explicitly overridden
R = 0.5; // Radius of circle

// Domain lengths
Lx = 20;
Ly = 20;
Lz = 6;

// Number of elements in each direction
Nelx = 30;
Nely = 30;
Nelz = 10;

// The four corners of the domain
p01 = newp ; Point(p01) = { 0.0 , 0.0 , 0.0 , MeshSize};
p02 = newp ; Point(p02) = {  Lx , 0.0 , 0.0 , MeshSize};
p03 = newp ; Point(p03) = {  Lx , Ly , 0.0 , MeshSize};
p04 = newp ; Point(p04) = {  0.0 , Ly , 0.0 , MeshSize};

// Circle center
p05 = newp ; Point(p05) = {  Lx/2.0 , Ly/2.0 , 0.0 , MeshSize};

// Circle end points
p06 = newp ; Point(p06) = {  Lx/2.0-R , Ly/2.0 , 0.0 , MeshSize};
p07 = newp ; Point(p07) = {  Lx/2.0 , Ly/2.0-R , 0.0 , MeshSize};
p08 = newp ; Point(p08) = {  Lx/2.0+R , Ly/2.0 , 0.0 , MeshSize};
p09 = newp ; Point(p09) = {  Lx/2.0 , Ly/2.0+R , 0.0 , MeshSize};

// All the lines of the domain
Line(1)   = {p01,p02};
Line(2)   = {p02,p03};
Line(3)   = {p03,p04};
Line(4)   = {p04,p01};
Circle(5) = {p06, p05, p07};
Circle(6) = {p07, p05, p08};
Circle(7) = {p08, p05, p09};
Circle(8) = {p06, p05, p09};

// Generate equally spaced nodes on all the lines
Transfinite Line(1) = Nelx Using Progression 1 ;
Transfinite Line(2) = Nely Using Progression 1 ;
Transfinite Line(3) = Nelx Using Progression 1 ;
Transfinite Line(4) = Nely Using Progression 1 ;
Transfinite Line(5) = 10 Using Progression 1 ;
Transfinite Line(6) = 10 Using Progression 1 ;
Transfinite Line(7) = 10 Using Progression 1 ;
Transfinite Line(8) = 10 Using Progression 1 ;

// Line loops and create the surface with a hole in it
Line Loop(1) = {1, 2, 3, 4};
Line Loop(2) = {5, 6, 7, -8};
Plane Surface(1) = { 1, 2 };
Plane Surface(2) = { 2 };

//Transfinite Surface(1) = {1, 2};
Mesh.Algorithm = 8;
Mesh.SubdivisionAlgorithm = 1; // To refine the mesh, otherwise comment. // Forces it to quads // Mesh subdivision algorithm (0=none, 1=all quadrangles, 2=all hexahedra)
Recombine Surface(1) = 20;
Recombine Surface(2) = 20;

v[] = Extrude {0,0,Lz} { Surface{1,2}; Layers{2}; Recombine;};
Physical Volume (1) = {v[1]};
