// Generate a mesh for a box which is periodic in x and y

// Parameters used when reading the mesh later
periodicBoundaryBaseNo=80;
periodicInX=0;
periodicInY=2;
periodicInZ=4;

// pitch
p = 1.0;
// lc mesh (bulk) parameter
lcm = 0.2;
// cell thickness
d = 1.5;

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

Line (1) = {1,2};
Line (2) = {2,3};
Line (3) = {3,4};
Line (4) = {4,1};
Line Loop (1) = {1, 2, 3, 4};
Plane Surface (1) = {1};
out[]=Extrude {0,p,0} {
    Surface{1}; Layers{5};
};

Printf("out[0] = %g", out[0]);
Printf("out[1] = %g", out[1]);
Printf("out[2] = %g", out[2]);
Printf("out[3] = %g", out[3]);
Printf("out[4] = %g", out[4]);
Printf("out[5] = %g", out[5]);

// Delete unwanted entities - leaving the two surfaces periodic in Y and lines 
// in the YoZ plane, which will now be used to form a new Plane Surface to be
// extruded.
Delete{
   Volume{1};
   Surface{13, 17, 21, 25};
   Line{12, 16};
}

// Form the new Plane Surface and extrude it in the X direction.
Line Loop(2) = {4, -20, -9, 11};
Plane Surface(2)={2};
out[]=Extrude {p,0,0} {
    Surface{2}; Layers{1};
};
Printf("out[0] = %g", out[0]);
Printf("out[1] = %g", out[1]);
Printf("out[2] = %g", out[2]);
Printf("out[3] = %g", out[3]);
Printf("out[4] = %g", out[4]);
Printf("out[5] = %g", out[5]);

// Delete the bottom and top ruled surfaces, these will be re-defined as Plane
// Surfaces so that they are meshed at the same resolution as the others. 
Delete{
   Volume{1};
   Surface{39, 47};
}

// Form the new bottom and top surfaces, in a real application these will be
// more complicated.
Line Loop(3) = {-3, -31, 8, -20};
Plane Surface(3) = {3};
Line Loop(4) = {6, -29, -1, 11};
Plane Surface(4) = {4};

// Now define the volume
Surface Loop(100) = {2, 26, 3, 1, 4, 48};
Volume(100) = {100};

// Define the physical entities
pbNo=periodicBoundaryBaseNo+periodicInX;
Physical Surface(pbNo) = {2};    // left periodic
Physical Surface(pbNo+1) = {48}; // right periodic
pbNo=periodicBoundaryBaseNo+periodicInY;
Physical Surface(pbNo) = {1};    // front periodic
Physical Surface(pbNo+1) = {26}; // back periodic

Physical Surface(0)={4};	// Bottom surface
Physical Surface(1)={3};	// Top surface

Physical Volume(0)={100};	// Single region

// Use the frontal algo to make sure that the periodic surface mesh is
// not modified when meshing the volume
Mesh.Algorithm3D = 4;

