// grid

OMesh = 1;

radiusZone1 = DefineNumber[ 5, Name "Foam/Geometry/radiusMesh", 
			   Label "Mesh radius [m]",Min 4, Max 10, Step 2];
nbElementsTangentialZone1 = DefineNumber[20, Name "Foam/Mesh/nbElementsT",
					 Label "Nb elements tang [-]"];
nbElementsRadialZone1 = DefineNumber[20, Name "Foam/Mesh/nbElementsR", 
				     Label "Nb elements rad [-]"];
progressionRadialZone1 = DefineNumber[1.05, Name "Foam/Mesh/progR", 
				      Label "Progression rad [-]"];

radiusBody = 0.5;
lengthZone2 = 2*2*radiusBody;
//nbElementsTangentialZone1 = 50;
//nbElementsRadialZone1 = 50; // has to be even
//progressionRadialZone1 = 1.05;
nbElementsLengthZone2 = 16;
progressionLengthZone2 = 1.01;
thirdDim = 1;

// body : definition
Point(1) = {0, 0, 0};
Point(11) = {0, -radiusBody, 0};
Point(12) = {-radiusBody, 0, 0};
Point(13) = {0, radiusBody, 0};
Point(14) = {Sqrt(2)/2*radiusBody, Sqrt(2)/2*radiusBody, 0};
Point(15) = {radiusBody, 0, 0};
Point(16) = {Sqrt(2)/2*radiusBody, -Sqrt(2)/2*radiusBody, 0};

Circle(11) = {11, 1, 12};
Circle(12) = {12, 1, 13};
Circle(13) = {13, 1, 14};
Circle(14) = {14, 1, 15};
Circle(15) = {15, 1, 16};
Circle(16) = {16, 1, 11};

// zone1 : definition
Point(101) = {0, -radiusZone1, 0};
Point(102) = {-radiusZone1, 0, 0};
Point(103) = {0, radiusZone1, 0};
Point(104) = {Sqrt(2)/2*radiusZone1, Sqrt(2)/2*radiusZone1, 0};
Point(105) = {radiusZone1, 0, 0};
Point(106) = {Sqrt(2)/2*radiusZone1, -Sqrt(2)/2*radiusZone1, 0};


Circle(101) = {101, 1, 102};
Circle(102) = {102, 1, 103};
Circle(103) = {103, 1, 104};
Circle(104) = {104, 1, 105};
Circle(105) = {105, 1, 106};
Circle(106) = {106, 1, 101};

Line(111) = {11,101};
Line(112) = {12,102};
Line(113) = {13,103};
Line(114) = {14,104};
Line(115) = {15,105};
Line(116) = {16,106};

Line Loop(101) = {101,-112,-11,111};
Plane Surface(101) = {101};
Line Loop(102) = {102,-113,-12,112};
Plane Surface(102) = {102};
Line Loop(103) = {103,-114,-13,113};
Plane Surface(103) = {103};
Line Loop(104) = {104,-115,-14,114};
Plane Surface(104) = {104};
Line Loop(105) = {105,-116,-15,115};
Plane Surface(105) = {105};
Line Loop(106) = {106,-111,-16,116};
Plane Surface(106) = {106};


If  (OMesh == 0)
    // zone2 : definition
    Point(201) = {radiusZone1+lengthZone2, Sqrt(2)/2*radiusZone1, 0};
    Point(202) = {radiusZone1+lengthZone2, 0, 0};
    Point(203) = {radiusZone1+lengthZone2, -Sqrt(2)/2*radiusZone1, 0};

    Line(201) = {104, 201};
    Line(202) = {105, 202};
    Line(203) = {106, 203};
    Line(204) = {201, 202};
    Line(205) = {202, 203};

    Line Loop(201) = {201,204,-202,-104};
    Plane Surface(201) = {201};
    Line Loop(202) = {202,205,-203,-105};
    Plane Surface(202) = {202};
EndIf

// zone1 : meshing
Transfinite Line {11} = nbElementsTangentialZone1+1 Using Progression 1;
Transfinite Line {101} = nbElementsTangentialZone1+1 Using Progression 1;
Transfinite Line {12} = nbElementsTangentialZone1+1 Using Progression 1;
Transfinite Line {102} = nbElementsTangentialZone1+1 Using Progression 1;
Transfinite Line {13} = nbElementsTangentialZone1/2+1 Using Progression 1;
Transfinite Line {103} = nbElementsTangentialZone1/2+1 Using Progression 1;
Transfinite Line {14} = nbElementsTangentialZone1/2+1 Using Progression 1;
Transfinite Line {104} = nbElementsTangentialZone1/2+1 Using Progression 1;
Transfinite Line {15} = nbElementsTangentialZone1/2+1 Using Progression 1;
Transfinite Line {105} = nbElementsTangentialZone1/2+1 Using Progression 1;
Transfinite Line {16} = nbElementsTangentialZone1/2+1 Using Progression 1;
Transfinite Line {106} = nbElementsTangentialZone1/2+1 Using Progression 1;

Transfinite Line {111} = nbElementsRadialZone1+1 Using Progression progressionRadialZone1;
Transfinite Line {112} = nbElementsRadialZone1+1 Using Progression progressionRadialZone1;
Transfinite Line {113} = nbElementsRadialZone1+1 Using Progression progressionRadialZone1;
Transfinite Line {114} = nbElementsRadialZone1+1 Using Progression progressionRadialZone1;
Transfinite Line {115} = nbElementsRadialZone1+1 Using Progression progressionRadialZone1;
Transfinite Line {116} = nbElementsRadialZone1+1 Using Progression progressionRadialZone1;

Transfinite Surface{101} = {};
Recombine Surface {101};
Transfinite Surface{102} = {};
Recombine Surface {102};
Transfinite Surface{103} = {};
Recombine Surface {103};
Transfinite Surface{104} = {};
Recombine Surface {104};
Transfinite Surface{105} = {};
Recombine Surface {105};
Transfinite Surface{106} = {};
Recombine Surface {106};


If  (OMesh == 0)
    // zone2 : meshing
    Transfinite Line {201} = nbElementsLengthZone2+1 Using Progression progressionLengthZone2;
    Transfinite Line {202} = nbElementsLengthZone2+1 Using Progression progressionLengthZone2;
    Transfinite Line {203} = nbElementsLengthZone2+1 Using Progression progressionLengthZone2;
    Transfinite Line {204} = nbElementsTangentialZone1/2+1 Using Progression 1;
    Transfinite Line {205} = nbElementsTangentialZone1/2+1 Using Progression 1;
    
    Transfinite Surface{201} = {};
    Recombine Surface {201};
    Transfinite Surface{202} = {};
    Recombine Surface {202};

    Extrude {0, 0, thirdDim} {Surface{101, 102, 103, 104, 105, 106, 201, 202};Layers{1};Recombine;}

    // patches boundaries : definition
    Physical Surface("frontAndBack") = {101, 227, 102, 249, 103, 271, 104, 293, 315, 105, 337, 106, 359, 201, 381, 202};
    Physical Surface("bodyBoundary") = {222, 244, 266, 288, 310, 332};
    Physical Surface("upStreamBoundary") = {214, 236};
    Physical Surface("downStreamBoundary") = {258, 346, 350, 372, 376, 324};
    Physical Volume("airFlow") = {1, 2, 3, 4, 5, 6, 7, 8};
EndIf

If  (OMesh == 1)
    Extrude {0, 0, thirdDim} {Surface{101, 102, 103, 104, 105, 106};Layers{1};Recombine;}

    // patches boundaries : definition
    Physical Surface("frontAndBack") = {101, 102, 103, 104, 105, 106, 138, 160, 182, 204, 226, 248};
    Physical Surface("bodyBoundary") = {133, 155, 177, 199, 221, 243};
    Physical Surface("upStreamBoundary") = {125, 147};
    Physical Surface("downStreamBoundary") = {169, 191, 213, 235};
    Physical Volume("airFlow") = {1, 2, 3, 4, 5, 6};
EndIf
